メビウス関数 μ(n) - fortran66のブログ

メビウス関数

数論入門とか組み合わせ論入門とかを読んでゆくと、メビウス関数というものがメビウスの反転公式とともに出てくるのですが、書いてあることがちゃお付録のおまじないブック並みに意味不明で、やる気を無くします。

http://kids.nifty.com/cs/kuchikomi/kids_bbs/list/aid_omajinai-love-01/1.htm

メビウス関数の定義

$\mu(n) \equiv \left\{ \begin{array}{ll} 0 & (n が平方因子をもつ[同じ素数で二回(以上)割りきれる]) \\ +1 & (n が相異なる偶数個の素因数に分解される) \\ -1 & (n が相異なる奇数個の素因数に分解される) \end{array} \right.$

メビウスの反転公式

このメビウス関数 $\mu(n)$ を使うと、
${\displaystyle g(n)=\sum_{d|n}f(d), }$
のとき、
${\displaystyle f(n)=\sum_{d|n}\mu({n\over d})g(d) }$
となる。

以下で、やる気を出して考えてみることにします。

整数の素因数分解と分割

整数をより小さな単位に分解するとき、

素因数分解（積による分解）
分割（和による分解）

の二つが考えられます。

分割の方は、オイラーの分割恒等式などで出てくるように、べき級数 ( $\sum_na_nx^n$ ) の母関数というものを利用してべき指数の和の問題として解くことができます。

これと同じように、母関数に指数的な級数 (ディリクレ級数 $\sum_n a_n/n^x$ ) をとることで、積の問題を解くことができ、これからメビウス関数と反転公式が出てきます。

ここで、ディリクレ級数の形が分数になっていて、ちょっとイライラするのですが、元々は $\sum_na_nx^n$ と $\sum_na_nn^x$ のように似た者同士でしたが、収束性を考えると、前者は $|x|<1$ で収束しますが、後者は $x<-1$ と負数でないと収束しないので、それならば初めから負数と思って、分母に置くことが考えられます。ディレクレが老婆心をだしてそうしてくれたと思うと、許してやろうという気にもなります。

さて、オイラー積を考えて、全ての素数に関する積をとると、
${\displaystyle \prod_{p_i}(1-p_i^{-s})^{-1}=\prod_{p_i}(1+p_i^{-s}+p_i^{-2s}+p_i^{-3s}+\cdots)=\sum_n n^{-s}=\zeta(s) }$
と、素因数分解のすべての組み合わせが出るので、全正整数に関する和となってゼータ関数が出てきます。

ところで、最初の積の逆数をとると、素数のダブりが許されない積になるのですべての整数は出てきません。また掛け算回数の偶奇性によって符号も変わることになります。これが丁度メビウス関数の定義に相当しています。つまり、
${\displaystyle \prod_{p_i}(1-p_i^{-s})=\sum_n\mu(n)n^{-s} }$
となります。

そうして、これらの式の左辺が丁度逆数の関係になっているので、右辺の級数同士も逆数になっているはずです。すると、ゼータ関数の逆関数は、メビウス関数が係数となったディリクレ級数であることが分かります。つまり、
${\displaystyle {\zeta(x)}^{-1}=\sum_n{\mu(n)\over n^x,} }$
となります。

また、一般のディリクレ級数
${\displaystyle A(x)=\sum_n{a_n\over n^x,} }$
${\displaystyle B(x)=\sum_n{b_n\over n^x} }$
に対して、
${\displaystyle A(x)\zeta(x)=B(x) }$
すなわち、
${\displaystyle A(x)=B(x)\zeta^{-1}(x) }$
の関係が成り立つとすると、
${\displaystyle A(s)\zeta(s)=(\sum_n{a_n\over n^x})(\sum_n{1\over n^x})={a_1\over1^s}+{a_2+a_1\over2^s}+{a_3+a_1\over3^s}+{a_4+a_2+a_1\over4^s}+\cdots }$
より、対応する指数の項を比較して、
${\displaystyle b_n=\sum_{d|n}a_n }$
となります。一方、逆関数の式からは、
${\displaystyle B(s)\zeta^{-1}(s)=(\sum_n{b_n\over n^x})(\sum_n{\mu(n)\over n^x})={b_1\mu_1\over1^s}+{b_2\mu(1)+b_1\mu(2)\over2^s}+{b_3\mu(1)+b_1\mu(3)\over3^s}+{b_4\mu(1)+b_2\mu(2)+b_1\mu(4)\over4^s}+\cdots }$
となるので、対応する指数の項の比較から、
${\displaystyle a_n=\sum_{d|n}b_d\mu({n\over d}) }$
または、
${\displaystyle a_n=\sum_{d|n}b_{n\over d}\mu(d) }$
が求まります。

以上により、母関数の方法によりメビウスの反転公式が示されました。

結局、分割の問題も、基本対称関数と同次対称関数の関係も、母関数とその逆関数の積およびそれらの級数展開の積の関係から求められており、同じ論理構造を持っていることが分かります。

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

作者:Metcalf, Michael,Reid, John,Cohen, Malcolm
Oxford University Press, USA

Amazon

Fortran でメビウス関数を使ってゼータ関数の逆数を求める

メビウス関数の求め方ですが、どの求め方が良いのかよく分からないので、
${\displaystyle \mu(n)=\sum^n_{k=1, (k,n)=1}\exp({2\pi i k\over n}) }$
で求めます。これがメビウス関数になることは、１の n 乗根の図を描いて考えると示せます。複素数の和ですが虚部は消えるので、 $\cos$ を使って実数だけで、
${\displaystyle \mu(n)=\sum^n_{k=1, (k,n)=1}\cos({2\pi k\over n}) }$
と書きかえられます。n が大きくなると結構重い計算になります。本来は、乗法的関数なのでもっと簡単に計算できるはずですが･･･めんどくさいので力づくでｗ

R4.7.24 追記:　素朴計算法
fortran66.hatenablog.com

ソース・プログラム

    module m_moebius
        implicit none
        real, parameter :: pi = 4 * atan(1.0)
    contains
        pure elemental integer function igcd(ia, ib)
            integer, value :: ia, ib
            igcd = ia
            do 
                if (ib == 0) exit
                igcd = ib
                ib = mod(ia, ib)
                ia = igcd
            end do
        end function igcd

        
        pure elemental integer function moebius(n) 
            integer, intent(in) :: n 
            integer :: k
            real :: x
            x = 0.0
            do k = 1, n
                if (igcd(k, n) == 1) x = x + cos(2 * pi * k / n)  
            end do
            moebius = nint(x)
        end function moebius
    end module m_moebius
    
    
    program zeta
        use m_moebius
        implicit none
        integer, allocatable :: m(:), moe(:)
        integer :: i, n
        real :: x, s
   ! zeta(2) = pi^2/6
        x = 0.0
        s = 0.0
        do i = 1, 10000
            x = x + moebius(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 + moebius(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 zeta

実行結果

定義通りに和を取っています。単精度なので、項の大きさが $10^{-8}$ 程度で打ち切っています。 $\zeta(2)$ は収束が遅いので、沢山足した割に誤差が大き目です。

  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

ゼータ関数とその逆関数の積は１になるはずです。