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

时间:2023-03-09 04:05:37
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;
}