SQUFOF：比 Pollard-Rho 更好的质因数分解算法

· · 算法·理论

Part 0：前言与算法简介

做数论题时，你遇到了许多需要 Pollard-Rho 分解质因数的题。

你想去学，不过看着那么多前置知识和优化，还是个随机化算法，甚至要卡常。你望而却步。

但是质因数分解算法还是需要的。你在网上检索，发现了一个名为 SQUFOF 的算法。

SQUFOF，全称 Shanks' SQUare FOrm Factorization，是二十纪七八十年代由 Daniel Shanks 发明的算法。有趣的是，他从未正式发表过该算法。

SQUFOF 相比于 Pollard-Rho 有如下优点：

  1. 在都没有卡常的情况下，SQUFOF 的常数远小于 Pollard-Rho。

  2. 可以证明 SQUFOF 运算次数上界为 t\times3\sqrt{2\sqrt n}=O(n^{0.25})，而非期望。其中 t 取决于你的实现，在 10^{18} 及以下时为 3。

  3. SQUFOF 的实现更简单，不需要倍增判环等等东西，只需递推。

可以说，SQUFOF 是分解 [10^{10},10^{18}] 中的数最快的算法。[1]

下面记需要分解的数为 n，SQUFOF 算法会返回一个 n 的非平凡因数（即不为 1,n）。

Part 1：算法原理

SQUFOF 的使用前提是 n 不为质数、偶数、完全平方数。这都是容易判掉的。

假设我现在找到了一对不相等的数 x,y，满足 x^2\equiv y^2\pmod n，那么，根据平方差公式，有 n|(x+y)(x-y)，即 \gcd(x\pm y,n)|n，且大概率 \gcd 不为 1。

为了找到这样的 x,y，我们将使用连分数。

设 \sqrt n=b_0+\dfrac{1}{b_1+\dfrac 1{\ddots+\dfrac {\vdots}{b_{i-1}+\dfrac{Q_i}{\sqrt n+P_i}}}}

容易得出递推式 (\star)：

b_{0}=\lfloor\sqrt{n}\rfloor, b_{i}=\left\lfloor\frac{\lfloor\sqrt{n}\rfloor+P_{i-1}}{Q_{i}}\right\rfloor P_{1}=b_{0}, P_{i}=b_{i} Q_{i}-P_{i-1} Q_{0}=1, Q_{1}=n-b_{0}^{2}, Q_{i+1}=Q_{i-1}-b_{i}\left(P_{i}-P_{i-1}\right)

将这些递推式带回原连分数容易验证，在此不过多赘述。

设 \dfrac{A_i}{B_i}=[b_0;b_1,\cdots,b_i]（记号参见 oi-wiki），我们有如下重要性质：

A_i^2-nB_i^2=(-1)^iQ_i

简单证明一下。

首先根据连分数写出 A_i,B_i 的递推式如下：

A_{0}=1, A_{1}=b_{0}, A_{i+1}=A_{i} b_{i}+A_{i-1} B_{0}=0, B_{1}=1, B_{i+1}=B_{i} b_{i}+B_{i-1}

既然是关于递推式的命题，考虑数学归纳法。当 i=0,1 时显然满足，假设 i\in[0,t] 时满足，下证当 i=t+1 时也满足。

写出等式左边为 A_{t+1}^2-nB_{t+1}^2。

代入递推式展开得 A_t^2b_t^2+A_{t-1}^2+2A_tA_{t-1}b_t-n(B_t^2b_t^2+B_{t-1}^2+2B_tB_{t-1}b_t)。

整理一下有 b_t^2(A_t^2-nB_t^2)+(A_{t-1}^2-nB_{t-1}^2)+2b_t(A_tA_{t-1}-nB_tB_{t-1})。

代入假设得 b_t^2(-1)^tQ_t+(-1)^{t-1}Q_{t-1}+2b_t(A_tA_{t-1}-nB_tB_{t-1})。

