C#:「ミラー・ラビン素数判定法」による素数判定メソッド

はじめに

「C#:エラトステネスの篩(ふるい)再び」で素数を列挙するコードを示しましたが、ある正の整数が素数かどうかを知りたいという時には、ちょっと使い勝手が悪いです。

ということで、素数かどうかを調べるメソッドも書いてみました。
さすがに、2,3,5,7,9,11,13,15,17と試し割りするのでは能が無さすぎるので、調べてみたところた、ミラー・ラビン素数判定法というものがあることを知りました。

最初のバージョン(ミラー・ラビン素数判定法)

いくつかのサイトを参考に書いたのが以下のクラスです。

IsPrimeMillarRrabin.cs

using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;

namespace Gushwell.Etude {
    public static class PrimeNumber {

        static long[] seedPrimes = {
          /*1,2,3,4, 5, 6, 7  8, 9,10,11,12,13,14,15,*/
            2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
        };

        public static bool IsPrimeMillarRrabin(long num) {
            if (num <= 1)
                return false;
            if ((num & 1) == 0)
                return num == 2;

            if (num < 100 && seedPrimes.Contains((int)num))
                return true;

            var WitnessMax = GetWitnessMax(num);

            long d = num - 1;
            long s = 0;
            while ((d & 1) == 0) {
                s++;
                d >>= 1;
            }
            foreach (var w in seedPrimes.Take(WitnessMax)) {
                if (!MillarRrabin(num, s, d, w))
                    return false;
            }
            return true;
        }

        private static int GetWitnessMax(long num) {
            if (num < 2047)
                return 1;
            if (num < 1373653)
                return 2;
            if (num < 25326001)
                return 3;
            if (num < 3215031751)
                return 4;
            if (num < 2152302898747)
                return 5;
            if (num < 3474749660383)
                return 6;
            if (num < 341550071728321)
                return 7;
            if (num < 3825123056546413051)
                return 9;
            return 12;
        }

        private static bool MillarRrabin(long num, long s, long d, long witness) {
            long x = ModPow(witness, d, num);
            if (x == 1)
                return true;
            for (long r = 0; r < s; r++) {
                if (x == num - 1)
                    return true;
                BigInteger rem;
                BigInteger.DivRem(BigInteger.Multiply(x,x), num, out rem);
                x = (long)(rem);
            }
            return false;
        }

        private static long ModPow(long baseValue, long exponent, long modulus) {
            return (long)BigInteger.ModPow(baseValue, exponent, modulus);
        }

    }
}

ModPowの底は、2,3,5,7,11,13,17,19,23,29,31,37まで調べれば、2^64 までの素数が決定的素数判定として利用できるようです。

ところどころ、BigInteger使ってますが、オーバーフローを回避するためと、long型のSqrtメソッドが用意されていないための苦肉の策です。もっと良い方法もあると思いますが、良い方法が思い浮かびません...

力技の方法と速度を比較してみる

このIsPrimeMillarRrabinメソッドの検証のために、以下のようなメソッドも書きました。
奇数すべてを試し割りするのではなく、できるだけ試し割りの数を減らす工夫をして速度を上げています。

IsPrimeBruteforce.cs

public static bool IsPrimeBruteforce(long num) {
    if (num == 1)
        return false;
    if (num != 2 && num % 2 == 0)
        return false;
    if (num != 3 && num % 3 == 0)
        return false;
    if (num != 5 && num % 5 == 0)
        return false;
    long i = 0;
    while (true) {
        foreach (var p in seedPrimes.Skip(3).Take(8)) {
            // 30m+2, 30m+3, 30m+4, 30m+5, 30m+6、30m+8、30m+9、30m+12... は割る必要はない。
            var primeCandidte = p + i;
            if (primeCandidte > Math.Sqrt(num))
                return true;
            if (num % (primeCandidte) == 0)
                return false;
        }
        i += 30;
    }
}

ふたつのメソッドの結果を比較するコードで確認したところ、問題なく動いているようです。
それでは、IsPrimeMillarRrabinメソッドは本当に速いんだろうかと調べてみました。

