計算法

以前メビウス関数の記事を書きましたが、Fortran での計算法を別のやり方でやってみます。

fortran66.hatenablog.com

ここでは素数の因数ごとに符号反転、素数の二乗以上のべきが混じるものは 0 にするという、愚直な方法でまず一覧表を作ってそれを利用することにします。

プログラム

算法：エラトステネスの篩で n 以下の素数の表を作り、これに基づき素数ごとにその倍数のメビウス関数の表の符号を反転し、さらに二乗三乗･･･の倍数のメビウス関数の表を 0 にクリアしてゆきます。

求まったメビウス関数の表を元に、メビウスの反転公式を使って素朴に計算したゼータ関数の逆数を 2,3 の場合に計算してみます。

    module moebius_m
        implicit none
    contains
        function iprime(n) result(ires) ! primes under n
            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)

        
   ! zeta(2) = pi^2/6
        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 *
        
    ! zeta(3) = Apery's constant
        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