分割数の直接計算

整数の分割数を Hardy-Ramanujan-Rademacher の公式で計算して見ます。

Fortran, Wolfram-Language 等でやったものの Julia 版です。

fortran66.hatenablog.com

Hardy-Ramanujan-Rademacher の公式

$\displaystyle{ p(n)=\sum_{k=1}^\infty{\sqrt{k}\over\pi\sqrt{2}}\sum_{0\lt h\leq k\\ (h,k)=1}e^{{-2\pi i n h\over k}+\pi i\sum_{j=1}^{k-1}{j\over k}({h j\over k}-\lfloor{h j\over k}\rfloor-{1\over2})} {d\over dz}{\sinh\left({\pi\over k}\sqrt{{2\over3}(z-{1\over24})}\right)\over \sqrt{z-{1\over24}} } }$

julia lang プログラム

実数は倍精度 Float64 がデフォルトで、p(250) から結果が狂います。４倍精度は別途モジュールを読み込む必要があるようで、よく分からないのでやめときました。Julia の FORMAT 文で、複素数が表現できないみたいで草。

型を指定しないと気楽に書けて便利ですが、Fortran になれていると、行儀が悪くて気持ち悪いです。型変換がちゃんとされているか不安になります。

Julia 言語版は、今まで書いてきた言語の内で一番軽く書けて、ぱっと見で何をやっているかもわかると思います。

簡潔性では Wolfram Language ですが癖が強すぎます。でも微分を別途人力でやらずに、プログラムに記述できるのは便利だと思います。

Haskell は分割数漸化式版の方は簡潔に気持ちよく書けましたが、こういう数式には向いてないと思いました。where 節で内部関数を記述できるのは便利な気がしました。無限級数を無限とカットオフで書けるのも気分は良かったです。

でも幼女なので Fortran の方がも～っと好きです。

1996年頃のCM 松本引越センターぞうさんのほうが、もっと好きです

function fs(k, z)
    xs = sqrt(z-1.0/24.0)
    fs = pi / k * sqrt(2.0/3.0)
    (fs * cosh(fs * xs) - sinh(fs * xs) / xs) / (2 * xs * xs)
end

function fb(k, h)
    s = 0.0
    for j = 1:k-1
        s += j * (h * j / k - floor(h * j / k) - 1/2)
    end 
    s / k
end

function fa(n, k)
    s = 0.0+0.0im
    for h = 1:k 
        if (gcd(h,k)==1) 
           z = -2*pi*n*h*im /k + pi*im*fb(k, h)
           s += exp(z)
        end    
    end 
    return s
end

function p(n)
    s = 0.0+0.0im
    for k=1:10
        s += sqrt(k) * fa(n, k) * fs(k, n)
    end 
    s = s / (pi * sqrt(2)) 
end

using Printf
for i = 1:10
    pn = p(i)
    @printf( "%10d %20d %20.7f %20.7f \n", i, round(Int64, real(pn)), real(pn), imag(pn)) 
end

         1                    1            1.0063342            0.0000000 
         2                    2            2.0044611            0.0000000 
         3                    3            3.0023039            0.0000000 
         4                    5            4.9960379            0.0000000 
         5                    7            6.9972722            0.0000000 
         6                   11           10.9886684            0.0000000 
         7                   15           14.9997219            0.0000000 
         8                   22           21.9944150            0.0000000 
         9                   30           30.0038956            0.0000000 
        10                   42           42.0078618            0.0000000

f:id:fortran66:20190728153657p:plain — Julia Partition

fortran66のブログ

fortran について書きます。

【メモ帳】Julia で分割数

分割数の直接計算

Hardy-Ramanujan-Rademacher の公式

julia lang プログラム