【メモ帳】分割数の漸化式

約数和を使った漸化式

Julia 言語は型を指定しないで書くと、色々な型に自動で対応してくれるというので、多倍長を試してみようと分割数を計算してみることにします。

まず Fortran で以前のプログラムを少し改良しておきます。Fortran のプログラムを移植します。

fortran66.hatenablog.com

fortran66.hatenablog.com

整数の分割数 と約数和関数 の関わる漸化式。

fortran

約数和 の計算を二重ループでの和の計算法に変えました。n 以下の全ての約数和が必要になるので都合良いです。

1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
-----------------------------------------------
1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
   2     2     2     2     2     2     2     2
      3        3        3        3        3  
         4           4           4           4
            5              5              5
               6                 6
                  7                    7
                     8                       8
                              ........

σ:　sum ↓
1  3  4  7  6 12  8 15 ....    

    module sub_m
        use, intrinsic :: iso_fortran_env
        implicit none
        integer, parameter :: ki = int64
    contains
        integer(ki) function p3(n)
            integer(ki), intent(in) :: n
            integer(ki) :: i, j, k, isigma(n), p(0:n)
            isigma = 0
            do i = 1, n
                isigma(i::i) = isigma(i::i) + i
            end do 
            !
            p(0) = 1
            do i = 1, n
                p(i) = sum( isigma(:i) * p(i - 1:0:-1) ) / i
            end do  
            p3 = p(n)
        end function p3
    end module sub_m
    
    program test
        use sub_m
        implicit none
        integer(ki) :: i
        do i = 1, 333
            print *, i, p3(i)
        end do
    end program test

Int64 では 317 でオーバーフローしてしまいます。

                   311     20017426762576945
                   312     21458096037352891
                   313     23000006655487337
                   314     24650106150830490
                   315     26415807633566326
                   316     28305020340996003
                   317    -27865439693783381
                   318    -25703328332159691
                   319    -23755504567687875
                   320    -21827600113664625
                   321    -20272407359078638
                   322    -18354279056910290
                   323    -17330349751739242
                   324    -15385735921416342
                   325    -14497628583594726
                   326    -13031122453577407
                   327    -12249745943421815
                   328    -10162837954485555
                   329    -10706307434721756
                   330     -8505358658431080
                   331     -7875939912433202
                   332     -6951402927955873
                   333     -7066008389406727

漸化式を求める都合、暗黙裡に n 以下の分割数を全て求めているので、このように列挙目的の場合は毎回計算するのはとても無駄な気もしますｗ

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Modern Fortran: Building efficient parallel applications

How to learn Applied Mathematics through modern FORTRAN

julia

よく文法を知らないので検索しながら書いてみましたが、まだ不明な点が多く、型を与えないと Float64 になってしまい、うまく Int64 や BigInt に出来ませんでした。また Julia の配列はそのままでは Fortran のように配列の始まりを指定できないようなので、漸化式をずらして書く羽目になりました。また総和関数 Σ の扱いの所で Fortran のように部分配列積にする方法が分かりませんでした。Julia では index をあらわに書くのでこっちの方が分かり良いかもしれません。（いつも数学書で総和の index の範囲が省略されているとあらわに求めるのに四苦八苦するのでｗ）

zeros(n) の所で型を指定しています。Int64 の時、p(229) で InexactError: Int64(4.413293488425501e13) を出して止まってしまいます。Overflow などへの安全策を取っているのでしょうか？よく分かりません。Fortran の p(316) は正しい結果になっています。

function p3(n)
    σ = zeros(BigInt, n)
    for i in 1:n
        for j in i:i:n
            σ[j] += i
        end
    end
    
    p = zeros(BigInt, n+1)
    p[1] = 1
    for i in 1:n
        p[i + 1] = sum( σ[j] * p[i + 1 - j] for j in 1:i ) / i
    end
    p[n + 1]
end

[p3(i) for i in 1:316]

316-element Array{BigInt,1}:
                 1
                 2
                 3
                 5
                 7
                11
                15
                22
                30
                42
                56
                77
               101
                 ⋮
 13162217895057704
 14118662665280005
 15142952738857194
 16239786535829663
 17414180133147295
 18671488299600364
 20017426762576945
 21458096037352891
 23000006655487337
 24650106150830490
 26415807633566326
 28305020340996003

BigInt を使うと p3(4895) で InexactError: が出ます。

p3(4894)

25107675557902586832045841006309659328628038105522320782464133131959268975

1から始める Juliaプログラミング

python

整数列辞典によれば python の sympy には数論関数が用意されていて分割数もあるようです。

oeis.org

from sympy.ntheory import npartitions

print([npartitions(i) for i in range(317)])  

The Encyclopedia of Integer Sequences