このプログラムの原型を書いたのは随分と昔のことです。当時素数を勉強したいと思って、このメビウス関数のコードを書いてみたのですが、数学の素養が欠けている僕には、このメビウス関数が一体何なのかが、イマイチわからない...でもせっかく書いたのだからと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までのメビウス関数の値を求めるコードを書いてみました。
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];
}
}
}
}
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...
...
の要素（素数の２乗で割りきれる数）には、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で公開したものを加筆・修正したものです。