6

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

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，但这是我朋友很久以前告诉我的一个编程任务。现在我想起了这一点，并想知道如何才能做到这一点。如何解决这个问题？

4

2 回答 2

4

有两个关键点，

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

让我们仔细看看这个案例q = 2。如果一个d-digit 数字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 <= 3or 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数字的 -th 次方的总和，边界将大约为

q - (q/20) + 1

所以对于更大q的，从中获得的可能解决方案的范围是巨大的。

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

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

由于单调数字d序列的数量相对较少，因此使用它的蛮力变得可行。特别是如果我们忽略对总和没有贡献的数字（0 表示所有指数，8 表示q >= 22，4 也表示q >= 32，所有偶数表示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，我们只需要考虑由奇数位组成的序列，并且长度单调序列的总数最多为 20使用 5 个符号是binom(20+5,20) - 1 = 53129. 现在，看起来不错。

概括

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

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

wheredigits映射n到它的数字列表。

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

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

然后我们正在寻找 的不动点G = F ∘ digits。Nown是 的一个不动点G当且仅当digits(n)是 的一个不动点H = digits ∘ F。因此，我们可以等效地寻找 的不动点H。

但是F是排列不变的，所以我们可以将自己限制在排序列表中并考虑K = sort ∘ digits ∘ F.

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

如果一些f(d)为 0（模 2 64），则可以进一步改进。让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。K. （并且H(clist)是 的不动点H，而compress ∘ sort分别H是 的不动点集之间的双射L。H）

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

所以策略是：

  1. 找到d要考虑的最大位数（由于模数而以 20 为界，对于 small 则更小q）。
  2. 生成最多d位数的压缩单调序列。
  3. 检查序列是否是 的不动点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]

它没有经过优化，但它2 <= q <= 10000在我的盒子上找到了不到 50 分钟的所有解决方案，从

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 开始有显着的加速。

于 2012-04-25T12:41:58.353 回答
0

这是一个蛮力解决方案，它将解决所有此类 n，包括 1 和任何其他大于您选择的任何范围内的第一个的 n（在这种情况下，我选择 base^q 作为我的范围限制）。您可以修改以忽略 1 的特殊情况并在第一个结果之后返回。它在 C# 中，但在带有 ** 幂运算符的语言中可能看起来更好。你也可以传入你的 q 和 base 作为参数。

int q = 5;
int radix = 10;
for (int input = 1; input < (int)Math.Pow(radix, q); input++)
{
    int sum = 0;
    for (int i = 1; i < (int)Math.Pow(radix, q); i *= radix)
    {
        int x = input / i % radix;       //get current digit
        sum += (int)Math.Pow(x, q);      //x**q;
    }
    if (sum == input)
    {
        Console.WriteLine("Hooray: {0}", input);
    }
}

因此，对于 q = 5，结果是：

Hooray: 1
Hooray: 4150
Hooray: 4151
Hooray: 54748
Hooray: 92727
Hooray: 93084
于 2012-04-24T13:20:57.330 回答