Edited at

C#:メビウス関数

More than 1 year has passed since last update.

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