• 2
    いいね
  • 0
    コメント
この記事は最終更新日から1年以上が経過しています。

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

メビウス関数

メビウス関数の定義は、以下のようなものです。

1) μ(1) = 1
2) μ(n) = 0 (n が平方因子を持つ(平方数で割り切れる)とき)
3) μ(n) = (-1)**k (n が相異なる k 個の素因数に分解されるとき)
      n が相異なる偶数個の素数の積ならば μ(n) = 1
      n が相異なる奇数個の素数の積ならば μ(n) = -1
      **は累乗を表す。

詳しくは、Wikipedia(メビウス関数)を参照してください。

C#のコード

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

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];
            }
        } 
    }
}
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で初期化。
それから、エラトステネスのふるいと似たやり方で、


2,4,6,8,10...
3,6,9,12,15...
5,15,20,25,30...
...

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

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

の部分を求めます。

さらに、


4,8,12,16,20...
9,18,27,36...  
25,50,75...
...

の要素(素数の2乗で割りきれる数)には、0を代入して、

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

を求めます。

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

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

var mebius = new Mebius(100);
int n = mebius[10];

のようにして値を得ます。

なお、素数を求めるPrimesクラスは、ここで示したメソッドを元に、クラスとして定義しなおしたものです。

実行結果の一部

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

μ(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で公開したものを加筆・修正したものです。