自动搬运

查看原文

来自洛谷，原作者为

	littlez_meow 这个小z很懒，连这个家伙很懒都没有留下

搬运于2025-08-24 21:33:34，当前版本为作者最后更新于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;
}