## Curious Numbers

### June 28, 2016

As in many of the Project Euler problems, the solution is easy if you know the trick and impossible if you don’t. A common pattern in solving these problems is to find the beginning of the solution sequence, then look up the sequence in the Online Encyclopedia of Integer Sequences. Here is a simple demonstration showing the beginning of the sequence of curious numbers:

```> (do ((n 1 (+ n 1))) (#f)
(when (= n (modulo (* n n) (expt 10 (+ (ilog 10 n) 1))))
(display n) (newline)))
1
5
6
25
76
376
625
9376
90625
109376
890625
2890625
7109376
12890625
^C```

That’s enough to plug in to the search engine at OEIS, where we find sequence A003226. Reading the comments there, and referring to A007185 and A016090, we quickly come up with the following function:

```(define (curious n)
(let ((xs (list 1)))
(do ((i 1 (+ i 1))) ((= i n) (apply + xs))
(let* ((x (expm 5 (expt 2 i) (expt 10 i)))
(y (- (expt 10 i) x -1)))
(display i) (display ": ")
(display x) (display " ")
(display y) (newline)
(when (not (member x xs))
(set! xs (cons x xs)))
(when (not (member y xs))
(set! xs (cons y xs)))))))```

And now it is easy to find the requested sum:

```> (time (curious 20))
1: 5 6
2: 25 76
3: 625 376
4: 625 9376
5: 90625 9376
6: 890625 109376
7: 2890625 7109376
8: 12890625 87109376
9: 212890625 787109376
10: 8212890625 1787109376
11: 18212890625 81787109376
12: 918212890625 81787109376
13: 9918212890625 81787109376
14: 59918212890625 40081787109376
15: 259918212890625 740081787109376
16: 6259918212890625 3740081787109376
17: 56259918212890625 43740081787109376
18: 256259918212890625 743740081787109376
19: 2256259918212890625 7743740081787109376
(time (curious 20))
no collections
0.000191847s elapsed cpu time
0.000190437s elapsed real time
10560 bytes allocated
11111110947536882377```

There are duplicate entries for 625, 9376 and 81787109376 (triplicate!) because the leading digit of the i-digit number is zero. Notice that the trailing digits, once calculated, never disappear from the sequence.

We used `ilog` and `expm` from the Standard Prelude. You can run the program at http://ideone.com/nkj3BC. You might also look at this article at John Cook’s blog, which inspired the exercise.

Pages: 1 2 3

### 5 Responses to “Curious Numbers”

1. Daniel said

Here’s a solution in Java. I noticed, as the solution mentions, that “trailing digits, once calculated, never disappear from the sequence”, and used that to generate subsequent numbers.

However, my solution, 103367370865749773002, includes a number that was seemingly omitted from the programmingpraxis solution sum (I haven’t tried figuring out why yet), 92256259918212890625.

```import java.math.BigInteger;

public class PraxisCuriousNumbers {
private static BigInteger getNextCurious(BigInteger integer) {
for (BigInteger i = BigInteger.ONE; true; i = i.add(BigInteger.ONE)) {
BigInteger next = new BigInteger(i.toString()+integer.toString());
if (next.multiply(next).toString().endsWith(next.toString())) {
return next;
}
}
}

public static void main(String[] args) {
BigInteger limit = new BigInteger("10").pow(20);
BigInteger sum = new BigInteger("12"); // 0 + 1 + 5 + 6

for (String start : new String[]{"5", "6"}) {
BigInteger current = new BigInteger(start);
while (true) {
current = getNextCurious(current);
if (current.compareTo(limit) < 0) {
} else {
break;
}
}
}

System.out.println("sum: " + sum);
}
}
```

Here’s the output:

```sum: 103367370865749773002
```
2. Daniel said

Here’s the output showing running time (0.133s) using a 2013 laptop with a 2.8 GHz Intel Core i7:

