「ウォルステンホルム素数…」を読む

をちゃんと計算してみたらできたのでまとめました。

素数生成

まずは素数を取り出す計算が必要。
といっても、AtCoderでも先の更新で入ったライブラリで

import Data.Numbers.Primes

とやれば、primes :: Integral int => [int] で手に入るので、今更何かをする必要はない…のだけど、何やら不穏なコードが貼り付けられている。

元ネタは、「エラトステネスの小型高速遠心分離機」というゴルフな人の記事。

はじまりのコード：

p=1#2;a#b|gcd a b<2=b:(a*b)#(b+1)|0<1=a#(b+1)

このままでは読めないのでdegolfする。gcdが2未満とは、1 つまり互いに疎だという意味。

primes = 1 # 2
  where
     a # b
       | gcd a b == 1 = b : (a*b) # succ b
       | otherwise    =      a    # succ b

動きを追ってみると：

1 # 2 => 2 : 2 # 3
2 # 3 => 3 : 6 # 5
6 # 5 => 5 : 30 # 6
30 # 6 => 30 # 7 => 7 : 210 # 8

つまり、aはこれまでに発見した素数全ての積で、これと互いに素な値が見つかるまでbを1ずつ増やす。
見つかったらそれは次の素数で、篩であるaはさらにbも掛け合わせて次に進む。

なるほど、ある数 $x$ が素数かどうかの判定は、$\sqrt x$ 未満の素数全てで割り切れないことで判断できる。これは$\sqrt x$ 未満の素数の個数だけの mod の実行を必要とする。
その代わりに、それら全ての積との $\gcd$ が1なら、どの素数とも互いに素であることが一回の gcd の呼び出しだけで計算できるということ。
これは一見凄いが、変数 a は急速に増加し、a に対する gcd は内部で何度も mod を呼ぶので、まぁまぁ大変なことになると思われる。また a は Int ではすぐ溢れるので、Integer が必須になる。

記事には、ここから発展した版

p=2:3:1#p;m#(a:b:x)=[n|n<-[a^2..b^2-2],odd n,gcd n m<2]++(m*b)#(b:x)

の、odd n の代わりに、「5以上の素数は6の倍数±1」を使った変種

p=2:3:5:5#p;m#(a:b:x)=[n|n<-[a^2..b^2-2],(mod (n+1) 6)*(mod (n-1) 6)==0,gcd n m<2]++(m*b)#(b:x)

が書かれている。
けど、(mod (n+1) 6)*(mod (n-1) 6)==0 と二度もmodを呼ぶより elem (mod n 6) [1,5] のが安いんじゃないかなぁという気もする。
このコードは説明しづらいので飛ばして、元記事の最終版

p=2:3:5#p;n#x@(m:p:y)=[n|gcd m n<2]++(n+2)#last(x:[m*p:y|p^2-3<n])

