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.
Programming Praxis 
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) # -> 8581146The 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) return totalAt 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
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; }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. :)
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)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
===
Run:
s secs=$p($h,”,”,2) x loop w !!,”sum = “,sum,!!,”Seconds = “,$p($h,”,”,2)-secs
iterations: 9999999
sum = 4006299
Seconds = 101
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:
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]; 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:
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