## Square Digit Chains

### October 19, 2018

Today’s task appears on Project Euler, Rosetta Code, and Sloane’s, so it is well worth a look; we even did a version of this task in a previous exercise. Here is the Project Euler version of the problem:

A number chain is created by continuously added the square of the digits of a number to form a new number until it has been seen before. For example:

44 → 32 → 13 → 10 → 1 → 1

85 → 89 → 145 → 42 → 20 → 4 → 16 → 37 → 58 → 89

Therefore any chain that arrives at 1 or 89 will become stuck in an endless loop. What is most amazing is that EVERY starting number will eventually arrive at 1 or 89.

How many starting numbers below ten million will arrive at 89?

Your task is to write a program to find the count of starting numbers below ten million that arrive at 89, following the usual Project Euler rule that your computation may take no more than a minute of computation on a recent-vintage personal computer. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

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([1]), set([89])
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[7] // reduce(mul, (FACS[f] for f in cnt.values()), 1)
```
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.

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
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] = 1;
lookup[89] = 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[0] = # zeros, a[1] = # 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[0] = 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[0] *= facts[a[0]]
return True
return False

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

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