C#:ポラード・ロー素因数分解法



## 力ずくで素因数分解

まずは、与えられた数(long型)を力ずくで強引に素因数分解するプログラムを書いてみました。

### C#のコード (Brute force版)

```csharp:PrimeFactorBruteForce.cs
public class PrimeFactorBruteForce {
    private List<long> _factors = new List<long>();
    public IEnumerable<long> Enumerate(long num) {
        _factors.Clear();
        if (num == 1)
            return new long[] { 1 };
        var m = num;
        long n = 2;
        while (m > 1) {
            // m を nで割り続ける。その結果を返す。
            // かつ、割り切れたら、 _factorsにAdd
            while (m % n == 0) {
                _factors.Add(n);
                m = m / n;
            }
            n++;
        }
        return _factors;
    }
}
```

Enumerateを呼び出すと、素因数分解した結果である素因数が列挙できます。

しかし、このコード、大きな素数を与えた場合、めちゃくちゃ遅いです。
そこで、少し改良したのが次のコード

### C#のコード (Brute force改良版)

```csharp:PrimeFactorBruteForce2.cs
public class PrimeFactorBruteForce {
    private List<long> _factors = new List<long>();
    public IEnumerable<long> Enumerate(long num) {
        _factors.Clear();
        if (num == 1)
            return new long[] { 1 };
        var n = DivideWith(num, 2);
        n = DivideWith(n, 3);
        n = DivideWith(n, 5);
        long i = 0;
        while (n > 1) {
            int[] ps = { 7, 11, 13, 17, 19, 23, 29, 31 };
            foreach (var p in ps) {
                // 30m+2, 30m+3, 30m+4, 30m+5, 30m+6、30m+8... は割る必要はない。
                n = DivideWith(n, i + p);
                if (n == 1)
                    break;
            }
            i += 30;
        }
        return _factors;
    }

    // m を nで割り続ける。割り切れたら、 _factorsにAdd。
    private long DivideWith(long m, long n) {
        while (m % n == 0) {
            _factors.Add(n);
            m = m / n;
        }
        return m;
    }
}
```

すこし工夫して、単に、2,3,4,5,6,7...と順番に割り算をするよりは効率よく素因数分解するようにしています。

まず、31までの素数で割り、それ以降は、30m+7, 30m+11, 30m+13, 30m+17, 30m+19, 30m+23, 30m+29, 30m+31 の数で割っています。mは整数です。

30m+2, 30m+3, ..., 30m+8, 30m+9, 30m+10, 30m+12 などは、素数でないことが明らかなので対象から除外しています。そのため、素数でない数でも割ることにはなりますが、単純なコードよりもすこし効率良く素因数分解できます。

ただ、大きな素数が含まれる数に対しては、無意味な試し割りが何べんも行われるため、ほんのわずかに速度向上しただけでした。

ちなみにこのロジックは、[ミラー・ラビン素数判定法による素数判定メソッド](http://qiita.com/gushwell/items/ff9ed83ba55350aaa369)の中で示した力技のメソッドIsPrimeBruteforceと同じ考えで書いたコードです。

## ポラード・ロー素因数分解法

もっと、速くする方法は無いのかと思って調べたら、ポラード・ロー素因数分解法というのがあることを知りました。

そこで、wikipediaの[ポラード・ロー素因数分解法](https://ja.wikipedia.org/wiki/%E3%83%9D%E3%83%A9%E3%83%BC%E3%83%89%E3%83%BB%E3%83%AD%E3%83%BC%E7%B4%A0%E5%9B%A0%E6%95%B0%E5%88%86%E8%A7%A3%E6%B3%95) のページに書いてある通りに書いてみたのですが、どうも結果が正しくありません。

良く良く読むと、与える値は素数以外が大前提で、素数以外を与えても失敗する(素数でない因数を見つけてしまう)場合があると書いてあります。

ということで、そのあたりを改善したコードが次のコードです。
前の２つの方法に比べると、圧倒的な速さです。


```csharp:PrimeFactor.cs
public class PrimeFactor {

    public IEnumerable<long> Enumerate(long n) {
        while (n > 1) {
            long factor = GetFactor(n);
            yield return factor;
            n = n / factor;
        }
    }

    private long GetFactor(long n, int seed = 1) {
        if (n % 2 == 0)
            return 2;
        if (IsPrime(n))
            return n;
        long x = 2;
        long y = 2;
        long d = 1;
        long count = 0;
        while (d == 1) {
            count++;
            x = f(x, n, seed);
            y = f(f(y, n, seed), n, seed);
            d = Gcd(Math.Abs(x - y), n);
        }
        if (d == n)
            // 見つからなかった、乱数発生のシードを変えて再挑戦。
            return GetFactor(n, seed+1);
        // 素数でない可能性もあるので、再度呼び出す
        return GetFactor(d);
    }

    private int[] seeds = new int[] { 3, 5, 7, 11, 13, 17 };
    private long f(long x, long n, int seed) {
        return (seeds[seed % 6] * x + seed) % n;
    }

    private static long Gcd(long a, long b) {
        if (a < b)
            return Gcd(b, a);  // 引数を入替えて自分を呼び出す
        if (b == 0)
            return a;
        long d = 0;
        do {
            d = a % b;
            a = b;
            b = d;
        } while (d != 0);
        return a;
    }

    // 効率は良くないが、これでも十分な速度がでたので、良しとする。
    private static bool IsPrime(long number) {
        long boundary = (long)Math.Floor(Math.Sqrt(number));

        if (number == 1)
            return false;
        if (number == 2)
            return true;

        for (long i = 2; i <= boundary; ++i) {
            if (number % i == 0)
                return false;
        }

        return true;
    }
}
```    

IsPrimeメソッドを、[ミラー・ラビン素数判定法による素数判定メソッド](http://qiita.com/gushwell/items/ff9ed83ba55350aaa369)で示した最終版のメソッドに置き換えれば、もっと速くなると思います。


使い方はこんな感じ

```csharp
var primeFactor = new PrimeFactor();
var line = Console.ReadLine();
var num = int.Parse(line);
var factors = primeFactor.Enumerate(num);
Console.Write(string.Join(" * ", factors));
```

## 実行例

以下、実行例です。

```
876402
2 * 3 * 3 * 181 * 269
```

```
1654620
2 * 2 * 3 * 5 * 11 * 23 * 109
```
-------------------
この記事は、Gushwell's C# Programming Pageで公開したものを加筆・修正したものです。