Miller-Rabin素性测试|Pollard's Rho算法

Miller-Rabin 素性测试

Miller-Rabin 素数测试

一本通上的M-R不透彻啊~

Miller-Rabin是利用随机化算法判断一个数是合数还是素数。

首先，如果一个数N是素数，那么他一定满足费马小定理。

\(a^{N-1}\equiv1\pmod N\)

我们可以任取数字a，计算这个式子的值来判断N是否为素数。

但是这么做不靠谱啊，有很多合数会被卡~

我们介绍一个相关的引理。

当p是素数且p大于2时，\(1\bmod p\)的平方根只有1和-1。

证明：

假设x是\(1\bmod p\)的平方根，则有\(x^2\equiv1\pmod p\)

则有\((x+1)(x-1)\equiv 0\pmod p\)

也就是说\((x+1)(x-1)\)能够整除素数p。根据欧几里得引理，(x+1)或(x-1)整除p，也就是\(x\equiv1\pmod p\)和\(x\equiv-1\pmod p\)

我们假设n是一个素数(n>2)，则n-1是一个偶数，且能被表示成\(2^s*d\)的形式，d是奇数。

那么对于任意的a和\(0\le r\le s-1\)，有\(a^d\equiv1\pmod n\)和\(a^{2^rd}\equiv -1\pmod n\)。

为什么呢？因为费马小定理。对于一个素数n，有\(a^{n-1}\equiv1\pmod n\)，由于上面的引理，如果我们不断对\(a^{n-1}\)取平方根，我们总会得到1或-1。如果我们得到了-1，则\(a^{2^rd}\equiv -1\pmod n\)成立，判定p是素数。如果我们从未得到-1，通过这个过程我们取遍了2的所有幂次，则\(a^d\equiv1\pmod n\)成立。

Miller-Rabin素性测试就在于上述原理的逆否，如果我们找到了一个a，使得对于任意的\(0\le r\le s-1\)，满足以下两个式子：\(a^d\not\equiv1\pmod n\)，\(a^{2^rd}\not\equiv-1\pmod n\)，n就不是一个素数。

如果一个数a通过了上述测试，就称a是n是合数的凭证。否则称a是n是素数的强伪证，也就是他可能出锅。（毕竟是随机化算法嘛）我们就需要多找几个a，使得被误判的概率更小。

Pollard's Rho算法

给出一个数N，快速提取N的[一个]因数。

生日悖论

在一群学生中，随机选一个学生，则他的生日为4 月 1 号的概率为\(\frac 1 {365}\)。

转化为古典概型，在[1,365]中随机选取一个数，该数为90的概率为\(\frac 1 {365}\)

然后我们现在需要随机选取2个人，他们两个人生日相同的概率也是\(\frac 1{365}\)

然而我们现在随机选取3个人，他们中任意两个人生日相同的概率是：

我们利用单步容斥，任意两人生日相同的概率=1-没有两人生日相同的概率。

没有两人生日相同的概率为\(\frac{365*364*363}{365*365*365}\)，约为0.9917958341152186（使用Windows7的计算器算出），则任意两人的生日相同概率为0.0082041658847814（同上）。

我们现在随机选取k个人，则容易得到的是，任意两个人生日相同的概率为：

\(1-\frac{365!}{(365-k)!*365^k}\)(k<=365)。

当k>365，P=1

当k=23时，P>50%，大概是\(\sqrt N\)

当k=42时，P>90%

当k=100是，P=99.99996928%

这个很大了~

因数分解

我们选取k个数，询问是否存在\(x_i-x_j\)能够整除N

当\(k=\sqrt N\)，可能性提高到了50%

我们将可能性从\(\frac 1 N\)提高到\(\sqrt{\frac1 N}\)。

对于一个\(10^{10}\)的大整数，我们只需要做\(10^5\)次比较，除法。这还是炸了。。。

我们选取k个数，询问是否存在\(gcd(x_i-x_y,N)>1\)。看着是棒了好多了，

下面我们介绍P-R

Pollard's Rho算法

这个Pollard's Rho算法只将两个数放在内存中。我们不一下子随机生成k个数，反而一个一个生成检查相邻的2个数。反复执行这个步骤。

我们有一个Interesting的函数\(f(x)=(x^2+a)\mod N\)。这个a可以自己指定，也可以rand()，这个不是重点。

