题意：给定非负整数\(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;
}