专栏文章

题解:P2063 二平方和定理

P2063题解参与者 3已保存评论 2

文章操作

快速查看文章及其快照的属性,并进行相关操作。

当前评论
2 条
当前快照
1 份
快照标识符
@mipjnj8x
此快照首次捕获于
2025/12/03 13:05
3 个月前
此快照最后确认于
2025/12/03 13:05
3 个月前
查看原文
我们把问题丢到复平面上,就是求 (x+yi)(xyi)=n(x+yi)(x-yi)=n 的所有解。
考虑在高斯整环 Z[1]\mathbb Z[\sqrt{-1}](就是实部和虚部都是整数的复数集合)分解质因数。我们接下来先讨论分解质因数的方式,再讨论如何得到答案。
根据费马平方和定理,对于一个奇质数 pp,若 p1(mod4)p\equiv 1\pmod 4,那么它可以表示成一对共轭复数之积;否则,pp 在高斯整环下仍是质数,即高斯素数。
也就是说,我们只需要把 nn 现在整环上质因数分解,再把 nn 的模 4411 的质因数分解成一对共轭复数之积即可。
设奇质数 p=zzˉp=z\bar z 满足 p1(mod4)p\equiv 1\pmod 4。我们若能找到整数 xx 满足 x21(modp)x^2\equiv -1\pmod p,即 p(x+i)(xi)p|(x+i)(x-i),则 z=gcd(x+i,p)z=\gcd(x+i,p)
根据二次剩余相关知识,xx 是存在的。我们随机一个 pp 的二次非剩余 yy,期望两次随到,那么 xyp14(modp)x\equiv y^{\frac{p-1}4}\pmod p
由于高斯整环是辗转相除环,这里的 gcd\gcd 可以用辗转相除法做,amodba\bmod b 定义为:将 ab\dfrac a b 的实部和虚部都四舍五入得到 qq,则 amodb=abqa\bmod b=a-bq。由于模长每次减半,复杂度仍是 O(logn)O(\log n)
那么我们已经可以把 nn 在高斯整环上质因数分解了。设 n=zzˉn=z\bar z,我们应该如何分配这些高斯素数到 z,zˉz,\bar z 中呢?
先处理模 4433 型质数。这些质数本身就是高斯素数,为了让两个因数共轭,必须均分到 z,zˉz,\bar z 中。如果不能均分,答案就是 00
然后是模 4411 型质数。设其中一个质数 p=xxˉp=x\bar x,它的指数为 cc,那么我们可以分配 ttxxctc-txˉ\bar xzz,剩下的给 zˉ\bar z。一共 cc 种分配方法。对每个这类质数枚举每一种分配方法,时间复杂度的上界是 O(τ(n))O(\tau(n)),其中 τ\tau 是约数个数函数。
最后是 22。你发现虽然 2=(1+i)(1i)2=(1+i)(1-i),但无论怎么分配,得到的解都只是旋转了 9090^\circ 的区别,我们可以在最后一步全部转到第一象限,所以怎么分都无所谓。
时间复杂度为 O(T(n14+ω(n)logn+τ(n)logn))O(T(n^{\frac 1 4}+\omega(n)\log n+\tau(n)\log n)),第一项是整环上分解质因数,第二项是在高斯整环上分解(需要 gcd\gcd),第三项是枚举分配方式(需要快速幂)。

代码

由于我自己写的 SQUFOF 假掉了,在这里粘了一个板子,就不放分解质因数的代码了。
CPP
#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;
}

评论

2 条评论,欢迎与作者交流。

正在加载评论...