bzoj4802 欧拉函数(附Millar-Rabin和Pollard-Rho讲解)

传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4802

【题解】

参考:http://www.matrix67.com/blog/archives/234

Millar-Rabin质数检验方法:

根据费马小定理,如果p是素数,a<p,那么有a^(p-1) mod p = 1。

直观想法我们直接取若干个a,如果都有一个不满足,那么p就是合数。

遗憾的是,存在Carmichael数:你无论取多少个a,有一个不满足,算我输。

比如:561 = 11*51就是一个Carmichael数。

bzoj4802 欧拉函数(附Millar-Rabin和Pollard-Rho讲解)

那么就很江了啊。。我们需要改进算法。

首先有:如果p是素数,x是小于p的正整数,且x^2 mod p = 1,那么要么x=1,要么x=p-1

(这个废话,x=p-1模意义下等于x=-1)

然后我们可以展示下341满足2^340 mod 341 = 1,却不是素数(341=31*11)的原因:

2^340 mod 341 = 1

2^170 mod 341 = 1

2^85 mod 341 = 32

(32这个数很py啊怎么不等于340也不等于1啊。。这明显有交易嘛)

那么就能说明这个数不是素数。

如果是素数,一定是从p-1变到1,或是把所有2的次幂去除完,本来就等于1(这样平方完就一直是1了)

所以要么把所有2的次幂去除完,本来就等于1,要么存在某一个次幂=p-1(这样就正常多了)

这就是Millar-Rabin素数验证的二次探测。

应该来说Millar-Rabin算法也是挺好写的

其中mul(a,b,c)表示a*b%c(因为a*b会爆longlong,所以用快速加)

namespace Millar_Rabin {
const int Prime[] = {, , , , , , , , , , , , , };
const int PN = ; inline bool witness(int pr, ll res, int times, ll n) {
ll p = pwr2((ll)pr, res, n);
if(p == ) return ;
for (int i=; i<times; ++i) {
if(p == n-) return false;
if(p == ) return false;
p = mul(p, p, n);
}
return true;
} inline bool main(ll n) {
for (int i=; i<=PN; ++i) {
if(n == Prime[i]) return ;
if(n % Prime[i] == ) return ;
}
ll p = n-;
int times = ;
while(!(p&)) {
++times;
p >>= ;
}
for (int i=; i<=PN; ++i)
if(witness(Prime[i], p, times, n)) return false;
return true;
}
}

然后我们会检验素数了,现在要质因数分解。

好了下一个是Pollard-Rho算法:

如果现在拆分的是n:Pollard-Rho(n)

主要流程:Millar-Rabin判断是否质数,是返回,否就试图找出其中一个因子d,然后递归做Pollard-Rho(d)和Pollard-Rho(n/d)。

我猜你会说废话这谁都会。问题在于:试图找出其中一个因子d

参考:https://wenku.baidu.com/view/3db5c7a6ad51f01dc381f156.html?from=search

参考文章讲的非常详细了。。我就不细讲了qwq

所以这题就是分解因数,按照欧拉函数定义式求解即可。

# include <stdio.h>
# include <string.h>
# include <iostream>
# include <algorithm>
// # include <bits/stdc++.h> using namespace std; typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
const int M = 5e5 + ;
const int mod = 1e9+; # define RG register
# define ST static inline ll mul(ll a, ll b, ll mod) {
ll ret = ;
a %= mod, b %= mod;
while(b) {
if(b&) {
ret = ret + a;
if(ret >= mod) ret -= mod;
}
a <<= ;
if(a >= mod) a -= mod;
b >>= ;
}
return ret;
} inline ll pwr2(ll a, ll b, ll mod) {
ll ret = ;
a %= mod;
while(b) {
if(b&) ret = mul(ret, a, mod);
a = mul(a, a, mod);
b >>= ;
}
return ret;
} inline ll gcd(ll a, ll b) {
return b== ? a : gcd(b, a%b);
} namespace Millar_Rabin {
const int Prime[] = {, , , , , , , , , , , , , };
const int PN = ; inline bool witness(int pr, ll res, int times, ll n) {
ll p = pwr2((ll)pr, res, n);
if(p == ) return ;
for (int i=; i<times; ++i) {
if(p == n-) return false;
if(p == ) return false;
p = mul(p, p, n);
}
return true;
} inline bool main(ll n) {
for (int i=; i<=PN; ++i) {
if(n == Prime[i]) return ;
if(n % Prime[i] == ) return ;
}
ll p = n-;
int times = ;
while(!(p&)) {
++times;
p >>= ;
}
for (int i=; i<=PN; ++i)
if(witness(Prime[i], p, times, n)) return false;
return true;
}
} namespace PollardRho {
const int N = ;
ll q[N]; int qn; inline void PR(ll n) {
if(Millar_Rabin::main(n)) {
q[++qn] = n;
return ;
}
ll a, b, c, del;
while() {
c = rand() % n;
a = b = rand() % n;
b = (mul(b, b, n) + c) % n;
while(a != b) {
del = a-b;
del = gcd(abs(del), n);
if(del > && del < n) {
PR(del); PR(n/del);
return ;
}
a = (mul(a, a, n) + c) % n;
b = (mul(b, b, n) + c) % n;
b = (mul(b, b, n) + c) % n;
}
}
} inline ll getphi(ll n) {
if(n == ) return 1ll;
sort(q+, q+qn+);
ll res = q[] - ;
for (int i=; i<=qn; ++i) {
if(q[i] != q[i-]) res = res * (q[i] - );
else res = res * q[i];
}
return res;
} inline void main(ll n) {
qn = ;
PR(n);
cout << getphi(n) << endl;
}
} int main() {
srand();
ll n; cin >> n;
if(n == ) {
puts("");
return ;
}
PollardRho::main(n);
return ;
}
上一篇:Cognos权限认证CJP方式之用户密码加密


下一篇:MySql 学习之路-聚合函数