圆上的整点 题解

题意：给定非负整数\(n\)，求方程\(x^2+y^2=n\)的所有整数解（\((x,y)\)与\((y,x)\)算一组解）
数据范围：\(n\le4\times 10^{18}\)
题解：
根据费马平方和定理，如果 \(n\) 分解质因数后形如 \(4p+3\) 的质数的幂次为奇数那么无解
否则将 \(n\) 分解为高斯素数并暴力 DFS 枚举分配方案求解
因为\(n\)非常大所以要用 Pollard-Rho 分解质因数
对于 \(2\) 显然可以分解为 \(\left(1+i\right)\left(1-i\right)\)
对于 \(4p+1\) 类型的数可以参考这篇博客
具体方法为：

• 用二次剩余求出 \(x_0^2\equiv-1\left(\bmod p\right)\) 的解
  （因为 \(p\) 为 \(4p+1\) 类型的数，所以 \((-1)^{\frac{p-1}{2}}=1\)，所以 \(-1\) 有二次剩余）
• 对 \(p\) 和 \(x_0\) 做辗转相除，当余数小于 \(\sqrt{p}\) 的时候停止
• 此时的 \(r\) 和 \(\sqrt{p-r^2}\) 即为方程 \(x^2+y^2=p\) 的非负整数解
• 证明我不会