に進む。degolfする。(#)の右項はリストだが、先頭はそれ以降と明らかに意味が違うので分けて、3引数関数 sv に置き換える。

primes = 2 : 3 : sv (head primes) (tail primes) 5
  where
    sv m (p:ps) n
      | gcd m n == 1 = n : next (n + 2)
      | otherwise    =     next (n + 2)
      where
        next | p^2 - 3 < n = sv (m * p) ps
             | otherwise   = sv  m   (p:ps)

n は5以降の奇数を順に試し、合格なら素数と認められる。
m は最初の版の a と同じで、n に対して調べるべき範囲の素数全ての総積を保持する。next のガードで、必要になったら次の素数を継ぎ足している。

最初の版よりは m の増加速度は小さいだろうけど、結局、巨大な素数の積を相手に gcd をするという方針は同じ。

コンビネーション

${}_nC_k = n! / \{k! \cdot (n-k)!\}$ を使うと

comb n k = div (product [succ k .. n]) (product [1 .. n - k])

ではあるが、これを個々の素数について毎回全ての乗算を計算し直すのは乱暴な気がする。

連続する素数 $p, q$ について、${}_{2q - 1}C_{q-1}$ を求めるのに、一つ前の値 ${}_{2p - 1}C_{p-1}$ を利用できないか。
${}_nC_k = \frac{n}{k} \cdot {}_{n-1}C_{k-1}$ の適用を繰り返すと
${}_nC_k = \frac{n}{k} \cdot \frac{n-1}{k-1} \cdot \ldots \cdot \frac{n-k+1}{1}$ となる。これを使って

${}_{2p-1}C_{p-1} = (p+1) (p+2) \dots (p + p-1) / 1 \cdot 2 \cdot \ldots \cdot (p-1)$
${}_{2q-1}C_{q-1} = (q+1) (q+2) \dots (q + q-1) / 1 \cdot 2 \cdot \ldots \cdot (q-1)$

$p < q$ であることに注意して二つを見比べる。
まずは、分母に $p \cdot \ldots \cdot (q-1)$ が増えているとわかる。(A)とする。
分子の前の方は $(p+1) \cdot \ldots \cdot q$ が余計なので、これは分母に掛ける。(B)とする。
また後ろの方、$(p + p) \cdot \ldots \cdot (q + q - 1)$ が増えている。
これは $\{ (p + p) \cdot \ldots \cdot (p + q - 1) \} \times \{(q + p) \cdot \ldots \cdot (q + q - 1)\}$ と前後に分けて(C),(D)とすると、(A)と(C)、(B)と(D)の個数が同じで、結局

$${}_{2q-1}C_{q-1} = {}_{2p-1}C_{p-1} \times \prod_{i=p}^{q-1} \frac{p+i}{i} \times \prod_{j=p+1}^{q} \frac{q-1+j}{j}$$

という漸化式が得られた。

毎回${}_{2p-1}C_{p-1}$ を計算する乗算の総回数 $\sum p \simeq O(N^2)$ に対して、$q - p$に比例する計算量で次の項が得られるということは、$O(N)$ への高速化が予想される。

なお、この分析は、分子の項に重なりがあるかせめて隔たりはないことを想定している。
$q + 1 \leq p + p - 1$ なら少なくとも一つの項は再利用される。
$q + 1 = p + p$ ならちょうど重なりがないが、間隙もないので上式は成立している。
これについては「ベルトランの仮説」から $q < 2p$ つまり $q + 1 \leq 2p$ は保証され、助かった。

プログラミング

Data.Ratio を使って以上の検討結果をコード化しよう。

素数 $p$ に対して ${}_{2p-1}C_{p-1}$ を求める関数

import Data.Ratio

comb2p1p1 :: Integer -> Rational
comb2p1p1 p = product [(p + i) % i | i <- [1 .. pred p]]

素数リストを使って、それらの ${}_{2p-1}C_{p-1}$ を順次生成する無限リスト

import Data.Numbers.Primes

comb2p1p1s :: [Integer]
comb2p1p1s = map numerator $ scanl step (comb2p1p1 2) $ zip primes $ tail primes
  where
    step c (p, q) = product $ c : [(p + i) % i | i <- [p .. pred q]] ++
                                  [(pred q + j) % j | j <- [succ p .. q]]

$10^9$以下の素数 $p$ について、${}_{2p-1}C_{p-1} \bmod p^4 = 1$ なもののリスト

whps :: [Integer]
whps = [p | (p, c) <- zip (takeWhile (10^9 >=) primes) comb2p1p1s, mod c (p^4) == 1]

時間をはかるメインルーチン

import System.CPUTime

main = do
  t1 <- getCPUTime
  print ("prime1", whps !! 0)
  t2 <- getCPUTime
  print ("Time", t2 - t1)
  print ("prime2", whps !! 1)
  t3 <- getCPUTime
  print ("Time", t3 - t2)
--  print ("10^9", whps)
--  t4 <- getCPUTime
--  print ("Time", t4 - t3)

実行

CPUTimeの単位はピコ$(10^{-12})$秒。

ghci> main
("prime1",16843)
("Time",62500000000)

インタプリタでも一つめが0.0625秒で得られた。
が、ふたつめは無理そうだったのでコンパイルして実行。

> ghc -O2 3.hs
> ./3
("prime1",16843)
("Time",31250000000)
("prime2",2124679)
("Time",652109375000000)
>

一つめは0.03125秒。
二つめは652秒≒11分。

どうせ既知の結果だとはいえ、自分の計算機／プログラムでちゃんと答えが算出できると嬉しい。