```\$ time java PraxisCuriousNumbers
sum: 103367370865749773002

real	0m0.133s
user	0m0.125s
sys	0m0.029s
```
3. programmingpraxis said

@Daniel: I mistakenly wrote (curious 20) instead of (curious 21). Your total is correct.

4. al said

We can build up the set of n-digit curious numbers recursively: if c is a curious n-digit number, let

c = x*(10^(n-1)) + y
where
* 10^(n-k-1) < y < 10^(n-k), for some k > 0, and
* x is in [1..9]

Then since
c = y (mod 10^(n-k)), and
c = c*c (mod 10^n)
we have
c = c*c = y*y = y (mod 10^(n – k))
so y is also curious

So we can build up the list l(n) of curious numbers of at most n digits from the list l(n-1) of curious numbers of at most (n-1) digits fairly easily by filtering the curious numbers from the following list of candidates:

candidates(n) = [x*(10^(n-1)) + y | x <- [0..9], y <- l(n-1)]

i.e. l(n) = filter isCurious candidates(n)

(This may generate repeats so we will need to take unique elements of the list.)

Since curious numbers are rare, the list of curious numbers at the previous step is small, so filtering from this candidate set is not terribly expensive.

In Haskell, runs in less than 0.1s on a core i5-6200U processor.

module Main where
import Data.List

isCurious :: Integer -> Bool
isCurious n = n == (mod (n * n) p) where
powersOfTen = map (\x -> 10^x) [1..]
p = head \$ dropWhile (< n) powersOfTen

curiousUpToNDigits :: Integer -> [Integer]
curiousUpToNDigits n
| n == 1 = [1,5,6]
| n > 1 = let l = (curiousUpToNDigits (n-1)) in
nub \$ filter isCurious [x*(10^(n-1)) + y | x <- [0..9], y <- l]
| otherwise = []

— Compute the sum of curious numbers of at most n digits

sumOfCurious :: Integer -> Integer
sumOfCurious n = sum \$ curiousUpToNDigits n

main :: IO ()
main = do
let sum = sumOfCurious 20
print “The sum is:”
print sum

5. matthew said

Nice problem. I like the connection with p-adic numbers.

Here’s another way of directly calculating curious numbers: k-digit A is curious if A*A = A mod 10ᵏ, ie. A*A-A = A*(A-1) = X*10ᵏ = X*5ᵏ*2ᵏ for some X. A and A-1 must be co-prime so 5ᵏ must divide one and 2ᵏ must divide the other, which means that A satisfies:

A mod 5ᵏ = 0 and A mod 2ᵏ = 1

or:

A mod 5ᵏ = 1 and A mod 2ᵏ = 0

and we can solve such modular equations with the Chinese Remainder Theorem: a = x mod p, a = y mod q has solution a = y*p*r + x*q*s where r and s are modular inverses of p and q, ie: p*r = 1 mod q and q*s = 1 mod p.

Here’s a Python program using the extended Euclidean algorithm to calculate the modular inverses:

```def egcd(a,b):
x = 1; y = 0; z = 0; w = 1
while b != 0:
q,r = divmod(a,b)
a,b =  b,r
x,y,z,w = z,w,x-q*z,y-q*w
return a,x,y

def curious(k):
p = 5**k; q = 2**k
(d,r,s) = egcd(p,q)
if r < 0: r += q
if s < 0: s += p
return (p*r,q*s)

def curioussum(n):
sum = 1
for k in range(1,n+1):
min = 10**(k-1)
(a,b) = curious(k)
if a >= min: sum += a
if b >= min: sum += b
return sum

>>> print(curioussum(20))
103367370865749773002
```

Interestingly, for larger k (> 100 or so) this approach is significantly faster than the suggested method – computation of the GCD is faster than exponentiation modulo 10ᵏ – we can compute eg. a 1 million digit curious numbers in a few minutes:

```>>> print("%d\n%d"%curious(1000000))
4770144151139899673799525600941150267423...
5229855848860100326200474399058849732576...
```