对 (\star) 变形，得到 Q_{i+1}=Q_{i-1}-b_i^2Q_i+2b_{i}P_{i-1}，带入上式有 (-1)^{t+1}Q_{t+1}+(-1)^t2b_tP_{t-1}+2b_t(A_tA_{t-1}-nB_tB_{t-1})。

对比等式右边，我们仅需证明 (-1)^{t-1}P_{t-1}=A_tA_{t-1}-nB_tB_{t-1}。

这可以再次使用数学归纳法，然后螺旋归纳得到结果。

综上，A_i^2-nB_i^2=(-1)^iQ_i 成立。

因此，A_i^2\equiv (-1)^i Q_i\pmod n。只要 (-1)^i Q_i 是完全平方数，我们就可以取 A_i,\sqrt{(-1)^iQ_i} 作为上述 x,y。

Part 2：算法流程

Part 3：代码

#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/hash_policy.hpp>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
using namespace std;
const int MR_Prime[]={2,325,9375,28178,450775,9780504,1795265022};
inline ll qpow(__int128 base,ll expo,const ll MOD){
    ll res(1);
    while(expo){
        (expo&1)&&(res=res*base%MOD);
        base=base*base%MOD,expo>>=1;
    }
    return res;
}
inline bool MR(ll n){
    if(n<=2) return n==2;
    if(!(n&1)) return 0;
    ll expo=n-1;
    int r=__builtin_ctzll(expo);
    expo>>=r;
    for(int i:MR_Prime){
        ll v=qpow(i,expo,n);
        if(v<=1||v==n-1) continue;
        F(j,1,r){
            v=__int128(v)*v%n;
            if(v==n-1&&j!=r){
                v=1;
                break;
            }
            if(v==1) return 0;
        }
        if(v!=1) return 0;
    }
    return 1;
}
const int SQUFOF_K[]={1,3,5};
inline ll SQUFOF(ll n){
    if(!(n&1)) return 2;
    ll qwq=sqrt(n);
    if(qwq*qwq==n) return qwq;
    for(int k:SQUFOF_K){
        if(k!=1&&!(n%k)) return k;
        ll kn=k*n,rn=sqrt(kn);
        auto trans=[&](ll&b,ll&P,ll&Q){
            swap(b,Q);
            ll pre=P,q=(rn+P)/b;
            P=q*b-P;
            Q+=q*(pre-P);
            return;
        };
        ll b=0,P=rn,Q=1;
        int lim=sqrt(rn<<1);
        lim<<=2;
        trans(b,P,Q);
        Q=(kn-P*P)/b;
        for(int i(2);i<lim;i+=2){
            trans(b,P,Q);
            qwq=sqrt(Q);
            if(qwq*qwq==Q){
                ll bb,PP,QQ;
                PP=-P,QQ=qwq;
                trans(bb,PP,QQ);
                QQ=(kn-PP*PP)/bb;
                do{
                    qwq=PP;
                    trans(bb,PP,QQ);
                }while(PP!=qwq);
                qwq=__gcd(n,qwq);
                if(qwq!=1) return qwq;
            }
            trans(b,P,Q);
        }
    }
    return 1;
}
ll ans;
inline void factorize(ll n){
    if(n==1) return;
    if(MR(n)) return ans=max(ans,n),void();
    ll qaq=SQUFOF(n);
    while(qaq==n) qaq=SQUFOF(n);
    while(!(n%qaq)) n/=qaq;
    factorize(n),factorize(qaq);
    return;
}
ll n;
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    int T;
    for(cin>>T;T;--T){
        cin>>n;
        ans=0;
        if(MR(n)) cout<<"Prime\n";
        else factorize(n),cout<<ans<<"\n";
    }
    return 0;
}

Part \infty：参考文献

[1] Jason E. Gower and Samuel S. Wagstaff, JR., Square form factorization, Mathematics of computation 77 (2008), no.261, 551-588.