algorithm - 数字等于其数字的幂的总和

Question

我还有另一个有趣的编程/数学问题。

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 q=3, n=153; for q=5, n=4150。

我不确定这个问题是否更适合math.se或stackoverflow，但这是我的 friend 很久以前告诉我的一项编程任务。现在，我想起了这一点，并想知道如何完成这些事情。如何处理呢？

Answer 1

有两个关键点，

可能的解决方案范围是有界的，

直到置换为止数字相同的任何一组数字最多包含一个解决方案。

让我们仔细看看 q = 2的情况。如果 d数字位数 n等于其位数的平方和，则

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

并且条件 10^(d-1) <= d*81可以轻松转换为 d <= 3或 n < 1000。要检查的数字并不多，对它们的暴力破解速度很快。对于 q = 3，条件 10^(d-1) <= d*729产生 d <= 4，但仍然没有很多要检查的数字。通过进一步分析，我们可以找到较小的边界，对于 q = 2，最多三个数字的平方和最多为243，因此解决方案必须小于244。在该范围内的最大数字平方和达到199 :1²+9²+9²= 163，继续，您可以轻松地找到一个解必须小于100。( q = 2的唯一解决方案是1。)对于 q = 3，四个数字立方体的最大和为4 * 729 = 2916，继续，我们可以看到 q = 3的所有解都小于1000。但是由于模数要求，这种边界的改进仅对小指数有用。当数字的幂的总和可以超过模数时，它将分解。因此，我停止寻找最大可能的数字位数。
现在，如果没有模数，对于数字的 q次幂的和，边界将近似为

q - (q/20) + 1

因此对于更大的 q，从中获得的可能解决方案的范围是巨大的。
但是这里有两点需要解决，首先是模数，​​它最多将解决方案空间限制为 2 <= n < 2^64，最多20位，其次是(模块化)数字功率和的排列不变性。
排列不变性意味着我们只需要构造 d数字的单调序列，计算 q -th次幂的和，并检查由此获得的数字是否具有正确的数字。
由于单调的 d位数字序列的数量相对较小，因此使用它的蛮力变得可行。特别是如果我们忽略对总和没有贡献的数字(所有指数为0， q >= 22为8， q >= 32也为4， q >= 64的所有偶数)。
使用 d符号的长度为 s的单调序列的数量为

binom(s+d-1, d)

s对于我们来说最多9个， d <= 20，从 d = 1到 d = 20求和，每个指数最多要考虑10015004个序列。还不算太多
尽管如此，对所考虑的所有 q这样做花费的时间很长，但如果考虑到 q >= 64，对于所有偶数位 x^q % 2^64 == 0，我们只需要考虑由奇数位组成的序列以及长度为单调的序列的总数最多使用5个符号的20个字符是 binom(20+5,20) - 1 = 53129。现在，看起来不错。
概括
我们考虑将数字映射到自然数的 f函数，并在寻找方程的解

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

其中 digits将 n映射到其数字列表。
从 f，我们从数字列表到自然数构建一个 F函数，

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

然后，我们在寻找 G = F ∘ digits的固定点。现在，且仅当 n是 G的固定点时， digits(n)才是 H = digits ∘ F的固定点。因此，我们可以等效地寻找 H的固定点。
但是 F是置换不变的，因此我们可以将自己限制在已排序列表中，并考虑 K = sort ∘ digits ∘ F。 H和 K的不动点是一一对应的。如果 list是 H的一个固定点，则 sort(list)是 K的一个固定点，如果 sortedList是 K的一个固定点，则 H(sortedList)是 sortedList的一个置换，因此 H(H(sortedList)) = H(sortedList)(即 H(sortedList))是 K的一个固定点。和 sort分别。 H是 H和 K的不动点集之间的双射。
如果某些 f(d)为0(模264)，则可能会进一步改善。令 compress为一个从数字列表中删除带有 f(d) mod 2^64 == 0的数字的函数，并考虑该函数 L = compress ∘ K。
从 F ∘ compress = F开始，如果 list是 K的固定点，则 compress(list)是 L的固定点。相反，如果 clist是 L的固定点，则 K(clist)是 K和 compress的固定点。 K是 L resp的不动点集之间的双射。 K。 (并且 H(clist)是 H的不动点，而 compress ∘ sort和 H是 L和 H的不动点集之间的双射。)
最多 d位的压缩排序列表的空间足够小，足以强行考虑所考虑的 f函数，即幂函数。
因此，策略是:

查找要考虑的最大位数d(由于模数而为20，对于较小的q较小)。

生成最多d位的压缩单调序列。

检查序列是否是L的不动点，如果是，F(sequence)是G的不动点，即问题的解决方案。

代码
幸运的是，您还没有指定语言，所以我选择了最简单的代码，即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]

它不是最佳化的，但从现在开始，它会在不到50分钟的时间内从我的盒子中找到 2 <= q <= 10000的所有解决方案，从

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]

并以

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]

从约10到63的指数花费的时间最长(个别地，不是累积的)，由于减少了搜索空间，因此从64指数开始的速度显着提高。