SQUFOF算法

来源：互联网发布：黑暗之魂3卡顿优化补丁编辑：程序博客网时间：2024/05/22 15:05

SQUFOF算法：
参考：维基百科 http://en.wikipedia.org/wiki/SQUFOF#Algorithm

思路：寻找一对整数x和y使其满足x^2=y^2(mod n),此时计算gcd(x-y,n)就会得到n的一个素因子，此时素因子在x和y之间

输入：输入一个n和k，要求n既不是完全平方数也不是素数，k是一个乘数
输出：n的一个素因子

初始化： $P_0=\lfloor\sqrt{kN}\rfloor,Q_0=1,Q_1=kN-P_0^2.$

重复执行如下操作：

$b_i=\left\lfloor\frac{\lfloor\sqrt{kN}\rfloor+P_{i-1}}{Q_i}\right\rfloor,P_i=b_iQ_i-P_{i-1},Q_{i+1}=Q_{i-1}+b_i(P_{i-1}-P_i)$

直到 $Q_i$ 是完全平方数.

初始化： $b_0=\left\lfloor\frac{\lfloor\sqrt{kN}\rfloor-P_{i-1}}{\sqrt{Q_i}}\right\rfloor,P_0=b_0\sqrt{Q_i}+P_{i-1},Q_0=\sqrt{Q_i},Q_1=\frac{kN-P_0^2}{Q_0}$

注意，此时的i是上一个循环终止时的i值

重复执行如下操作：

$b_i=\left\lfloor\frac{\lfloor\sqrt{kN}\rfloor+P_{i-1}}{Q_i}\right\rfloor,P_i=b_iQ_i-P_{i-1},Q_{i+1}=Q_{i-1}+b_i(P_{i-1}-P_i)$

直到 $P_{i+1}=P_i.$

T如果 $f=\gcd(N,P_i)$ ！= $1$ || $f=\gcd(N,P_i)$ != $N$ , 此时 $f$ 就是 $N$ 的一个素因子. 否则通过改变 $k$ 的值的大小来做这些操作.

注意，在这儿的 $f=\gcd(N,P_i)$ 中的Pi是最后 $P_{i+1}=P_i.$ 相等时的Pi

一个例子

N = 11111, k = 1

P₀ = 105 Q₀ = 1 Q₁ = 86

P₁ = 67 Q₁ = 86 Q₂ = 77

P₂ = 87 Q₂ = 77 Q₃ = 46

P₃ = 97 Q₃ = 46 Q₄ = 37

P₄ = 88 Q₄ = 37 Q₅ = 91

P₅ = 94 Q₅ = 91 Q₆ = 25

Here Q₆ is a perfect square

P₀ = 104 Q₀ = 5 Q₁ = 59

P₁ = 73 Q₁ = 59 Q₂ = 98

P₂ = 25 Q₂ = 98 Q₃ = 107

P₃ = 82 Q₃ = 107 Q₄ = 41

P₄ = 82

Here P₃ = P₄

gcd(11111, 82) = 41, which is a factor of 11111.

const int maxn=100;long long gcd(long long a,long long b){    return b==0?a:gcd(b,a%b);}long long b[maxn],p[maxn],q[maxn];long long SQUFOF(long long n,long long k){    long long a=(long long )(sqrt(k*n*1.0));    long long t;    int i;    k=1;    p[0]=a;    q[0]=1;    q[1]=k*n-p[0]*p[0];    for(i=1;i<maxn;i++)    {        b[i]=(long long)((a+p[i-1])/q[i]);        p[i]=b[i]*q[i]-p[i-1];        q[i+1]=q[i-1]+b[i]*(p[i-1]-p[i]);        t=(long long)sqrt(q[i]*1.0);        if(t*t==q[i])break;    }    b[0]=(long long)((a-p[i-1])/t);    p[0]=b[0]*t+p[i-1];    q[0]=t;    q[1]=(k*n-p[0]*p[0])/q[0];    for(i=1;i<maxn;i++)    {        b[i]=(long long)((a+p[i-1])/q[i]);        p[i]=b[i]*q[i]-p[i-1];        q[i+1]=q[i-1]+b[i]*(p[i-1]-p[i]);        if(p[i+1]==p[i])break;    }    long long f=gcd(n,p[i]);    if(f!=1||f!=n)return f;    else k++;}