Miller-Rabin 素性测试
Miller-Rabin 素数测试
一本通上的M-R不透彻啊~
Miller-Rabin是利用随机化算法判断一个数是合数还是素数。
首先，如果一个数N是素数，那么他一定满足费马小定理。
(a^{N-1}equiv1pmod N)
我们可以任取数字a，计算这个式子的值来判断N是否为素数。
但是这么做不靠谱啊，有很多合数会被卡~
我们介绍一个相关的引理。
当p是素数且p大于2时，(1mod p)的平方根只有1和-1。
证明：
假设x是(1mod p)的平方根，则有(x^2equiv1pmod p)
则有((x+1)(x-1)equiv 0pmod p)
也就是说((x+1)(x-1))能够整除素数p。根据欧几里得引理，(x+1)或(x-1)整除p，也就是(xequiv1pmod p)和(xequiv-1pmod p)
我们假设n是一个素数(n>2)，则n-1是一个偶数，且能被表示成(2^s*d)的形式，d是奇数。
那么对于任意的a和(0le rle s-1)，有(a^dequiv1pmod n)和(a^{2^rd}equiv -1pmod n)。
为什么呢？因为费马小定理。对于一个素数n，有(a^{n-1}equiv1pmod n)，由于上面的引理，如果我们不断对(a^{n-1})取平方根，我们总会得到1或-1。如果我们得到了-1，则(a^{2^rd}equiv -1pmod n)成立，判定p是素数。如果我们从未得到-1，通过这个过程我们取遍了2的所有幂次，则(a^dequiv1pmod n)成立。
Miller-Rabin素性测试就在于上述原理的逆否，如果我们找到了一个a，使得对于任意的(0le rle s-1)，满足以下两个式子：(a^d otequiv1pmod n)，(a^{2^rd} otequiv-1pmod 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
", 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
");
else
printf("%lld
", ans);
}
return 0;
}