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.
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