题解：P2063 二平方和定理

littlez_meow · 2025-04-20 19:52:55 · 题解

我们把问题丢到复平面上，就是求 (x+yi)(x-yi)=n 的所有解。

考虑在高斯整环 \mathbb Z[\sqrt{-1}]（就是实部和虚部都是整数的复数集合）分解质因数。我们接下来先讨论分解质因数的方式，再讨论如何得到答案。

根据费马平方和定理，对于一个奇质数 p，若 p\equiv 1\pmod 4，那么它可以表示成一对共轭复数之积；否则，p 在高斯整环下仍是质数，即高斯素数。

也就是说，我们只需要把 n 现在整环上质因数分解，再把 n 的模 4 余 1 的质因数分解成一对共轭复数之积即可。

设奇质数 p=z\bar z 满足 p\equiv 1\pmod 4。我们若能找到整数 x 满足 x^2\equiv -1\pmod p，即 p|(x+i)(x-i)，则 z=\gcd(x+i,p)。

根据二次剩余相关知识，x 是存在的。我们随机一个 p 的二次非剩余 y，期望两次随到，那么 x\equiv y^{\frac{p-1}4}\pmod p。

由于高斯整环是辗转相除环，这里的 \gcd 可以用辗转相除法做，a\bmod b 定义为：将 \dfrac a b 的实部和虚部都四舍五入得到 q，则 a\bmod b=a-bq。由于模长每次减半，复杂度仍是 O(\log n)。

那么我们已经可以把 n 在高斯整环上质因数分解了。设 n=z\bar z，我们应该如何分配这些高斯素数到 z,\bar z 中呢？

先处理模 4 余 3 型质数。这些质数本身就是高斯素数，为了让两个因数共轭，必须均分到 z,\bar z 中。如果不能均分，答案就是 0。

然后是模 4 余 1 型质数。设其中一个质数 p=x\bar x，它的指数为 c，那么我们可以分配 t 个 x 和 c-t 个 \bar x 给 z，剩下的给 \bar z。一共 c 种分配方法。对每个这类质数枚举每一种分配方法，时间复杂度的上界是 O(\tau(n))，其中 \tau 是约数个数函数。

最后是 2。你发现虽然 2=(1+i)(1-i)，但无论怎么分配，得到的解都只是旋转了 90^\circ 的区别，我们可以在最后一步全部转到第一象限，所以怎么分都无所谓。

时间复杂度为 O(T(n^{\frac 1 4}+\omega(n)\log n+\tau(n)\log n))，第一项是整环上分解质因数，第二项是在高斯整环上分解（需要 \gcd），第三项是枚举分配方式（需要快速幂）。

代码

由于我自己写的 SQUFOF 假掉了，在这里粘了一个板子，就不放分解质因数的代码了。

#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
#define int ll
#define fi first
#define se second
using namespace std;
__gnu_pbds::gp_hash_table<ll,int>e;
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;
}
namespace Factorize{
    inline void factorize(ll n,ll cnt=1){

    }
} 
struct Cpx{
    __int128 a,b;//a+bi;
    Cpx(const ll&x=0,const ll&y=0):a(x),b(y){}
    Cpx operator+(const Cpx&x)const{return Cpx(a+x.a,b+x.b);}
    Cpx operator-(const Cpx&x)const{return Cpx(a-x.a,b-x.b);}
    Cpx operator*(const Cpx&x)const{return Cpx(a*x.a-b*x.b,a*x.b+b*x.a);}
    Cpx operator/(const Cpx&x)const{
        double qwq=x.a*x.a+x.b*x.b;
        return Cpx(round((a*x.a+b*x.b)/qwq),round((b*x.a-a*x.b)/qwq));
    }
    Cpx operator%(const Cpx&x)const{
        return *this-x*(*this/x);
    }
    Cpx operator^(ll ex)const{
        Cpx res(1),b=*this;
        while(ex) (ex&1)&&(res=res*b,1),b=b*b,ex>>=1;
        return res;
    } 
    bool operator==(const Cpx&x)const{return a==x.a&&b==x.b;}
    __int128 abs()const{return a*a+b*b;}
};
Cpx gcd(const Cpx&x,const Cpx&y){
    if(x.abs()<y.abs()) return gcd(y,x);
    return y==Cpx()?x:gcd(y,x%y);
}
int cnt,ex[100];
ll pr[100];
Cpx coef[100];
mt19937 gen(114514);
vector<pair<ll,ll> >ans;
void dfs(int step,const Cpx&now){
    if(step==cnt+1){
        ll x=now.a,y=now.b;
        if(x>=0&&y>=0) ans.emplace_back(x,y);
        if(x<=0&&y<=0) ans.emplace_back(-x,-y);
        if(y<=0&&x>=0) ans.emplace_back(-y,x);
        if(y>=0&&x<=0) ans.emplace_back(y,-x);
        return;
    }
    if((pr[step]&3)!=1) dfs(step+1,now*coef[step]);
    else{
        Cpx bar(coef[step].a,-coef[step].b),qwq=bar^ex[step];
        F(i,0,ex[step]) dfs(step+1,now*qwq),qwq=qwq/bar*coef[step];
    }
    return;
}
inline void solve(){
    cnt=0;
    for(auto i:e){
        pr[++cnt]=i.fi,ex[cnt]=i.se;
        if(i.fi==2) coef[cnt]=Cpx(1,1)^i.se;
        else if((i.fi&3)==1){
            ll qwq;
            while(1){
                qwq=gen()%(i.fi-1)+1;
                if(qpow(qwq,i.fi>>1,i.fi)==i.fi-1)break;
            }
            qwq=qpow(qwq,i.fi>>2,i.fi);
            coef[cnt]=gcd(Cpx(i.fi),Cpx(qwq,1));
        }else if(i.se&1) return cout<<"0\n",void();
        else coef[cnt]=Cpx(qpow(i.fi,i.se>>1,1e18));
    }
    vector<pair<ll,ll> >().swap(ans);
    dfs(1,Cpx(1,0));
    sort(ans.begin(),ans.end());
    cout<<ans.size()<<"\n";
    for(auto i:ans) cout<<i.fi<<" "<<i.se<<"\n";
    return;
} 
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    ll T,n;
    for(cin>>T;T;--T){
        e.clear();
        cin>>n;
        if(n<=1e10){
            for(int i=2;i*i<=n;++i) if(n%i==0) while(n%i==0) n/=i,++e[i];
            if(n!=1) ++e[n];
        }else Factorize::factorize(n);
        solve();
    }
    return 0;
}