別の興味深いプログラミング/数学の問題があります。
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 のどちらに当てはまるかはわかりませんでしたが、これはかなり前に友人から言われたプログラミング作業でした。今、私はそれを思い出し、そのようなことができる方法を知りたいと思っています. これにアプローチする方法は？
2つの重要なポイントがあります。
ケースを詳しく見てみましょうq = 2。d1桁の数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、最大3桁の二乗和は最大243であるため、解は244未満である必要があります。その範囲の桁の二乗和の最大値は199に達します。 1²+9²+9²=163、続けて、解は100未満でなければならないことが簡単にわかります。（の唯一の解q = 2は1です。）の場合q = 3、4桁の立方体の最大合計は4 * 729 = 2916であり、続き、すべてのソリューションがq = 3は1000未満です。しかし、この種の境界の改善は、モジュラス要件のために小さな指数に対してのみ役立ちます。桁の累乗の合計がモジュラスを超える可能性がある場合、それは崩壊します。したがって、可能な最大桁数を見つけるのをやめます。
ここで、モジュラスがないq場合、桁の5乗の合計の場合、境界はほぼ次のようになります。
q - (q/20) + 1
したがって、より大きなq場合、そこから得られる可能な解決策の範囲は膨大です。
しかし、ここで2つのポイントが救い出されます。1つは解空間を最大20桁に制限するモジュラス、もう1つ2 <= n < 2^64は（モジュラー）デジタルパワー和の順列不変性です。
順列不変性は、単調な数字のシーケンスを作成し、 5乗dの合計を計算して、このようにして得られた数字が正しい数字であるかどうかを確認するだけでよいことを意味します。q
単調な桁のdシーケンスの数は比較的少ないので、それを使用するブルートフォースが実行可能になります。特に、合計に寄与しない数字を無視する場合（すべての指数の場合は0、の場合は8 q >= 22、またの場合は4、のq >= 32場合はすべて偶数の数字q >= 64）。
シンボルdを使用した長さの単調なシーケンスの数は次のとおりです。s
binom(s+d-1, d)
sは最大で9であり、からd <= 20の合計で、各指数について考慮する必要のあるシーケンスは最大で10015004です。そんなに多くはありません。d = 1d = 20
それでも、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)
ここで、その桁のリストにdigitsマップされます。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。
Hとの不動点Kは1対1で対応しています。listがの不動点である場合、はの不動点でありH、がの不動点であるsort(list)場合K、はの順列であり、つまり、は、の不動点であり、respです。との不動点のセット間の全単射です。sortedListKH(sortedList)sortedListH(H(sortedList)) = H(sortedList)H(sortedList)KsortHHK
f(d)一部が0（モジュロ2 64 ）の場合、さらに改善することができます。数字のリストからmodを使用してcompress数字を削除する関数とし、関数を検討します。f(d) 2^64 == 0L = compress ∘ K
F ∘ compress = F、の場合listはの不動点なのでK、compress(list)はの不動点ですL。逆に、clistがの不動点である場合L、K(clist)は、の不動点でありK、compressrespです。Kそれぞれの不動点のセット間の全単射Lです。K。（そして、はのH(clist)不動点であり、それぞれは、それぞれの不動点のセット間の全単射です。）Hcompress ∘ sortHLH
最大で桁数の圧縮されたソート済みリストのスペースは、検討中の関数、つまり電力関数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]
最適化されていませんが、現状では、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以降で大幅な高速化が実現します。
これは、1 と、選択した範囲内の最初の n より大きいその他の n を含む、そのようなすべての 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