对于 \(4p+3\) 类型的数不需要继续分解（因为无法分解）
代码如下：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
#define M 1000005
using namespace std;
long long P[M],num[M],R;
long long read(){
	char c=getchar();long long ans=0;
	while (c<'0'||c>'9') c=getchar();
	while (c>='0'&&c<='9') ans=ans*10+c-'0',c=getchar();
	return ans;
}
void Write(long long x){
	if (x<10) putchar(x^48);
	else Write(x/10),putchar((x%10)^48);
	return;
}
long long add(long long u,long long v,long long mod){u+=v-mod;return u+((u>>63)&mod);}
long long gcd(long long x,long long y){return y?gcd(y,x%y):x;}
long long Abs(long long x){return x>0?x:-x;}
long long mult(long long x,long long y,long long mod){
	if (x<=1000000000&&y<=1000000000) return x*y%mod;
	long long ans=x*y-(long long)((long double)(x)*y/mod)*mod;
	if (ans<0) return ans+mod;
	if (ans>=mod) return ans-mod;
	return ans;
}
long long power(long long x,long long y,long long mod){
	long long ans=1,now=x%mod;
	for (register long long i=y;i;i>>=1,now=mult(now,now,mod))
		if (i&1) ans=mult(ans,now,mod);
	return ans;
}
namespace Pollard_Rho{
	bool check(long long x,long long p){
		long long d=0,tmp=x-1;
		while (!(tmp&1)) tmp>>=1,++d;
		long long cur=power(p,tmp,x);
		if (cur==1||cur==x-1) return 1;
		while (--d) if ((cur=mult(cur,cur,x))==x-1) return 1;
		return 0;
	}
	bool is_prime(long long x){
		if (x==2ll||x==3ll||x==7ll||x==61ll||x==24251ll) return 1;
		if (x<2ll||x==46856248255981ll||!(x&1)) return 0;
		return check(x,2ll)&&check(x,3ll)&&check(x,7ll)&&check(x,61ll)&&check(x,24251ll);
	}
	long long F(long long x,long long c,long long mod){return add(mult(x,x,mod),c,mod);}
	long long get_divide(long long x){
		long long s=0,t=0,c=rand()%(x-1)+1,val=1;
		for (register int goal=1;;goal<<=1,s=t,val=1){
			for (register int step=1;step<=goal;++step){
				t=F(t,c,x),val=mult(val,Abs(t-s),x);
				if (step%127==0){
					long long d=gcd(val,x);
					if (d>1) return d;
				}
			}
			long long d=gcd(val,x);
			if (d>1) return d;
		}
		return 0;
	}
	void divide(long long x){
		if (x<2) return;
		if (is_prime(x)){P[++P[0]]=x;return;}
		long long now=get_divide(x);
		while (x%now==0) x/=now;
		divide(x),divide(now);
		return;
	}
}
namespace Sqrt{
	long long mod,I;
	bool check(long long x){return power(x,(mod-1)>>1,mod)==1;}
	struct Sqrt{
		long long x,y;
		Sqrt operator* (const Sqrt b) const{
			long long X=add(mult(x,b.x,mod),mult(mult(y,b.y,mod),I,mod),mod);
			long long Y=add(mult(x,b.y,mod),mult(y,b.x,mod),mod);
			return (Sqrt){X,Y};
		}
		Sqrt operator*= (const Sqrt &b){return *this=*this*b;}
		Sqrt operator^ (const long long b) const{
			Sqrt ans=(Sqrt){1,0},now=*this;
			for (register long long i=b;i;i>>=1,now*=now)
				if (i&1) ans*=now;
			return ans;
		}
	};
	long long get_sqrt(long long n,long long Mod){;
		mod=Mod;long long now=rand()%mod;
		while (check((mult(now,now,mod)+mod-n)%mod)) now=rand()%mod;
		I=(mult(now,now,mod)-n+mod)%mod;
		return ((Sqrt){now,1}^((mod+1)>>1)).x;
	}
}
namespace Get_ans{
	pair <int,int> ans[M<<2];
	int ans_num;
	struct Complex{
		long long X,Y;
		Complex operator* (const Complex b) const{
			return (Complex){X*b.X-Y*b.Y,X*b.Y+Y*b.X};
		}
		Complex operator*= (const Complex &b){return *this=*this*b;} 
	}B[M];
	Complex work(Complex x){return (Complex){x.X,-x.Y};}
	Complex divide_num(long long x){
		long long tmp=x,X0=Sqrt::get_sqrt(x-1,x),r=x%X0,sqrtx=(long long)(sqrt(x));
		while (r>sqrtx+3||r*r>tmp) x=X0,X0=r,r=x%X0;
		return (Complex){r,(long long)(sqrt(tmp-r*r+0.000001))};
	}
	void DFS(long long x,Complex now){
		if (x==P[0]+1){ans[++ans_num]=make_pair(min(Abs(now.X),Abs(now.Y)),max(Abs(now.X),Abs(now.Y)));return;}
		if (P[x]==2){
			for (register int i=1;i<=num[x];i++) now*=B[x];
			DFS(x+1,now);return;
		}
		if (P[x]%4==3){
			for (register int i=1;i<=(num[x]>>1);i++) now*=B[x];
			DFS(x+1,now);return;
		}
		for (register int i=0;i<=num[x];i++){
			Complex nxt=now;
			for (register int j=1;j<=i;j++) nxt*=B[x];
			for (register int j=1;j<=num[x]-i;j++) nxt*=work(B[x]);
			DFS(x+1,nxt);
		}
		return;
	}
	void solve(){
		for (register long long i=1;i<=P[0];i++){
			long long tmp=R;
			while (tmp%P[i]==0) tmp/=P[i],++num[i];
			if (P[i]==2) B[i]=(Complex){1,1};
			else if (P[i]%4==1) B[i]=divide_num(P[i]);
			else {
				if (P[i]%4==3&&(num[i]&1)){printf("No Solution");return;}
				else B[i]=(Complex){P[i],0};
			}
		}
		DFS(1,(Complex){1,0});
		sort(ans+1,ans+1+ans_num);ans_num=unique(ans+1,ans+1+ans_num)-(ans+1);Write(ans_num),putchar('\n');
		for (register int i=1;i<=ans_num;i++) Write(ans[i].first),putchar(' '),Write(ans[i].second),putchar('\n');
		return;
	}
}
int main(){
	if (!(R=read())){printf("1\n0 0");return 0;} 
	Pollard_Rho::divide(R);sort(P+1,P+1+P[0]);
	P[0]=unique(P+1,P+1+P[0])-(P+1);Get_ans::solve();
	return 0;
}