Meertens数
Posted feedbacks - Haskell
Haskell の勉強にと移植してみました。GHCi で33秒とあんまり速くはありませんが、わりとすっきりしたコードになったので投稿してみます。ついでに中身について少し説明でも。
上の桁から順に決めていくのですが、その際に桁数を固定することでゲーデル数化する前と後の数(sum, prod)の両方を部分的に計算でき、枝刈りがやりやすくなります。
実際にやっているのは次の二つですが、他にも末尾の 0 の数で刈れたりするかもしれません。
- 積>和 になったら打ち切り(積のほうが速く増えるので)
- bound にまだ決まっていない部分の積の上限を保持し、それを超えたら打ち切り
後者のは実質「積が 10^i を超えたら」なのですが、上のようにした方が速いようです。(扱う数が fixnum に収まりやすい？)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | primes = flip take [2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71]
basis n = map (\x -> 10^x) [n-1,n-2..0]
f :: [Integer] -> [Integer] -> Integer -> Integer -> Integer -> Integer -> [Integer]
f [] [] _ prod sum _
| prod == sum = [sum]
| otherwise = []
f (p:ps) (b:bs) bound prod sum start
| start == 0 && prod > sum = []
| otherwise =
[x | i <- [start..9], p^i <= bound,
x <- f ps bs (div bound $ p^i) (prod * p^i) (sum + b * i) 0]
meertens n = [x | i <- [1..n], x <- f (primes i) (basis i) (10^i) 1 0 1]
|
ghcでコンパイルすると 桁数 : 時間 10 : 0.6s 11 : 1.9s 12 : 6.1s 13 : 18.0s 14 : 52.4s 15 : 150.1s 1桁増えるごとに約3倍 20桁なら10時間かな? Meertens数は81312000しか知られていないらしいです。 他にあるかどうかも未知らしい。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | module Main (main) where
import Data.Tree
import System.Environment
primes :: [Integer]
primes = sieve [2..] where sieve (p:xs) = p:sieve [x | x <- xs, x `mod` p /= 0]
type Label = (Integer, Integer)
powers :: Integer -> [Integer]
powers = flip iterate 1 . (*)
mktree :: Int -> Label -> [[Integer]] -> Tree Label
mktree k l pss = unfoldTree phi (l,pss)
where phi (x,[]) = (x,[])
phi (x,ps:pps) = (x,zip (labels x ps) (repeat pps))
labels (n,g) ps = zip (map (10*n +) [0..9]) (takeWhile (10^k >) (map (g*) ps))
gtree :: Int -> Tree Label
gtree k = snip $ mktree k (0,1) $ map powers $ take k primes
where snip (Node x ts) = Node x (tail ts)
meertens :: Int -> [Integer]
meertens k = [ n | (n,g) <- flatten $ gtree k, n == g ]
main :: IO ()
main = putStr . unlines . map show . meertens . read . head =<< getArgs
|
nobsun
#4710()
Rating4/4=1.00
お題#100「正整数のゲーデル数化？」で定義した goedel を適用すると自分自身になるような数，すなわち goedel (n) == n となるような正整数 n を見つける関数を定義してください．
このような数のことをMeertens数と言うそうです．
32bitsの符号なし整数(あるいは10進10桁整数)までの範囲で探すのにどのくらい計算時間がかかったかをCPUのスペックとともに教えてください．また，その実装で64bit符号なし整数(あるいは10進20桁整数)までの範囲で探すのにどのくらい計算時間がかかりそうか見積ってください(もちろん実際に計算して計算時間を示していただくのでもかまいません)．
1 reply [ reply ]