Number equal to the sum of powers of its digits

Question

I've got another interesing programming/mathematical problem.

For a given natural number q from interval [2; 10000] find the number n
which is equal to sum of q-th powers of its digits modulo 2^64.

for example: for q=3, n=153; for q=5, n=4150.

I wasn't sure if this problem fits more to math.se or stackoverflow, but this was a programming task which my friend told me quite a long time ago. Now I remembered that and would like to know how such things can be done. How to approach this?

Answer 1

There are two key points,

• the range of possible solutions is bounded,
• any group of numbers whose digits are the same up to permutation con contain at most one solution.

Let us take a closer look at the case q = 2. If a d-digit number n is equal to the sum of the squares of its digits, then

n >= 10^(d-1) // because it's a d-digit number
n <= d*9^2    // because each digit is at most 9

and the condition 10^(d-1) <= d*81 is easily translated to d <= 3 or n < 1000. That's not many numbers to check, a brute-force for those is fast. For q = 3, the condition 10^(d-1) <= d*729 yields d <= 4, still not many numbers to check. We could find smaller bounds by analysing further, for q = 2, the sum of the squares of at most three digits is at most 243, so a solution must be less than 244. The maximal sum of squares of digits in that range is reached for 199: 1² + 9² + 9² = 163, continuing, one can easily find that a solution must be less than 100. (The only solution for q = 2 is 1.) For q = 3, the maximal sum of four cubes of digits is 4*729 = 2916, continuing, we can see that all solutions for q = 3 are less than 1000. But that sort of improvement of the bound is only useful for small exponents due to the modulus requirement. When the sum of the powers of the digits can exceed the modulus, it breaks down. Therefore I stop at finding the maximal possible number of digits.

Now, without the modulus, for the sum of the q-th powers of the digits, the bound would be approximately

q - (q/20) + 1

so for larger q, the range of possible solutions obtained from that is huge.

But two points come to the rescue here, first the modulus, which limits the solution space to 2 <= n < 2^64, at most 20 digits, and second, the permutation-invariance of the (modular) digital power sum.

The permutation invariance means that we only need to construct monotonous sequences of d digits, calculate the sum of the q-th powers and check whether the number thus obtained has the correct digits.

Since the number of monotonous d-digit sequences is comparably small, a brute-force using that becomes feasible. In particular if we ignore digits not contributing to the sum (0 for all exponents, 8 for q >= 22, also 4 for q >= 32, all even digits for q >= 64).

The number of monotonous sequences of length d using s symbols is

binom(s+d-1, d)

s is for us at most 9, d <= 20, summing from d = 1 to d = 20, there are at most 10015004 sequences to consider for each exponent. That's not too much.

Still, doing that for all q under consideration amounts to a long time, but if we take into account that for q >= 64, for all even digits x^q % 2^64 == 0, we need only consider sequences composed of odd digits, and the total number of monotonous sequences of length at most 20 using 5 symbols is binom(20+5,20) - 1 = 53129. Now, that looks good.

Summary

We consider a function f mapping digits to natural numbers and are looking for solutions of the equation

n == (sum [f(d) | d <- digits(n)] `mod` 2^64)

where digits maps n to the list of its digits.

From f, we build a function F from lists of digits to natural numbers,

F(list) = sum [f(d) | d <- list] `mod` 2^64

Then we are looking for fixed points of G = F ∘ digits. Now n is a fixed point of G if and only if digits(n) is a fixed point of H = digits ∘ F. Hence we may equivalently look for fixed points of H.

But F is permutation-invariant, so we can restrict ourselves to sorted lists and consider K = sort ∘ digits ∘ F.

Fixed points of H and of K are in one-to-one correspondence. If list is a fixed point of H, then sort(list) is a fixed point of K, and if sortedList is a fixed point of K, then H(sortedList) is a permutation of sortedList, hence H(H(sortedList)) = H(sortedList), in other words, H(sortedList) is a fixed point of K, and sort resp. H are bijections between the set of fixed points of H and K.

A further improvement is possible if some f(d) are 0 (modulo 2⁶⁴). Let compress be a function that removes digits with f(d) mod 2^64 == 0 from a list of digits and consider the function L = compress ∘ K.

