Help us understand the problem. What is going on with this article?

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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした