Pythonで素因数分解をする（競プロ用）

# 素因数分解完全に理解した
嘘です先達の皆様マサカリください:bow:

素因数分解アルゴリズムを、競技プログラミング用に実装しました。

Pythonで競プロは素人？まあそれはそれとして。

### 素因数分解
素因数分解は一般に難しい計算であり、この難しさが現代の暗号の安全性を支えていたりします。

素因数分解には様々なアルゴリズムがあり、とても書ききれない（理解しきれない）ので、下記リンク参考。

[素因数分解アルゴリズム(特にSQUFOF)のこと](http://d.hatena.ne.jp/lemniscus/20130226/1361874593)
[素数判定いろいろ - フェルマーテスト・ミラーラビン素数判定法](https://qiita.com/srtk86/items/609737d50c9ef5f5dc59)


# 実装
実装に当たって、下記のQiita記事を参考にしました。
[C#:「ミラー・ラビン素数判定法」による素数判定メソッド](https://qiita.com/gushwell/items/ff9ed83ba55350aaa369)
[C#:ポラード・ロー素因数分解法](https://qiita.com/gushwell/items/561afde2e00bf3380c98)

コードで語れと言わんばかりにコードしか書いてないですが、ちゃんとした解説は参考記事を読んでください。

ごく簡単に言えば、
1. `find_factors`が`find_factor`を繰り返し適用し、`n`の全ての素因数を見つける。
2. `find_factor`は、再帰的＆確率的に`n`の素因数を一つ見つける。
3. `is_prime`は、`find_factor`で素数判定に使われる。
4. `is_prime_miller_rabin`は素数判定を確率的に行える手法（ミラーラビン素数判定法）を繰り返し、精度が十分に高くなるようにする。
 
という感じです。

```python
import math
import random


class Prime:
    seed_primes = [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]

    def is_prime(self, n):
        """
        prime test (hybrid)

        see also: https://qiita.com/gushwell/items/ff9ed83ba55350aaa369

        :param n:
        :return: boolean
        """
        is_prime_common = self.is_prime_common(n)
        if is_prime_common is not None:
            return is_prime_common

        if n < 2000000:
            return self.is_prime_brute_force(n)
        else:
            return self.is_prime_miller_rabin(n)

    def is_prime_common(self, n):
        if n == 1: return False
        if n in Prime.seed_primes: return True
        if any(map(lambda x: n % x == 0, self.seed_primes)): return False

    def is_prime_brute_force(self, n):
        """
        brute force prime test
        use with is_prime_common if you want to skip seed_primes

        :param n:
        :return: boolean
        """
        for k in range(2, int(math.sqrt(n)) + 1):
            if n % k == 0:
                return False
        return True

    def is_prime_miller_rabin(self, n):
        """
        miller rabin prime test
        use with is_prime_common if you want to skip seed_primes

        see also
            algorithm: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
            implementation: https://qiita.com/srtk86/items/609737d50c9ef5f5dc59
            improvement: https://qiita.com/gushwell/items/ff9ed83ba55350aaa369

        :param n:
        :return: boolean
        """

        d = n - 1
        while d & 1 == 0:
            d >>= 1

        # use one of these lines / upper is more efficient.
        witnesses = self.get_witnesses(n)
        # witnesses = [random.randint(1, n - 1) for _ in range(100)]

        for w in witnesses:
            y = pow(w, d, n)

            while d != n - 1 and y != 1 and y != n - 1:
                y = (y * y) % n
                d <<= 1

            if y != n - 1 and d & 1 == 0:
                return False

        return True

    def get_witnesses(self, num):
        def _get_range(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

        return self.seed_primes[:_get_range(num)]

    def gcd(self, a, b):
        if a < b:
            return self.gcd(b, a)
        if b == 0:
            return a
        while b:
            a, b = b, a % b
        return a

    @staticmethod
    def f(x, n, seed):
        """
        pseudo prime generator
        :param x:
        :param n:
        :param seed:
        :return: pseudo prime
        """
        p = Prime.seed_primes[seed % len(Prime.seed_primes)]
        return (p * x + seed) % n

    def find_factor(self, n, seed=1):
        """
        find one of factor of n
        this function is based to Pollard's rho algorithm

        see also
            algorithm: https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
            implementation: https://qiita.com/gushwell/items/561afde2e00bf3380c98

        :param n:
        :param seed:
        :return: factor
        """
        if self.is_prime(n):
            return n

        x, y, d = 2, 2, 1
        count = 0
        while d == 1:
            count += 1
            x = self.f(x, n, seed)
            y = self.f(self.f(y, n, seed), n, seed)
            d = self.gcd(abs(x - y), n)

        if d == n:
            return self.find_factor(n, seed+1)
        return self.find_factor(d)

    def find_factors(self, n):
        primes = {}
        if self.is_prime(n):
            primes[n] = 1
            return primes

        while n > 1:
            factor = self.find_factor(n)

            primes.setdefault(factor, 0)
            primes[factor] += 1

            n //= factor

        return primes


prime = Prime()

if __name__ == '__main__':
    prime = Prime()
    assert prime.find_factors(36610051291281) == {653: 1, 593783: 1, 3: 3, 13: 1, 269: 1}
```