我们从\(x_1=2\)开始，让\(x_2=f(x_1)\)....\(x_{i+1}=f(x_i)\)。如果相邻的两个数的差与N的gcd大于1，返回他们的gcd。作为N的因子。但是这个可能会出现环(详见模域平方)。

有一个Floyed发明的神奇的算法。小明和小红在操场上跑圈，小明的速度恰好是小红的2倍，最后小明跑2圈时，小红恰好跑了1圈，小明从小红的后面追上小红，但是实际上小明比小红快整整一圈。

我们用floyed判环，我们如果失败了我们只需要重新选择一个a即可。

这就是P R算法。

upd 9.6

然后洛谷板子题真心坑人各种卡常，瞬间变身卡同学

然后这是代码（抄了最优解Rank1好多。。。）

外加Miller-Rabin素数测试，有空更新一下

（这份代码细节不太透彻，留个坑，NOIP后重新透彻一遍）

#include <bits/stdc++.h>

using namespace std;

long long const_A = 5;

unsigned long long myrand()
{
    static unsigned long long seed = 3589425289U;
    return seed * 74658973 + 348633;
}

void read(long long &x)
{
    x = 0;
    static char ch;
    ch = getchar();
    while (!isdigit(ch))
        ch = getchar();
    while (isdigit(ch))
    {
        x = x * 10 + ch - 48;
        ch = getchar();
    }
}

inline long long qpow(long long x, long long y, long long p)
{
    long long ans = 1;
    while (y > 0)
    {
        if (y & 1)
            ans = ((__int128)ans * (__int128) x) % p;
        x = ((__int128)x * (__int128)x) % p;
        y >>= 1;
    }
    return ans;
}

bool Miller_Rabin(long long N, long long a)
{
    if (N == a)
    	return true;
    if (N == 46856248255981LL || N < 2 || N % a == 0)
        return false; 
    long long d = N - 1;
    while (!(d & 1))
        d >>= 1;
    long long t = qpow(a, d, N);
    while (d != N - 1 && t != 1 && t != N - 1)
    {
        t = ((__int128)t * (__int128)t) % N;
        d <<= 1;
    }
    return (t == N - 1) || ((d & 1) == 1);
}

inline bool prime(long long x)
{
    return (Miller_Rabin(x, 2) && Miller_Rabin(x, 3) && Miller_Rabin(x, 7) && Miller_Rabin(x, 61) && Miller_Rabin(x, 24251));
}

inline long long gcd(long long a, long long b)
{
    if (!a) return b;
    if (!b) return a;
    int t = __builtin_ctzll(a | b);
    a >>= __builtin_ctzll(a);
    do
    {
        b >>= __builtin_ctzll(b);
        if (a > b)
        {
            long long tt = a;
            a = b;
            b = tt;
        }
        b -= a;
    }
    while(b != 0);
    return a << t;
}

long long Pollards_Rho(long long n)
{
    if (n % 2 == 0)
    	return 2;
    if (n % 3 == 0)
    	return 3;
    long long x = 0, y = x, t = 1, q = 1, c = (rand() % (n - 1)) + 1;
    for (register int k = 2; ; k <<= 1, y = x, q = 1)
    {
    	for (register int i = 1; i <= k; ++i)
    	{
    		x = (((__int128)x * (__int128)x) % n + c) % n;
    		q = ((__int128)q * (__int128)abs(x - y)) % n;
    		if ((i & 127) == 0)
            {
                t = gcd(q, n);
                if (t > 1)
                    break;
            }
    	}
    	if (t > 1 || (t = gcd(q, n)) > 1)
    		break;
    }
    return t;
}

void find(long long x, long long &ans)
{
//	printf("find %lld %lld\n", x, ans);
    if (x <= ans || x == 1)
        return;
    if (prime(x))
    {
        ans = ans > x ? ans : x;
        return;
    }
    long long t = x;
    while (t == x)
        t = Pollards_Rho(x);
    while (x % t == 0)
        x /= t;
    find(t, ans);
    find(x, ans);
}

int main()
{//freopen("d.txt", "r", stdin);
//	freopen("w.txt", "w", stdout);
    long long T;
    read(T);
    while (T--)
    {
        long long x;
        read(x);
        long long ans = 0;
        find(x, ans);
        if (x == ans)
            printf("Prime\n");
        else
            printf("%lld\n", ans);	
    }
    return 0;
}