挑战性题目DSCT501:大整数因子分解
问题描述
大整数质因子分解。
题解
笑死,数据结构课背包都不讲,却给学生出泼辣肉,我的心里只有感恩。
因为本题的数据范围巨大,传统的 O ( n ) O(\sqrt n) O(n)的质因数分解已经不适用,需要另外的算法。
首先是质数探测算法——Miller_rabin算法,其基于费马小定理,若p是质数,
a
p
−
1
≡
1
(
m
o
d
p
)
a^{p-1}\equiv1\left(\mod p\right)
ap−1≡1(modp)。同时,若
p
p
p为质数,
x
2
≡
1
(
m
o
d
p
)
x^2\equiv1\left(\mod p\right)
x2≡1(modp)可推出
x
=
1
o
r
p
−
1
x=1\ or\ p-1
x=1 or p−1。所以,我们要判断一个质数是否是质数时,先套费马小定理判断,如果符合费马小定理,再判断
(
n
−
1
)
/
2
(n-1)/2
(n−1)/2是否满足上述第二个性质,如果对于一个预先设定的质数集合p[]
(本人使用的集合为
2
,
325
,
9375
,
28178
,
450775
,
9780504
,
1795265022
2, 325, 9375, 28178, 450775, 9780504, 1795265022
2,325,9375,28178,450775,9780504,1795265022),数字
n
n
n均通过上述检验,则认为其是质数。
之后,再利用大整数因数探测算法Pollar_Rho算法,若正整数 n n n至少有一个因子不超过 n \sqrt n n,设它任意一个因子为 p p p,若随机生成 p p p以内的数,期望 O ( p ) O\left(\sqrt p\right) O(p)次就可以随机到两个一样的数。
随机生成 [ 1 , n ] [1,n] [1,n]的正整数,它们 m o d p mod\ p mod p下可近似看做在 [ 1 , p ] [1,p] [1,p]内随机,随机生成两个数 a , b ( b > a ) a,b(b>a) a,b(b>a),若 b = a ( m o d p ) → p ∣ ( b − a ) → g c d ( n , b − a ) ≥ p b=a(mod\ p)\rightarrow p|(b-a)\rightarrow gcd(n,\ b-a)\geq p b=a(mod p)→p∣(b−a)→gcd(n, b−a)≥p
那么 g c d ( n , b − a ) gcd\left(n,b-a\right) gcd(n,b−a)必为 n n n的一个大于 1 1 1的因数。令随机数这样生成: a 0 = C , f ( x ) = x 2 + C , a i = f ( a i − 1 ) a_0=C,f\left(x\right)=x^2+C,a_i=f\left(a_{i-1}\right) a0=C,f(x)=x2+C,ai=f(ai−1)即跳一步就套用一次 f ( x ) f\left(x\right) f(x)。这样做的好处是:若 a x = a y ( m o d p ) → a x + 1 − a y + 1 = f ( a x ) − f ( a y ) = ( a x + a y ) ( a x − a y ) = 0 ( m o d p ) a_x=a_y\left(\mod p\right)\rightarrow a_{x+1}-a_{y+1}=f\left(a_x\right)-f\left(a_y\right)=\left(a_x+a_y\right)\left(a_x-a_y\right)=0(\mod p) ax=ay(modp)→ax+1−ay+1=f(ax)−f(ay)=(ax+ay)(ax−ay)=0(modp) 即任一距离为 d d d的随机数 m o d p \mod p modp相同,所有距离为 d d d的随机数 m o d p \mod p modp都相同,这样只需要检查其中一个即可。
免除记忆化:维护两个数,一个数每次调用两次 f f f,一个每次调用一次,如果走完了环,最后两个数一定会相遇。同时也枚举了所有距离。
优化:瓶颈在于 g c d gcd gcd,可以减少取 g c d ( n , ∣ a − b ∣ ) gcd\left(n,\left|a-b\right|\right) gcd(n,∣a−b∣)的次数,将所有的 ∣ a − b ∣ \left|a-b\right| ∣a−b∣乘起来,累计 128 128 128次后再与 n n n取 g c d gcd gcd(前提是 a , b a,b a,b还没有在环上相遇),复杂度约为 O ( n 1 / 4 ) O(n^{1/4}) O(n1/4)。
代码
感谢我的超人——Rockdu提供的泼辣肉模板,说实话我也没怎么搞懂,以后有缘再学。
#pragma GCC optimize(2)
#include<bits/stdc++.h>
#define LL long long
using namespace std;
LL n;
namespace NT {
LL gcd(LL a, LL b) {
return b ? gcd(b, a % b) : a;
}
LL mul(LL a, LL b, LL m) {
LL s = a * b - (LL)((long double)a * b / m + 0.5) * m;
return s < 0 ? s + m : s;
}
LL fpow(LL x, LL a, LL m) {
LL ans = 1;
while(a) {
if(a & 1) ans = mul(ans, x, m);
x = mul(x, x, m), a >>= 1;
}
return ans;
}
}
namespace Miller_Rabin {
using namespace NT;
LL p[15] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
int detect(LL n, LL a) {
if(n == a) return 1;
if(a % n == 0) return 1;
LL now = n - 1, lst = 1;
if(fpow(a, now, n) != 1)
return 0;
while(!(now & 1)) {
now /= 2;
LL p = fpow(a, now, n);
if(lst == 1 && p != 1 && p != n - 1)
return 0;
lst = p;
}
return 1;
}
LL MR(LL n) {
if(n < 2) return 0;
for(int i = 0; i < 7; ++i) {
if(!detect(n, p[i]))
return 0;
}
return 1;
}
}
namespace Pollard_Rho {
using namespace NT;
using namespace Miller_Rabin;
LL f(LL x, LL C, LL p) {
return (mul(x, x, p) + C) % p;
}
LL Rand() {return (rand() << 15) + rand();}
LL Randll() {return (Rand() << 31) + Rand();}
LL PR(LL n) {
if(n == 4) return 2;
if(MR(n)) return n;
while(1) {
LL C = Randll() % (n - 1) + 1;
LL s = 0, t = 0;
LL acc = 1;
do {
for(int i = 1; i <= 128; ++i) {
t = f(f(t, C, n), C, n);
s = f(s, C, n);
LL tmp = mul(acc, abs(t - s), n);
if(s == t || tmp == 0)
break;
acc = tmp;
}
LL d = gcd(n, acc);
if(d > 1) return d;
} while(s != t);
}
}
typedef pair<LL, int > pii;
vector<pii > DCOMP(LL n) {
vector<pii > ret;
while(n != 1) {
LL p = PR(n);
while(!MR(p))
p = PR(p);
int c = 0;
while(n % p == 0) n /= p, ++c;
ret.push_back(make_pair(p, c));
}
return ret;
}
}
void out(long long x){if(x>9)out(x/10);putchar(x%10+'0');}
int main(int argc,char* argv[]) {
srand(19260817);
LL n=atoll(argv[1]);
if(n==1)return 0;
auto ans = Pollard_Rho::DCOMP(n);
out(ans[0].first);
ans[0].second--;
for(auto i : ans)
while(i.second--) putchar('*'),out(i.first);
return 0;
}