Since F ∘ compress = F, if list is a fixed point of K, then compress(list) is a fixed point of L. Conversely, if clist is a fixed point of L, then K(clist) is a fixed point of K, and compress resp. K are bijections between the sets of fixed points of L resp. K. (And H(clist) is a fixed point of H, and compress ∘ sort resp. H are bijections between the sets of fixed points of L resp. H.)

The space of compressed sorted lists of at most d digits is small enough to brute-force for the functions f under consideration, namely power functions.

So the strategy is:

1. Find the maximal number d of digits to consider (bounded by 20 due to the modulus, smaller for small q).
2. Generate the compressed monotonic sequences of up to d digits.
3. Check whether the sequence is a fixed point of L, if it is, F(sequence) is a fixed point of G, i.e. a solution of the problem.

Code

Fortunately, you haven't specified a language, so I went for the option of simplest code, i.e. Haskell:

{-# LANGUAGE CPP #-}
module Main (main) where

import Data.List
import Data.Array.Unboxed
import Data.Word
import Text.Printf

#include "MachDeps.h"

#if WORD_SIZE_IN_BITS == 64

type UINT64 = Word

#else

type UINT64 = Word64

#endif

maxDigits :: UINT64 -> Int
maxDigits mx = min 20 $ go d0 (10^(d0-1)) start
  where
    d0 = floor (log (fromIntegral mx) / log 10) + 1
    mxi :: Integer
    mxi = fromIntegral mx
    start = mxi * fromIntegral d0
    go d p10 mmx
        | p10 > mmx = d-1
        | otherwise = go (d+1) (p10*10) (mmx+mxi)

sortedDigits :: UINT64 -> [UINT64]
sortedDigits = sort . digs
  where
    digs 0 = []
    digs n = case n `quotRem` 10 of
               (q,r) -> r : digs q

generateSequences :: Int -> [a] -> [[a]]
generateSequences 0 _ 
    = [[]]
generateSequences d [x] 
    = [replicate d x]
generateSequences d (x:xs)
    = [replicate k x ++ tl | k <- [d,d-1 .. 0], tl <- generateSequences (d-k) xs]
generateSequences _ _ = []

fixedPoints :: (UINT64 -> UINT64) -> [UINT64]
fixedPoints digFun = sort . map listNum . filter okSeq $ 
                        [ds | d <- [1 .. mxdigs], ds <- generateSequences d contDigs]
  where
    funArr :: UArray UINT64 UINT64
    funArr = array (0,9) [(i,digFun i) | i <- [0 .. 9]]
    mxval = maximum (elems funArr)
    contDigs = filter ((/= 0) . (funArr !)) [0 .. 9]
    mxdigs = maxDigits mxval
    listNum = sum . map (funArr !)
    numFun = listNum . sortedDigits
    listFun = inter . sortedDigits . listNum
    inter = go contDigs
      where
        go cds@(c:cs) dds@(d:ds)
            | c < d     = go cs dds
            | c == d    = c : go cds ds
            | otherwise = go cds ds
        go _ _ = []
    okSeq ds = ds == listFun ds


solve :: Int -> IO ()
solve q = do
    printf "%d:\n    " q
    print (fixedPoints (^q))

main :: IO ()
main = mapM_ solve [2 .. 10000]

It's not optimised, but as is, it finds all solutions for 2 <= q <= 10000 in a little below 50 minutes on my box, starting with

2:
    [1]
3:
    [1,153,370,371,407]
4:
    [1,1634,8208,9474]
5:
    [1,4150,4151,54748,92727,93084,194979]
6:
    [1,548834]
7:
    [1,1741725,4210818,9800817,9926315,14459929]
8:
    [1,24678050,24678051,88593477]
9:
    [1,146511208,472335975,534494836,912985153]
10:
    [1,4679307774]
11:
    [1,32164049650,32164049651,40028394225,42678290603,44708635679,49388550606,82693916578,94204591914]

And ending with

9990:
    [1,12937422361297403387,15382453639294074274]
9991:
    [1,16950879977792502812]
9992:
    [1,2034101383512968938]
9993:
    [1]
9994:
    [1,9204092726570951194,10131851145684339988]
9995:
    [1]
9996:
    [1,10606560191089577674,17895866689572679819]
9997:
    [1,8809232686506786849]
9998:
    [1]
9999:
    [1]
10000:
    [1,11792005616768216715]

The exponents from about 10 to 63 take longest (individually, not cumulative), there's a remarkable speedup from exponent 64 on due to the reduced search space.