素数判断和因子分解

这里记录的是素数判断 Miller Rabbin 的快速判定算法，而分解使用的算法是 Pollard Rho 的玄学分解算法，这个算法的理论依据是生日悖论，这里不具体讨论这个数学问题，具体讨论看前面链接就可以了。

以下我们分别简要介绍这两个算法的原理

首先这是Miller Rabbin的维基百科介绍页，由于页面是中文，这里不做重复介绍。

下面简要解释一下Pollard’s Rho算法的原理，这是它的维基百科页面

首先我们来定义这样一个函数

$$f(x)=(x*x+c)%n$$

其中n是需要被分解的数，c是一个任意小于n的正整数常数。那么我们给一个初始值x1，通过计算 $f(x1)$ 得到x2，然后去求 $|x1-x2|$ 与 n 的最大公约数，只要结果大于1就找到了分解结果。那如果最大公约数仍然是1，那么就通过刚才的x2求出 $x3=f(x2)$ 并计算 $|x1-x3|$ 与 n 的最大公约数。也就是说，我们可以通过这个公式得到一个数列

$$ x1,x2,x3,…,xn,… $$

但这个数列总是不会无限的，总是会循环的，那怎么去判断它有没有陷入循环呢，这用到的办法是Floyd判环算法，两个数在数列里以不同速度前进，如果环存在，那它们总会有相遇的时刻。在代码里，我们令变量x是正常速度，变量y是在x每迭代两次后y再迭代一次，所以如果环的长度是k，那么迭代k次就能发现环的存在。

而在具体实现中，如果发现环的存在，我们需要变换常数c，再次进行尝试。那么我们下面不废话直接上代码，作为笔记。为了避免部分环境的输入输出问题，main函数使用了C++的输入输出，其它部分可以作为C函数使用。

#include<iostream>
#include<cstdint>

int64_t gcd(int64_t a, int64_t b)
{
    return b == 0 ? a : gcd(b, a%b);
}

int64_t mul(int64_t a, int64_t b, int64_t mod)
{
    a %= mod;
    int64_t res = 0;
    while (b) {
        if (b & 1) res = (res + a) % mod;
        b >>= 1;
        a = (a << 1) % mod;
    }
    return res;
}

int64_t quick_pow(int64_t a, int64_t n, int64_t mod)
{
    int64_t res = 1;
    while (n)
    {
        if (n & 1) res = mul(res, a, mod);
        a = mul(a, a, mod);
        n >>= 1;
    }
    return res;
}

int miller_rabbin(int64_t n, int a)
{
    int r = 0;
    int64_t s = n - 1;
    if (!(n % a))
        return 0;
    while (!(s & 1))
    {
        s >>= 1;
        r++;
    }
    int64_t k = quick_pow(a, s, n);
    if (k == 1)
        return 1;
    for (int j = 0; j < r; j++, k = mul(k, k, n))
        if (k == n - 1)
            return 1;
    return 0;
}

int is_prime(int64_t n)
{
    int test_large[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };
    int test_small[] = { 2, 7, 61 };
    int test_size = 9, *test = test_large;
    if (n < 2)
        return 0;
    if (n < 4759123141)
    {
        test = test_small;
        test_size = 3;
    }
    for (int i = 0; i < test_size; ++i)
    {
        if (n == test[i])
            return 1;
        if (miller_rabbin(n, test[i]) == 0)
            return 0;
    }
    return 1;
}

int64_t pollard_rho(int64_t n, int c)
{
    int64_t x = rand() % (n - 1) + 1, y;
    y = x = (mul(x, x, n) + c) % n;
    for (int n_try = 3, i = 0; n_try; )
    {
        x = (mul(x, x, n) + c) % n;
        if (y == x)
        {
            --n_try;
            c = rand() % (n - 1) + 1;
            continue;
        }
        int64_t d = gcd((y - x + n) % n, n);
        if (d != 1 && d != n)
            return d;
        if (++i == 2)
        {
            y = (mul(y, y, n) + c) % n;
            i = 0;
        }
    }
    return n;
}

int main(void)
{
    int64_t n;
    std::cout << "质数计算器，请输入任意小于1e18的正整数" << std::endl;
    while (std::cin >> n)
    {
        if (is_prime(n))
        {
            std::cout << n << "是质数" << std::endl;
        }
        else if (n > 2)
        {
            int64_t p = pollard_rho(n, rand() % (n - 1) + 1);
            std::cout << n << "不是质数，能被" << p << "整除" << std::endl;
        }
        else
        {
            std::cout << "请输入大于1小于1e18的整数" << std::endl;
        }
    }
    return 0;
}

以下为python版本实现

import random

def _try_composite(a, d, n, s):
    if pow(a, d, n) == 1:
        return False
    for i in range(s):
        if pow(a, 2**i * d, n) == n-1:
            return False
    return True

def _is_prime(n):
    if n in _known_primes:
        return True
    if n < 2 or any((n % p) == 0 for p in _known_primes):
        return False
    return True

def is_prime(n, _precision_for_huge_n = 32):
    if n in _known_primes:
        return True
    d, s = n - 1, 0
    while not d % 2:
        d, s = d >> 1, s + 1
    if n < 4759123141: 
        return not any(_try_composite(a, d, n, s) for a in (2, 7, 61))
    if n < 1122004669633: 
        return not any(_try_composite(a, d, n, s) for a in (2, 13, 23, 1662803))
    if n < 3317044064679887385961981: 
        return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41))
    return not any(_try_composite(a, d, n, s) 
                   for a in _test_primes[:_precision_for_huge_n])
 
_known_primes = [2, 3]
_known_primes += [x for x in range(5, 1000, 2) if _is_prime(x)]
_test_primes = [2,3,307] + _known_primes[2:]

def gcd(a, b):
    if b == 0:
        return a
    return gcd(b, a%b)

def pollard_rho(n):
    x, c, 0 = random.randint(1, n - 1), random.randint(1, n - 1), 0
    y = x = (pow(x, 2, n) + c) % n
    for n_try in range(3):
        while True:
            x = (pow(x, 2, n) + c) % n
            if x == y:
                c = random.randint(1, n - 1)
                break
            d = gcd(y - x + n, n)
            if d > 1 and d < n:
                return d
            a += 1
            if a == 2:
                y = (pow(y, 2, n) + c) % n
                a = 0
    return n

number = 3451973391686190983
for i in range(number, number + 100, 2):
    if is_prime(i):
        print(i)
        break

if is_prime(number):
    print(number, "prime number")
else:
    a = pollard_rho(number)
    print(number, "=", a, "*", number // a)