1-100万までの整数で素数がいくつあるかをカウントするコードを書いて調べたところ、あれ、IsPrimeBruteforceの方が速いです。

メソッド	結果	速度
IsPrimeBruteforce	78498	671 ミリ秒
IsPrimeMillarRrabin	78498	921 ミリ秒

どこかで速さが逆転するはずと思い調べたら、だいたい2,000,000くらいの値の素数判定で速度が逆転するようです。

ちなみに、50,000,000,000から50,000,100,000までの整数で素数がいくつあるかをカウントするコードを書いて調べたところ、僕のPCでは、以下のような結果で圧倒的にミラー・ラビン素数判定法の方が速いという結果がでました。

メソッド	結果	速度
IsPrimeBruteforce	4097	115930 ミリ秒
IsPrimeMillarRrabin	4097	665 ミリ秒

ということで、最終的には、IsPrimeメソッドは、以下のようなハイブリッド型としました。

最終版IsPrimeメソッド　(ハイブリッド版)

IsPrime.cs

public static class PrimeNumber {

    static long[] seedPrimes = {
      /*1,2,3,4, 5, 6, 7  8, 9,10,11,12,13,14,15,*/
        2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
    };

    public static bool IsPrime(long num) {
        if (num == 1)
            return false;
        if (seedPrimes.Contains(num))
            return true;
        if (seedPrimes.Any(x => num % x == 0))
            return false;

        return (num < 2000000) ? IsPrimeBruteforce(num) : IsPrimeMillarRrabin(num);
    }

    private static bool IsPrimeBruteforce(long num) {
        if (num == 1)
            return false;
        if (num != 2 && num % 2 == 0)
            return false;
        if (num != 3 && num % 3 == 0)
            return false;
        if (num != 5 && num % 5 == 0)
            return false;
        long i = 0;
        while (true) {
            foreach (var p in seedPrimes.Skip(3).Take(8)) {
                // 30m+2, 30m+3, 30m+4, 30m+5, 30m+6、30m+8、30m+9、30m+12... は割る必要はない。
                var primeCandidte = p + i;
                if (primeCandidte > Math.Sqrt(num))
                    return true;
                if (num % (primeCandidte) == 0)
                    return false;
            }
            i += 30;
        }
    }

    private static bool IsPrimeMillarRrabin(long num) {
        if (num <= 1)
            return false;
        if ((num & 1) == 0)
            return num == 2;

        if (num < 100 && seedPrimes.Contains((int)num))
            return true;

        var WitnessMax = GetWitnessMax(num);

        long d = num - 1;
        long s = 0;
        while ((d & 1) == 0) {
            s++;
            d >>= 1;
        }
        foreach (var w in seedPrimes.Take(WitnessMax)) {
            if (!MillarRrabin(num, s, d, w))
                return false;
        }
        return true;
    }


    private static int GetWitnessMax(long num) {
        if (num < 2047)
            return 1;
        if (num < 1373653)
            return 2;
        if (num < 25326001)
            return 3;
        if (num < 3215031751)
            return 4;
        if (num < 2152302898747)
            return 5;
        if (num < 3474749660383)
            return 6;
        if (num < 341550071728321)
            return 7;
        if (num < 3825123056546413051)
            return 9;
        return 12;
    }


    private static bool MillarRrabin(long num, long s, long d, long witness) {
        long x = ModPow(witness, d, num);
        if (x == 1)
            return true;
        for (long r = 0; r < s; r++) {
            if (x == num - 1)
                return true;
            BigInteger rem;
            BigInteger.DivRem(BigInteger.Multiply(x, x), num, out rem);
            x = (long)(rem);
        }
        return false;
    }

    private static long ModPow(long baseValue, long exponent, long modulus) {
        return (long)BigInteger.ModPow(baseValue, exponent, modulus);
    }

}

IsPrimeの前半でも素数判定の一部を行っているので、IsPrimeBruteforce、IsPrimeMillarRrabinでは一部不要な判断をしているのですが、IsPrimeBruteforce、IsPrimeMillarRrabinは、それ単独で素数判定が成り立つようにしたかったので、そのままとしています。