約数和を使った漸化式
Julia 言語は型を指定しないで書くと、色々な型に自動で対応してくれるというので、多倍長を試してみようと分割数を計算してみることにします。
まず Fortran で以前のプログラムを少し改良しておきます。Fortran のプログラムを移植します。
整数の分割数 と約数和関数 の関わる漸化式。
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 以下の分割数を全て求めているので、このように列挙目的の場合は毎回計算するのはとても無駄な気もしますw
Modern Fortran: Building efficient parallel applications
- 作者:Curcic, Milan
- 発売日: 2020/11/24
- メディア: ペーパーバック
How to learn Applied Mathematics through modern FORTRAN
- 作者:Hernandez Ramos, Juan Antonio,Escoto Lopez, Javier
- 発売日: 2020/02/23
- メディア: ペーパーバック
julia
よく文法を知らないので検索しながら書いてみましたが、まだ不明な点が多く、型を与えないと Float64 になってしまい、うまく Int64 や BigInt に出来ませんでした。また Julia の配列はそのままでは Fortran のように配列の始まりを指定できないようなので、漸化式をずらして書く羽目になりました。また総和関数 Σ の扱いの所で Fortran のように部分配列積にする方法が分かりませんでした。Julia では index をあらわに書くのでこっちの方が分かり良いかもしれません。(いつも数学書で総和の index の範囲が省略されているとあらわに求めるのに四苦八苦するのでw)
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
python
整数列辞典によれば python の sympy には数論関数が用意されていて分割数もあるようです。
from sympy.ntheory import npartitions print([npartitions(i) for i in range(317)])
The Encyclopedia of Integer Sequences
- 作者:Sloane, N. J.A.,Plouffe, Simon
- 発売日: 1995/01/15
- メディア: ハードカバー