さんぽみち

なにか思いついたときとかに気まぐれに更新されます。

C#で Project Euler P3

はじめに

 P2はフィボナッチ数列級数に関する問題だけど、問題に魅力を感じなかったので飛ばしました。
 で、僕としては好みなP3をやります。

 そういや問題には直接関係ないけど、C#って末尾再帰の最適化してくれないのね。。。
 再帰的に書こうと思ってたけど再帰はなるべく減らしてます。

問題

Problem 3 「最大の素因数」
13195 の素因数は 5, 7, 13, 29 である.
600851475143 の素因数のうち最大のものを求めよ.
Problem 3 - PukiWiki

解答

P3.cs

using System;

public static class P3
{
    public static void Main()
    {
        Console.WriteLine(Calc());
    }
    
    private static ulong Calc()
    {
        ulong n = 600851475143;
        
        ulong index;
        for(index = 0; n != 1; index++)
        {
            while(n % PrimeGenerator.GetPrime(index) == 0)
            {
                n /= PrimeGenerator.GetPrime(index);
            }
        }
        
        return PrimeGenerator.GetPrime(index - 1);
    }
}

PrimeGenerator

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

public static class PrimeGenerator
{
    private static ulong primeIndex;    // 素数の添え字
    private static ulong primeCounter;  // 素数を探すために回すカウンタ
    
    private static readonly Predicate<ulong> TryDivision;   // 試し割りを行うデリゲート
    
    private static Dictionary<ulong, ulong> primeList = new Dictionary<ulong, ulong>(); // 見つけた素数をバッファする空間
    
    // 初期化
    static PrimeGenerator()
    {
        primeList.Add(0, 2);    // 2の登録
        primeIndex = 1;
        primeCounter = 1;
        ulong multi = 1;        // 素数の積、偶数はとばす為2は登録しない
        while(primeIndex < 10)
        {
            while(MyMath.EuclideanAlgorithem(primeCounter += 2, multi) != 1){ }
            primeList.Add(primeIndex++, primeCounter);
            multi *= primeCounter;
        }
        
        // 試し割りの登録、素数の積との最大公約数が1であるかで判定する
        TryDivision = n => MyMath.EuclideanAlgorithem(multi, n) == 1;
    }
    
    // 指定した添え字の素数を取得する
    public static ulong GetPrime(ulong idx)
    {
        while(primeIndex <= idx)
        {
            while(!PrimeTest(primeCounter += 2)){ }
            primeList.Add(primeIndex++, primeCounter);
        }
        return primeList[idx];
    }
    
    // 素数判定を行う
    private static bool PrimeTest(ulong n)
    {
        // 試し割り、O(log n)
        if(!TryDivision(n))
        {
            return false;
        }
        // フェルマーテスト、O(log n)
        if(MyMath.ModPow(2, n - 1, n) != 1)
        {
            return false;
        }
        // ミラーラビンテスト、O(m * (log n)^2)
        if(!MillerRabinTest(n))
        {
            return false;
        }
        // 全てのテストに合格
        return true;
    }
    
    private const int BigIntegersTestCount = 15;   // 大きな素数判定の時は適当に15回くらいテストする
    private static readonly ulong[] MRTestPattern4759123141 = new ulong[]{ 2, 7 , 61 };
    private static readonly ulong[] MRTestPattern341550071728321 = new ulong[]{ 2, 3, 5, 7, 11, 13, 17 };
    // ミラーラビンテスト
    private static bool MillerRabinTest(ulong n)
    {
        // テストパターンの選択
        IEnumerable<ulong> testPattern = MRTestPattern4759123141;
        if(n < 4759123141)
        {
            // 初期化で代入した為処理なし
        }
        else if(n < 341550071728321)
        {
            testPattern = MRTestPattern341550071728321;
        }
        else
        {
            // これより大きい数字のテストパターンを知らないので
            // 適当な素数を乱択する
            testPattern = GetRandomTestPattern().Take(BigIntegersTestCount);
        }
        
        // 全てのパターンでテストする
        foreach(var testNum in testPattern)
        {
            if(testNum >= n)
            {
                break;
            }
            if(!MillerRabinTestSub(n, testNum))
            {
                return false;
            }
        }
        return true;
    }
    private static Random rnd = new Random();
    static List<ulong> usedPrimes = new List<ulong>();
    // 素数を乱択して返すシーケンス
    private static IEnumerable<ulong> GetRandomTestPattern()
    {
        usedPrimes.Clear();
        ulong index;
        while(true)
        {
            index = (ulong)rnd.Next(primeList.Count);
            if(usedPrimes.Contains(index))
            {
                continue;
            }
            usedPrimes.Add(index);
            yield return primeList[index];
        }
    }
    // ミラーラビンテストを実際に行う
    private static bool MillerRabinTestSub(ulong n, ulong a)
    {
        ulong m1 = n - 1;           // n - 1
        // n - 1 を r * d へ分解、ただし r = 2^s, d ≡ 1 (mod 2)
        ulong r = m1 & (~m1 + 1);   // 2のべき乗要素
        ulong d = m1 / r;           // 奇数要素
        
        ulong pow2 = 0;
        ulong mp = MyMath.ModPow(a, d, n);
        // 2の0乗で 1 または -1 の場合は合格
        if(mp == 1 || mp == m1)
        {
            return true;
        }
        for(pow2 = 1; pow2 <= r; pow2 <<= 1)
        {
            mp = MyMath.ModPow(a, pow2 * d, n);
            // 1 より先に -1 が出た場合は合格
            if(mp == m1)
            {
                return true;
            }
            // 1 が -1 より先に出た場合は不合格
            if(mp == 1)
            {
                return false;
            }
        }
        // 合格して抜けなければ不合格
        return false;
    }
    
    public static IEnumerable<ulong> GetPrimeList()
    {
        ulong index = 0;
        while(true)
        {
            yield return GetPrime(index);
            index++;
        }
    }
    public static IEnumerable<ulong> GetPrimeList(Func<ulong, ulong, bool> predicator)
    {
        for(ulong index = 0; predicator(index, GetPrime(index)); index++)
        {
            yield return GetPrime(index);
        }
    }
}

MyMath.cs

public static class MyMath
{
    // O(log n)
    public static ulong ModPow(ulong n, ulong p, ulong m)
    {
        ulong res = 1;
        for(; p > 0; p >>= 1)
        {
            res = (p & 1) == 1 ? res * n % m : res;
            n = (n * n) % m;
        }
        return res;
    }
    // O(log n)
    public static ulong EuclideanAlgorithem(ulong a, ulong b)
    {
        ulong work;
        while(b != 0)
        {
            work = b;
            b = a % b;
            a = work;
        }
        return a;
    }
}

計算量 $O(m\varphi(n)(\log n)^2)$
ただし、m はミラーラビンテストのテストパターン数、φはオイラーのトーシェント関数。
計算量あってるか怪しいけど多分これくらい!
計算量でφ関数とか見たことないけどこうとしか書けなかった。
最大の素因数を求めるということで、素直に素因数分解した。
総当たりの計算量$O(\sqrt{n})$くらいになるのかな?に比べればφ関数で母数をぐっと減らせるから速いと思う。

↑違いました、素因数でない素数の判定も行っているので$O(\sqrt{n}m(\log n)^2)$。
下限を出そうと思ったけどある自然数以下の素数を数える方法がわからずに表現できませんでした。。。

でも、使う素数の範囲があらかじめわかるので、素数判定するよりもエラトステネスの篩で素数表を作った方が速いかも。