計算法
以前メビウス関数の記事を書きましたが、Fortran での計算法を別のやり方でやってみます。
fortran66.hatenablog.com
ここでは素数の因数ごとに符号反転、素数の二乗以上のべきが混じるものは 0 にするという、愚直な方法でまず一覧表を作ってそれを利用することにします。
プログラム
算法:エラトステネスの篩で n 以下の素数の表を作り、これに基づき素数ごとにその倍数のメビウス関数の表の符号を反転し、さらに二乗三乗・・・の倍数のメビウス関数の表を 0 にクリアしてゆきます。
求まったメビウス関数の表を元に、メビウスの反転公式を使って素朴に計算したゼータ関数の逆数を 2,3 の場合に計算してみます。
module moebius_m
implicit none
contains
function iprime(n) result(ires)
integer, intent(in) :: n
integer, allocatable :: ires(:)
integer :: itab(n), i
itab(1) = 0
forall(i = 2:n) itab(i) = i
do i = 2, int(n**0.5)
if (itab(i) /= 0) itab(i**2::i) = 0
end do
ires = pack(itab, itab /= 0)
end function iprime
subroutine moebius_table(moebius)
integer, intent(out) :: moebius(:)
integer, allocatable :: iprime_table(:)
integer :: i, j, k, n
n = size(moebius)
iprime_table = iprime(n)
moebius = 1
do i = 1, size(iprime_table)
k = iprime_table(i)
moebius(k::k) = -moebius(k::k)
j = k**2
do while(j <= n)
moebius(j::j) = 0
j = j * k
end do
end do
end subroutine moebius_table
end module moebius_m
program MoebiusP
use moebius_m
implicit none
real, parameter :: pi = 4 * atan(1.0)
integer, allocatable :: m(:)
integer :: i, n
real :: x, s
integer :: mu(10**5)
call moebius_table(mu)
x = 0.0
s = 0.0
do i = 1, 10000
x = x + mu(i) / real(i)**2
s = s + 1.0 / real(i)**2
end do
print *, x, 6.0 / pi**2, s, pi**2 / 6
print *, 'zeta(2)^-1 * zeta(2) =', x * s
print *
x = 0.0
s = 0.0
do i = 1, 1000
x = x + mu(i) / real(i)**3
s = s + 1.0 / real(i)**3
end do
print *, x, 1.0 / 1.20205690, s, 1.20205690
print *, 'zeta(3)^-1 * zeta(3) =', x * s
end program MoebiusP
計算結果
以前の結果を再現しています。一気にテーブルを作るので前より計算結果が早く出てきます。
0.6079294 0.6079271 1.644725 1.644934
zeta(2)^-1 * zeta(2) = 0.9998769
0.8319076 0.8319074 1.202051 1.202057
zeta(3)^-1 * zeta(3) = 0.9999951
オンライン整数列辞典
oeis.org