C#:メビウス関数


このプログラムの原型を書いたのは随分と昔のことです。当時素数を勉強したいと思って、このメビウス関数のコードを書いてみたのですが、数学の素養が欠けている僕には、このメビウス関数が一体何なのかが、イマイチわからない...でもせっかく書いたのだからとWebサイトに公開したのが約7年前。すこし手直ししてここに再掲載することとします。


##メビウス関数

メビウス関数の定義は、以下のようなものです。
<pre>1) μ(1) = 1
2) μ(n) = 0 （n が平方因子を持つ（平方数で割り切れる）とき）
3) μ(n) = (-1)\*\*k （n が相異なる k 個の素因数に分解されるとき）
　　　　　　n が相異なる偶数個の素数の積ならば μ(n) = 1
　　　　　　n が相異なる奇数個の素数の積ならば μ(n) = -1
　　　　　　**は累乗を表す。</pre>

詳しくは、Wikipedia([メビウス関数](https://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%93%E3%82%A6%E3%82%B9%E9%96%A2%E6%95%B0))を参照してください。



## C#のコード

上の定義を元に、1から100までのメビウス関数の値を求めるコードを書いてみました。


```csharp:Mebius.cs
using System;
using System.Linq;

namespace Gushwell.Etude {
    class Program {
        static void Main(string[] args) {
            // 1 - 100 までのメビウス関数を求める
            int upper = 100;
            Mebius mebius = new Mebius(upper);
            for (int i = 1; i <= upper; i++) {
                Console.WriteLine("μ({0}) = {1}", i, mebius[i]);
            }
            Console.ReadLine();
        }
    }
    class Mebius {
        private int[] mebius;
        public Mebius(int maxnum) {
            mebius = Enumerable.Repeat(1, maxnum + 1).ToArray();
            foreach (int p in PrimeNumber.Enumerate().TakeWhile(n => n <= maxnum)) {
                for (int i = p; i <= maxnum; i += p)
                    mebius[i] *= -1;
                int p2 = p * p;
                for (int pp = p2; pp <= maxnum; pp += p2)
                    mebius[pp] = 0;
            }
        }
        public int this[int n] {
            get {
                return mebius[n];
            }
        } 
    }
}
```

```csharp:PrimeNumber.cs
using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;

namespace Gushwell.Etude {
    public static class PrimeNumber {

        public static IEnumerable<int> Enumerate() {
            // 2,3は既知の素数とする
            var primes = new List<int>() { 2, 3 };
            foreach (var p in primes)
                yield return p;

            // 4以上の整数から素数を列挙する。int.MaxValueを超えたときには対処していない
            int ix = 0;
            while (true) {
                int prime1st = primes[ix];
                int prime2nd = primes[++ix];
                // ふるい用の配列の下限、上限を求め、配列を確保する。
                var lower = prime1st * prime1st;
                var upper = prime2nd * prime2nd - 1;
                // ふるいは、[4:8], [9:24], [25:48], [49:120]... と変化する。
                // []内の数値は、配列の下限と上限
                var sieve = new BoundedBoolArray(lower, upper);

                // 求まっている素数を使い、ふるいに掛ける 
                foreach (var prime in primes.Take(ix)) {
                    var start = (int)Math.Ceiling((double)lower / prime) * prime;
                    for (int index = start; index <= upper; index += prime)
                        sieve[index] = true;
                }

                // ふるいに掛けられて残った値が素数。これを列挙する。
                // 併せて、求まった素数は、primesリストに記憶していく。
                // この素数が次にふるいに掛ける際に利用される。
                for (int i = lower; i <= upper; i++) {
                    if (sieve[i] == false) {
                        primes.Add(i);
                        yield return i;
                    }
                }
            }
        }
    }

    // 下限、上限が指定できるbool型配列
    class BoundedBoolArray {
        private BitArray _array;
        private int _lower;

        public BoundedBoolArray(int lower, int upper) {
            _array = new BitArray(upper - lower + 1);
            _lower = lower;
        }

        public bool this[int index] {
            get {
                return _array[index - _lower];
            }
            set {
                _array[index - _lower] = value;
            }
        }
    }
}

```

##簡単な解説

まず、Mebius関数の値を格納する配列を用意し、1で初期化。 
それから、エラトステネスのふるいと似たやり方で、

<pre>
2,4,6,8,10...
3,6,9,12,15...
5,15,20,25,30...
...
</pre>

の位置の要素に、-1を掛けていきます。これで、

　n が相異なる偶数個の素数の積ならば μ(n) = 1
　n が相異なる奇数個の素数の積ならば μ(n) = -1  

の部分を求めます。

さらに、

<pre>
4,8,12,16,20...
9,18,27,36...  
25,50,75...
...
</pre>

の要素（素数の２乗で割りきれる数）には、0を代入して、

　n が平方因子を持つ（平方数で割り切れる）ときμ(n) = 0  

を求めます。

2以上の数で、上記どれにも該当しなかいのはあり得ないので、配列のすべてにメビウス関数の値が設定されることになります。

このC#のコードでは、 μ(1), μ(2) μ(3) ... の値を列挙するのはなく、整数nを与えると、メビウス関数の値が戻ってくるような、読み取り専用のインデクサを定義しました。
範囲チェックはやってません。


```csharp
var mebius = new Mebius(100);
int n = mebius[10];
```
のようにして値を得ます。

なお、素数を求めるPrimesクラスは、[ここ](http://qiita.com/gushwell/items/6370e32371e983358519)で示したメソッドを元に、クラスとして定義しなおしたものです。


## 実行結果の一部

実行結果の一部を掲載します。

    μ(1) = 1
    μ(2) = -1
    μ(3) = -1
    μ(4) = 0
    μ(5) = -1
    μ(6) = 1
    μ(7) = -1
    μ(8) = 0
    μ(9) = 0
    μ(10) = 1
    μ(11) = -1
    μ(12) = 0
    μ(13) = -1
    μ(14) = 1
    μ(15) = 1
    μ(16) = 0
    μ(17) = -1
    μ(18) = 0
    μ(19) = -1
    μ(20) = 0
    μ(21) = 1

------------------
この記事は、Gushwell's C# Programming Pageで公開したものを加筆・修正したものです。