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

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

More than 1 year has passed since last update.

はじめに

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は、それ単独で素数判定が成り立つようにしたかったので、そのままとしています。

gushwell
株式会社ジード / Microsoft MVP for Developer Technologies 2005-2019 / 著書『実戦で役立つ C#プログラミングのイディオム/定石&パターン』『新・標準プログラマーズライブラリ なるほどなっとく C#入門』『C#プログラミング入門―オブジェクト指向のプログラミング手法を基礎から解説』
https://github.com/gushwell
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
ユーザーは見つかりませんでした