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.

Advertisements

Pages: 1 2

10 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)
        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] = 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];
        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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: