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

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,但这是我朋友很久以前告诉我的编程任务.现在我记得那个,并想知道如何做到这一点.怎么解决这个问题？

Answer 1

有两个关键点，

\n

\n
• 可能的解决方案的范围是有限的，
\n
• 任意一组数字在排列 con 之前都相同的数字至多包含一个解。
\n

\n

让我们仔细看看这个案例q = 2。如果一个d位数n等于其各位数字的平方和，则

\n

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

\n

并且条件10^(d-1) <= d*81很容易转换为d <= 3or n < 1000。需要检查的数字并不多，对这些数字进行暴力破解很快。对于q = 3，条件10^(d-1) <= d*729产生d <= 4，但仍然没有多少数字需要检查。通过进一步分析，我们可以找到更小的界限，对于q = 2，最多三位数字的平方和最多为 243，因此解必须小于 244。该范围内的数字平方和最大达到 199： 1\xc2\xb2 + 9\xc2\xb2 + 9\xc2\xb2 = 163，继续下去，我们很容易发现解一定小于100。（ 的唯一解是q = 21。）对于q = 3， 的最大和四个数字的立方是 4*729 = 2916，继续，我们可以看到 的所有解都q = 3小于 1000。但是由于模数要求，这种边界的改进仅对小指数有用。当数字的幂之和超过模数时，它就会崩溃。因此我停止寻找最大可能的位数。

\n

现在，如果没有模数，对于q数字的 - 次幂之和，界限将大约为

\n

q - (q/20) + 1\n

\n

因此q，对于更大的问题，从中获得的可能解决方案的范围是巨大的。

\n

但这里有两点可以解决问题，首先是模数，​​它将解空间限制为2 <= n < 2^64最多 20 位数字，其次是（模）数字幂和的​​排列不变性。

\n

排列不变性意味着我们只需要构造单调的d数字序列，计算q第-次幂的和，并检查由此得到的数字是否具有正确的数字。

\n

由于单调数字序列的数量d相对较小，因此使用暴力破解变得可行。特别是如果我们忽略对总和没有贡献的数字（0 表示所有指数，8 表示q >= 22，4 表示q >= 32，所有偶数表示q >= 64）。

\n

d使用符号的长度单调序列的数量s是

\n

binom(s+d-1, d)\n

\n

s对我们来说最多 9, d <= 20，从d = 1到求和d = 20，每个指数最多有 10015004 个序列需要考虑。这还不算太多。

\n

尽管如此，对所有q考虑的数字这样做仍然需要很长时间，但是如果我们考虑到对于q >= 64，对于所有偶数数字x^q % 2^64 == 0，我们只需要考虑由奇数数字组成的序列，并且长度最多为 20 的单调序列总数使用 5 个符号是binom(20+5,20) - 1 = 53129. 现在，看起来不错。

\n

概括

\n

我们考虑一个f将数字映射到自然数的函数，并寻找方程的解

\n

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

\n

其中digits映射n到其数字列表。

\n

从 开始f，我们构建了一个F从数字列表到自然数的函数，

\n

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

\n

然后我们寻找 的不动点G = F \xe2\x88\x98 digits。Nown是 的不动点G当且仅当digits(n)是 的不动点H = digits \xe2\x88\x98 F。因此我们可以等效地寻找 的不动点H。

\n

但F是排列不变的，因此我们可以将自己限制在排序列表并考虑K = sort \xe2\x88\x98 digits \xe2\x88\x98 F。

\n

的不动点H和 ofK是一一对应的。如果list是 的不动点H，则sort(list)是 的不动点K，如果sortedList是 的不动点K，则H(sortedList)是 的排列sortedList，因此H(H(sortedList)) = H(sortedList)，换句话说，H(sortedList)是 的不动点K，并且sort分别。是和H的一组不动点之间的双射。HK

\n

f(d)如果某些值为 0（模 2 ⁶⁴），则可以进一步改进。设是一个从数字列表中compress删除带有f(d) mod 的数字的函数，并考虑该函数。 2^64 == 0L = compress \xe2\x88\x98 K

\n

因为F \xe2\x88\x98 compress = F，如果list是 的不动点K，则compress(list)是 的不动点L。相反，如果clist是 的不动点L，则K(clist)是 的不动点K，并且compress分别是。K分别是固定点集之间的双射L。K。（H(clist)是 的不动点H， 是compress \xe2\x88\x98 sort的H不动点集之间的双射L。H）

\n

最多d位数的压缩排序列表的空间足够小，足以对f所考虑的函数（即幂函数）进行暴力破解。

\n

所以策略是：

\n

\n
1. 找到要考虑的最大位数d（由于模数的限制，以 20 为界，小 则更小q）。
\n
2. 生成最多d位数的压缩单调序列。
\n
3. 检查序列是否为 的不动点L，如果是，F(sequence)则为 的不动点G，即问题的解。
\n

\n

代码

\n

幸运的是，你没有指定一种语言，所以我选择了最简单的代码，即 Haskell：

\n

{-# LANGUAGE CPP #-}\nmodule Main (main) where\n\nimport Data.List\nimport Data.Array.Unboxed\nimport Data.Word\nimport Text.Printf\n\n#include "MachDeps.h"\n\n#if WORD_SIZE_IN_BITS == 64\n\ntype UINT64 = Word\n\n#else\n\ntype UINT64 = Word64\n\n#endif\n\nmaxDigits :: UINT64 -> Int\nmaxDigits mx = min 20 $ go d0 (10^(d0-1)) start\n  where\n    d0 = floor (log (fromIntegral mx) / log 10) + 1\n    mxi :: Integer\n    mxi = fromIntegral mx\n    start = mxi * fromIntegral d0\n    go d p10 mmx\n        | p10 > mmx = d-1\n        | otherwise = go (d+1) (p10*10) (mmx+mxi)\n\nsortedDigits :: UINT64 -> [UINT64]\nsortedDigits = sort . digs\n  where\n    digs 0 = []\n    digs n = case n `quotRem` 10 of\n               (q,r) -> r : digs q\n\ngenerateSequences :: Int -> [a] -> [[a]]\ngenerateSequences 0 _ \n    = [[]]\ngenerateSequences d [x] \n    = [replicate d x]\ngenerateSequences d (x:xs)\n    = [replicate k x ++ tl | k <- [d,d-1 .. 0], tl <- generateSequences (d-k) xs]\ngenerateSequences _ _ = []\n\nfixedPoints :: (UINT64 -> UINT64) -> [UINT64]\nfixedPoints digFun = sort . map listNum . filter okSeq $ \n                        [ds | d <- [1 .. mxdigs], ds <- generateSequences d contDigs]\n  where\n    funArr :: UArray UINT64 UINT64\n    funArr = array (0,9) [(i,digFun i) | i <- [0 .. 9]]\n    mxval = maximum (elems funArr)\n    contDigs = filter ((/= 0) . (funArr !)) [0 .. 9]\n    mxdigs = maxDigits mxval\n    listNum = sum . map (funArr !)\n    numFun = listNum . sortedDigits\n    listFun = inter . sortedDigits . listNum\n    inter = go contDigs\n      where\n        go cds@(c:cs) dds@(d:ds)\n            | c < d     = go cs dds\n            | c == d    = c : go cds ds\n            | otherwise = go cds ds\n        go _ _ = []\n    okSeq ds = ds == listFun ds\n\n\nsolve :: Int -> IO ()\nsolve q = do\n    printf "%d:\\n    " q\n    print (fixedPoints (^q))\n\nmain :: IO ()\nmain = mapM_ solve [2 .. 10000]\n

\n

它没有经过优化，但就目前情况而言，它2 <= q <= 10000在我的盒子上在不到 50 分钟的时间内找到了所有解决方案，从

\n

2:\n    [1]\n3:\n    [1,153,370,371,407]\n4:\n    [1,1634,8208,9474]\n5:\n    [1,4150,4151,54748,92727,93084,194979]\n6:\n    [1,548834]\n7:\n    [1,1741725,4210818,9800817,9926315,14459929]\n8:\n    [1,24678050,24678051,88593477]\n9:\n    [1,146511208,472335975,534494836,912985153]\n10:\n    [1,4679307774]\n11:\n    [1,32164049650,32164049651,40028394225,42678290603,44708635679,49388550606,82693916578,94204591914]\n

\n

并以

\n

9990:\n    [1,12937422361297403387,15382453639294074274]\n9991:\n    [1,16950879977792502812]\n9992:\n    [1,2034101383512968938]\n9993:\n    [1]\n9994:\n    [1,9204092726570951194,10131851145684339988]\n9995:\n    [1]\n9996:\n    [1,10606560191089577674,17895866689572679819]\n9997:\n    [1,8809232686506786849]\n9998:\n    [1]\n9999:\n    [1]\n10000:\n    [1,11792005616768216715]\n

\n

从大约 10 到 63 的指数花费的时间最长（单独计算，而不是累积），由于搜索空间减少，从指数 64 开始有显着的加速。

\n