前言
玄学的概率测素数, 这是一种随机的算法, 出错率大约在$($1/4$)^n$, 所以当n > 50以后正确率就接近于1了(反正比计算机本身的正确率要高)
前置芝士
费马小定理:
当a于p互质时,$(a^p) \equiv a (\bmod p)$
二次探测定理
内容:如果 φ(p)=p−1,p>1,p>X,且$X^2 ≡ 1(\bmod p)$,那么 $X = 1 || p - 1$
证明:
∵$X^2≡1 \pmod p$
∴ $p|X^2 - 1$
∴$p|(X - 1)(X+ 1)$
∵p是大于X的质数
∴ $p = X + 1 || p \equiv X-1\pmod p$ 即$X = 1 || p - 1$
证毕
原理
好, 那么接下来让我们走进它的原理。
首先,如果当满足$(a^p-1) \equiv a \pmod p$时p都为质数那这个问题是不是就迎刃而解了呢?(就可以和愚蠢的二次探测定理说再见了)
但愿望是美好的,现实是骨感的,当你自信地将这个代码交上去以后,你会发现:你Wa了。
我们来看一下这个数据:$2^{340} \equiv 1 \pmod {341}$,然而 341=31∗11
……
实现
- 1.将$p - 1$提出所有2的因数, 我们假设有t个。然后将剩余的部分定义为res(用于二次探测定理)
- 2.枚举一个数a,并定义一个数$x = a^{res} \pmod p $
3.如果 $\forall y = x^2 \pmod p, y != p - 1$ 那么p就不是一个质数
4.当我们的底数已经足够多了时,就可以跳出了
主函数及其他流程
牛B的大佬就可以去刷题了(
别来找茬了)我们先判断当前数是否为素数
如果不是素数的话,就找它的因子
递归对该因子和约出来的另一个因子进行分解(直到为质数)
因子从何而来
如果不嫌弃的话,我们可以一个一个试咳,我们来说正解。
我们假设要找的因子为a
我们运用随机的艺术find一个$x,y$,并不停随机$x$, 具体的法子一般是$x = x^2 + c$(c就是随机艺术的产物)
使$a = gcd(y - x, n) && a \in(1, n)$则我们就找到了一个因子(
至于为什么……)那如果$x = y$出现了呢,这就说明出现了循环,所以我们就要判环,运用倍增的思想,让$y$记住$x$,当$x = {x_0}^2$时,y = x,所以当x跑到y时,已经跑完一个圈。
定义一个$i = 1,j = 2$,当执行$gcd$时$i ++$,如果i == j, 则$y = x, j <<= 1$
代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
using namespace std;
int T;
LL n, ans;
LL a[30] = {2,3,5,7,11,13,17,41,61,24251};
LL Gcd(LL x, LL y)
{
return y ? Gcd(y, x % y) : x;
}
LL mul(LL x, LL k, LL mod)
{
LL res = 0;
while(k)
{
if(k & 1) res += x, res %= mod;
x += x;
x %= mod;
k >>= 1;
}
return res % mod;
}
LL qpow(LL x, LL k, LL mod)
{
LL res = 1;
while(k)
{
if(k & 1) res = mul(res, x, mod);
x = mul(x, x, mod);
k >>= 1;
}
return res % mod;
}
bool Robin(LL p)
{
if(p < 2) return false;
if(p == 2) return true;
if(p % 2 == 0) return false;
LL d = p - 1;
int s = 0;
while(!(d & 1)) d >>= 1, s ++;
for(int i = 0;i <= 9 && a[i] < p;i ++)
{
LL x = qpow(a[i], d, p), y = 0;
for(int j = 1;j <= s;j ++)
{
y = mul(x, x, p);
if(y == 1 && x != 1 && x != (p - 1)) return false;
x = y;
}
if(y != 1) return false;
}
return true;
}
inline LL pollard(LL p)
{
while(1)
{
LL x = rand() % (p - 2) + 1;
LL y = rand() % (p - 2) + 1;
LL c = rand() % (p - 2) + 1;
LL i = 0, j = 1;
while(++ i)
{
x = (mul(x, x, p) + c) % p;
if(x == y) break;
LL d = Gcd(abs(y - x), p);
if(d > 1) return d;
if(i == j)
{
y = x;
j <<= 1;
}
}
}
}
void Find(LL p)
{
if(p == 1) return ;
if(Robin(p))
{
ans = max(ans, p);
return ;
}
LL k = p;
while(p <= k) k = pollard(k);
Find(k), Find(p / k);
}
void print(LL x)
{
if(x < 0) x = -x, putchar('-');
if(x > 9) print(x / 10);
putchar(x % 10 + '0');
}
int main()
{
scanf("%d", &T);
srand(time(0));
while(T --)
{
ans = 0;
scanf("%lld", &n);
Find(n);
if(ans == n) puts("Prime");
else print(ans),puts("");
}
return 0;
}
我们会发现:咦怎么T了
这时就要卡常(打表)了,经过被巨佬一顿嘲讽,我们终于问到了卡常的办法
快速乘
1 | typedef long double ld; |
连简单的$\log(n)$的龟速乘都不放过……
然后我们就可以得到优秀的93分
最后的优化
93分了,冲鸭
向哪冲呢?
还记得我们的$Gcd$吗,当然,$Gcd$函数是不能再优化了,但我们调用$Gcd$的次数是可以再优化的。
- 我们都知道, 龟速乘的模数都是质数, 但我们的模数可能并不是质数
- 所以我们可以根据取模的性质:如果模数和被模的数都含有一个公约数$k$,那么这次模运算的结果必然也会是这个公约数$k$的倍数。所以如果我们将若干个$(y - x)$ 相乘,因为模数是$n$ ,所以如果若干个$(y - x)$中有一个与$n$有公约数,最后的结果定然也会含有这个公约数。
- 所以我们可以多算几次$(y - x)$, 再来求$Gcd$(63次擦不多吧)
- 记得一边欧拉$(y - x)$,一边倍增,判环哟
AC代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
using namespace std;
int T;
LL n, ans;
LL a[30] = {2,3,5,7,11,13,17,41,61,24251};
LL Gcd(LL x, LL y)
{
return y ? Gcd(y, x % y) : x;
}
typedef long double ld;
typedef long long ll;
typedef unsigned long long ull;
inline ull mul(ull a,ull b, ull mod){
ll c=(ll)(a)*b-(ll)((ull)((ld)(a)*b/mod)*mod);
return c<0? c+mod:(c<mod? c:c-mod);
}
LL qpow(LL x, LL k, LL mod)
{
LL res = 1;
while(k)
{
if(k & 1) res = mul(res, x, mod);
x = mul(x, x, mod);
k >>= 1;
}
return res % mod;
}
bool Robin(LL p)
{
if(p < 2) return false;
if(p == 2) return true;
if(p % 2 == 0) return false;
LL d = p - 1;
int s = 0;
while(!(d & 1)) d >>= 1, s ++;
for(int i = 0;i <= 9 && a[i] < p;i ++)
{
LL x = qpow(a[i], d, p), y = 0;
for(int j = 1;j <= s;j ++)
{
y = mul(x, x, p);
if(y == 1 && x != 1 && x != (p - 1)) return false;
x = y;
}
if(y != 1) return false;
}
return true;
}
inline LL pollard(LL p)
{
while(1)
{
LL x = rand() % (p - 2) + 1;
LL y = rand() % (p - 2) + 1;
LL c = rand() % (p - 2) + 1;
LL i = 0, j = 1, b = 1;
while(++ i)
{
x = (mul(x * 1ull, x * 1ull, p * 1ull) + c) % p;
b = mul(b * 1ull, abs(y - x) * 1ull, p * 1ull);
if(!b || x == y) break;
if(i == j || !(i % 63))
{
LL d = Gcd(b, p);
if(d > 1) return d;
if(i == j) y = x, j <<= 1;
}
}
}
}
void Find(LL p)
{
if(p <= ans || p == 1) return ;//小优化
if(Robin(p))
{
ans = max(ans, p);
return ;
}
LL k = pollard(p);
while(p % k == 0) p /= k;
Find(k), Find(p);
}
void print(LL x)
{
if(x < 0) x = -x, putchar('-');
if(x > 9) print(x / 10);
putchar(x % 10 + '0');
}
int main()
{
scanf("%d", &T);
srand(time(0));
while(T --)
{
ans = 0;
scanf("%lld", &n);
Find(n);
if(ans == n) puts("Prime");
else print(ans),puts("");
}
return 0;
}