## Square Digit Chains

### October 19, 2018

We begin with a simple program to display one chain:

```(define (sdc n)
(cond ((member n '(1 89)) n)
(else (display n) (display #\space)
(sdc (sum (map square (digits n)))))))

> (sdc 44)
44 32 13 10 1
> (sdc 145)
145 42 20 4 16 37 58 89```

We could invoke that program ten million times, but it wouldn’t finish in any reasonable time. Instead, we observe that 1234567 and 7654321, and indeed any permutation of those seven digits, will lead to the same square digit chain, so we only have to compute the canonical form of all the permutations less than ten million, compute the number of permutations for each canonical form that leads to 89, and add them up:

```(define (euler92)
(define (sdc n)
(if (member n '(0 1 89)) n
(sdc (sum (map square (digits n))))))
(define (perms xs)
(define (fact n) (apply * (range n 0)))
(let loop ((p (fact (length xs))) (freq (uniq-c = xs)))
(if (null? freq) p
(loop (/ p (fact (cdar freq))) (cdr freq)))))
(fold-of + 0 p
(a range 0 10)
(b range a 10)
(c range b 10)
(d range c 10)
(e range d 10)
(f range e 10)
(g range f 10)
(n is (list a b c d e f g))
(= (sdc (undigits n)) 89)
(p is (perms n))))

> (time (euler92))
(time (euler92))
2 collections
0.038801653s elapsed cpu time, including 0.000115820s collecting
0.038834982s elapsed real time, including 0.000122825s collecting
20564512 bytes allocated, including 16846320 bytes reclaimed
8581146```

There are 11,440 canonical forms, of which 9814 lead to 89. You can run the program at https://ideone.com/w1QtD8.

Pages: 1 2

### 11 Responses to “Square Digit Chains”

1. Paul said

In Python.

```SQUARES = {str(i): i**2 for i in range(10)}

def program(limit=10_000_000):
S1, S89 = set(), set()
for i in range(2, limit):
j = i
js = [j]
while j not in S1 and j not in S89:
j = sum(SQUARES[nstr] for nstr in str(j))
js.append(j)
if j in S89:
S89 |= set(js)
else:
S1 |= set(js)
return len(S89)    # -> 8581146
```
2. Paul said

The solution of ProgrammingPraxis in Python. This is indeed a lot faster than my solution.

```import collections
from functools import reduce
from itertools import combinations_with_replacement
from math import factorial
from operator import mul

SQUARES = {str(i): i**2 for i in range(10)}
FACS = list(map(factorial, range(8)))

def program1(limit=10_000_000):
total = 0
for nlist in combinations_with_replacement("0123456789", 7):
j = int("".join(nlist))
while j not in (0, 1, 89):
j = sum(SQUARES[nstr] for nstr in str(j))
if j == 89:
cnt = collections.Counter(nlist)
total += FACS // reduce(mul, (FACS[f] for f in cnt.values()), 1)
return total
```
3. Ernie said

At most we have 7 digits. The range of the sumsq of the digits will be 1…567 (7.9^2). Thus every
sumsq yields a value in [1, 567]. So first we calculate the chains for 1…567 storing the result
in table[]. Thereafter we calculate sumsq(n), n > 567, and increment count if table[sumsq[n]) = 89.

Answer 8581146 0.35 sec

4. Noddysevens said

A brute force solution in cpp

```#include <iostream>

int get_square_sum(int a);
const int MAX_ROUNDS = 10000000;

int main()
{
int count = 0;
for (int i = 1; i < MAX_ROUNDS; i++)
{
int a = i;
while(a != 0 && a != 1 && a != 89)
{
a = get_square_sum(a); if (a == 89) count++;
}
}
std::cout << "Found " << count << " out of " << MAX_ROUNDS << " starting numbers arriving at 89\n";
return 0;
}

int get_square_sum(int a)
{
int factor = 1, sum = 0;
while((a / factor) > 0)
{
int digit = a / factor % 10;
sum += digit * digit;
factor *= 10;
}
return sum;
}
```
5. Russell Glasser said

I’d just like to mention that, although it’s of course admirable to simplify by developing a permutations based solution, doing a full check on 10,000,000 numbers took about half a second on my current desktop. Maybe I’m not clear on what a “recent-vintage personal computer” is meant to encompass, but it doesn’t need to be optimized if it meets the requirements of the problem. :)

6. matthew said

Here’s a nice way to solve the problem (I took the basic idea from the first C solution on the Rosetta code page): define a function count(n,m) that gives the number of n digit numbers (including leading zeros) that have a square digit sum of m – this has a fairly simple recursive formulation that by itself is rather inefficient but can be improved dramatically by memoizing its arguments (alternatively we can use a dynamic programming approach, see the above mentioned C solution for details). With Python’s bignums we can easily count all the 20 digit (or more) numbers with iterated square sum of 89:

```from functools import lru_cache

@lru_cache(maxsize=None)
def sd(n):
while n not in [1,89]: n = sum(int(i)**2 for i in str(n))
return n

@lru_cache(maxsize=None)
def count(n,m):
if m < 0: return 0
if n == 0: return int(m == 0)
return sum(count(n-1,m-i*i) for i in range(10))

for n in range(1,21):
total = sum(count(n,i) for i in range(1,n*81+1) if sd(i) == 89)
print(n,total)
```
```\$ time python3 sdsum.py
1 7
2 80
3 857
4 8558
5 85623
6 856929
7 8581146
8 85744333
9 854325192
10 8507390852
11 84908800643
12 850878696414
13 8556721999130
14 86229146720315
15 869339034137667
16 8754780882739336
17 87975303595231975
18 881773944919974509
19 8816770037940618762
20 87994965555707002706

real	0m0.408s
user	0m0.348s
sys	0m0.008s
```
7. Klong and Mumps solutions

