试除法是最显而易见的一种质因数分解的方法.我们可以从小到大枚举质因子,然后将它一步步从原数中剔除.由于唯一分解定理,我们可知这个分解结果是唯一的.
对于质因数分解我们可以显然地估计出一个上界:$\mathtt{O}(\sqrt{n})$.因为一个合数必定有不大于$\sqrt{n}$的质因子,反证法易证,那么我们只需要当枚举的这个数大于当前整数的$\sqrt{n}$时跳出,剩下的那个数一定是一个质数.
对于试除法我们可以用更好的生成素数的方法进行常数优化.对于大整数这个方法的效率极端低下.
struct prime_exponent{
long prime,exponent;
prime_exponent(long prime=2,long exponent=0):prime(prime),exponent(exponent){}
};
struct trial_division_factor{
long length;
prime_exponent p[40];
inline void clear(){
/* you don't need to clear manually,
// it will be automatically called when a
// new factorization session is started. */
for(long i=0;i<40;++i) p[i]=prime_exponent(2,0);
}
inline long sqr(long p){ return p*p; }
inline long operator()(long inp){
long factor_now=2;
length=0;
while(inp>=sqr(factor_now)){
if(inp%factor_now){
++factor_now;
continue;
}else{
p[length].prime=factor_now;
p[length].exponent=0;
while(!(inp%factor_now)) ++p[length].exponent,inp/=factor_now;
++length;
}
}
if(inp!=1){
p[length].prime=factor_now;
p[length].exponent=1;
++length;
}
return length;
}
};
一个不含因子2的合数$q$可以写成$q=x^2-y^2$,其中$x,y$都是整数.证明:假设$q=q_1q_2$,因为$q$为奇数所以$q_1,q_2$都为奇数,那么$q1+q2$与$q1-q2$都为偶数.这样,令$x=\frac{q_1+q_2}{2},y=\frac{q_1-q_2}{2}$,显然有$x+y=q_1,x-y=q_2$,那么此时$(x+y)(x-y)=x^2-y^2$,且$x,y$都是偶数.
费马法的思想就是这么来的:构造出这里的$x$和$y$.我们使用暴力枚举法.这个算法对于两个大质因数比较接近的数是比较快的,它的最坏复杂度是$O(\sqrt{n})$.其实这相当于另外一种试除法,即从$\sqrt{n}$向0试除.
有最快的一些现代的分解算法(QS,GNFS等)是基于费马法的
简单优化区别与复杂优化的重要一点是,时间增长级别的优化.比如指数级相对于亚指数级.费马法有一些不错的优化方法.
这个非常显然,只需要设置一个试除法处理到的上界(或费马法的下界)
注意到若是$x=n^2$,则$x\equiv n^2\pmod{\mathtt{[some~number]}}$.
而我们需要找的,实际上是一个$a^2-N=b^2$.
那么有一个性质,就是$(a+p)^2\equiv a^2\pmod{p}$.这个性质有什么用呢?
其实这启发了我们可以在搜索时跳过一些不必要的$a$值.假设$a^2-N\not= b^2$,那$(a+p)^2-N\not=b^2$.
那么我们就可以设计出一个算法,它可以跳过一些不必要搜索的值.我们从小到大选取一堆素数,先对最小的那个素数$p_1$,找出一个$a^2-N\equiv b^2\pmod{p_1}$,再以$p_1$一步跳跃着找$a^2-N\equiv b^2\pmod{p+2}$,再以$p_1p_2$一步跳着找...这个过程需要递归,因为二次剩余不只一个,我们需要对所有的都找过去.
//pseudo code
long fermat_with_sieve(long k,long a,long a_end,long a_step,long modulus){
long p=a,q=0;
while(q<modulus){
if(p>a_end) return -1;
long b2=a*a-k;
if(is_quadratic_residue(b2%modulus,modulus)){
if(no_next_prime()) return a;
long pp=fermat_with_sieve(k,a,a_end,a_step*modulus,next_prime(modulus));
if(~pp) return pp;
}
a+=a_step;
++q;
}
return -1;
}
John Pollard于1974年发表了第一个质因数分解算法,即这个Pollard p-1算法.
基本思想:著名的费马小定理说明$a^{k(p-1)}\equiv 1\pmod{p}$,也就是说$a^{k(p-1)}\equiv 0\pmod{p}$,那么若$N=tp$,则$a^{k(p-1)}\bmod N\mid t$.那么,如果一个给定整数$N$拥有一个素因子$p$,其$p-1$可以分解为一些小素数之积,那我们就可以猜测一个以$p-1$为一个因子的数$q$,那么$\gcd{a^q\bmod N,N}$($a$可以任意选定,但是要和$N$互质)就是$N$的一个因子,很可能是非平凡因子(不是1或$N$).
但是在不知道$N$的因子时如何猜测一个数$q$使得$p-q\mid q$呢?大多数$p-1$可以分解为一些小素数的积,那么我们只需要将一些较小的素数都乘起来就可以得到$q$了.此时我们选择一个值$B$,使得$q=\prod_{t=\mathtt{prime}(i)\le B}t^{\lfloor \log_t{B}\rfloor}$就可以完成任务了.当然,这个算法真正的效果是不确定的,甚至不一定能完成分解,而取的值$B$越大就约可能完成分解.
注意到在实现中,我们并不需要将$B$值完整地计算出来;我们只需要使用快速幂的技巧来在模$N$下一边幂一边运算.
这个算法的复杂度是$\mathtt{O}(B\log{B}\log^2{N})$.在实际使用中,有时会取$B=\mathtt{O}(n^{^1/_6})$,这有$\frac{1}{27}$的机率能分解成功.对于$B=\mathtt{O}(n^{^1/_{2k}})$有$k^{-k}$概率分解成功,所以这个算法是指数级的.
在OI/ACM当中,Pollard $\rho$是一种常用的快速分解质因数的算法.这个算法由John Pollard于1975年发明的.这个算法基于Floyd找环算法,其基本思想就是生成一个$[0,N-1]$间的随机序列,在模$N$的一个约数下找到环.一个合数$N$必定有不大于$\sqrt{N}$的因数,而要使一个数字范围为$[0,N-1]$间随机的序列中有大于50%的可能性出现相同的数的序列长度是$\mathtt{O}(\sqrt{N})$的(生日攻击),那么我们只需要期望$\mathtt{O}(\sqrt{p})$的数就可以找到$N$的一个小约数$p$,即时间复杂度最大期望$\mathtt{O}\left(n^{^1/_4}\right)$.更好的估计是$\mathtt{O}(\sqrt{p})$,其中$p$为这个数最小的约数.这个方法对于有着小约数的数更有效.
对于随机数列的生成,我们使用多项式作为生成工具,如$G(n)=n^2+1$.这样,每次数与下一个数连边,可以构成一个类似于$\rho$形状的图.
注意:随着$N$的增大,分解的成功率是增大的,这个可以通过一些很简单的计算验证.
算法过程:
这里找环可以有一些技巧,比如说,用Brent找环算法代替Floyd找环算法.下面是一个应用了Brent找环算法的Pollard $\rho$算法实例程序
//author: zball no rights reserved.
typedef long long ll;//long long当然可以改成高精度数
ll Pollard_rho(ll n,ll c,ll u){
ll i=1,k=2;
ll y=u,x0=u;
while(1){
++i;
x0=(modmul(x0,x0,n)+c)%n;//注意modmul需要使用快速乘或(各种奇怪的方法)实现
ll d=gcd(y-x0,n);//求最大公约数随便你是欧几里得还是Stein还是HGCD都可以
if(d!=1 && d!=n) return d;
if(y==x0) return n;
if(i==k) y=x0,k+=k;
}
}