```
Notes:
1) The times for both solutions are quite long, due in part at least by running under a VM.
2) Both of the solutions predicted 4,006,299 chains ending in 89.  Not sure why ProgrammingPraxis' solution indicated another number.

---
Klong version
:" Converts number into list of digits"
cnv::{{1:\$1\$x}'\$x}

:" Sums the squares of a list of digits"
sumsq::{+/x^2}

:" Given a starting number, return a chain"
mkch::{[last lst]; lst::,x; {x; 1=#lst?last::*|lst}{lst::lst,{+/x^2}({{1:\$1\$x}'\$x}(last)); x+1}:~1; lst}

:" Analyze x chains"
chsumsq::{[limit lst sum]; limit::x; sum::0; limit{[last]; lst::mkch(x); :[89=*|lst; sum::sum+1; sum::sum]; x+1}:*1; .p("Starting numbers arriving at 89: ",\$sum)}
```

Run:
Klong 20180213
]l lib/steve2.kg
loading “./lib/steve2.kg”
num::9999999; .sys(“date”); chsumsq(num); .sys(“date”)
Tue Oct 23 15:03:11 PDT 2018
Starting numbers arriving at 89: 4006299
Tue Oct 23 15:35:58 PDT 2018
0

===

```
Mumps version
s cnv="s sumsq=0 f  s sumsq=sumsq+(\$e(num)*\$e(num)),num=\$e(num,2,\$l(num)) q:'\$l(num)"
s gench="k lst s lst=num,lst(num)=1 f  x cnv s num=sumsq,lst=lst_"" ""_num,lst(num)=\$g(lst(num))+1 i lst(num)=2 s:num=89 sum=sum+1 q"
s loop="r !!,""# iterations: "",iter s sum=0 f i=1:1:iter s num=i x gench"
```

Run:
s secs=\$p(\$h,”,”,2) x loop w !!,”sum = “,sum,!!,”Seconds = “,\$p(\$h,”,”,2)-secs

iterations: 9999999

sum = 4006299

Seconds = 101

8. Daniel said

Here’s a solution in C. On my desktop it takes about 0.30 seconds without compiler optimizations and 0.10 seconds with -O2.

```#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

#define MAX_INPUT 9999999
#define MAX_DIGIT_SQ_SUM 567

int digit_sq_sum(int x) {
int result = 0;
while (x != 0) {
int mod = x % 10;
result += mod * mod;
x /= 10;
}
return result;
}

int main(int argc, char* argv[]) {
uint8_t* lookup = calloc(MAX_DIGIT_SQ_SUM + 1, sizeof(uint8_t));
lookup = 1;
lookup = 89;
for (int x = 1; x <= MAX_DIGIT_SQ_SUM; ++x) {
int y = x;
int* path = NULL;
int n_path = 0;
while (true) {
uint8_t state = lookup[y];
if (state != 0) {
for (int i = 0; i < n_path; ++i) {
lookup[path[i]] = state;
}
break;
}
++n_path;
path = realloc(path, n_path * sizeof(int));
path[n_path - 1] = y;
y = digit_sq_sum(y);
}
free(path);
}
int count = 0;
for (int x = 1; x <= MAX_INPUT; ++x) {
int y = x;
if (y > MAX_DIGIT_SQ_SUM) y = digit_sq_sum(y);
count += lookup[y] == 89;
}
printf("%d\n", count);
free(lookup);
return EXIT_SUCCESS;
}
```

Output:

```8581146
```
9. matthew said

Here’s a nice variation on the combinations method: generate the combinations explicitly, but at the same time keep track of the digit sum and factorials, adjusting them as we go:

```from math import factorial

def sd(n):
while True:
if n in (0,1,89): return n
n = sum(int(d)*int(d) for d in str(n))

facts = [factorial(i) for i in range(30)]
sds = [sd(i) for i in range(2000)]

def nextcombo(N,a,b,f):
"""
Return the next combination in a
N is number of digits
a is list of digit counts, a = # zeros, a = # ones, etc.
b is cumulative (from the right) square digit sum
f is cumulative (from the right) product of factorials of a
"""
for i in range(N-1):
# Find a non-zero entry
if a[i] > 0:
# Move 1 to right and the rest to the beginning
# eg. [0,0,3,2,...] -> [2,0,0,3,...]
n = a[i]
a[i] = 0
a = n-1;
a[i+1] += 1
# Adjust digit sum and factorials
b[i+1] += (i+1)*(i+1)
f[i+1] *= a[i+1]
for j in range(i+1):
b[j] = b[i+1]
f[j] = f[i+1]
f *= facts[a]
return True
return False

def check(N,M):
a,b,f = *N,*N,*N
a = M; f = facts[a]
total = 0
while nextcombo(N,a,b,f):
# b is square digit sum
# f is product of factorials
if sds[b] == 89: total += facts[M]//f;
return total

for i in range(1,21): print(i,check(10,i))
```

Not as fast as previous solution, but gets up to 20 digits in a minute or two:

```\$ time python3 sb-pub.py
1 7
2 80
3 857
4 8558
5 85623
6 856929
7 8581146
8 85744333
9 854325192
10 8507390852
11 84908800643
12 850878696414
13 8556721999130
14 86229146720315
15 869339034137667
16 8754780882739336
17 87975303595231975
18 881773944919974509
19 8816770037940618762
20 87994965555707002706

real	1m42.021s
user	1m42.000s
sys	0m0.016s
```
10. Eric said

Wow. The official solution put in a lot of effort. I just slapped some memoization into a simplistic approach and called it acceptable. The added restriction of one minute compute time seemed a bit silly. https://ideone.com/9xi2dd