NOD 1171 x^2+y^2=N

(1)先找出同余方程的最小正整数解。

(2)对和进行欧几里德辗转相除运算，记每次的余数为，当满足条件时停止运算，此时的就是x

这样就得到了x，那么y可以通过得到，问题解决。本方法是目前为止发现解此问题最快的方法。

所以将问题转化为求同余方程的最小正整数解。

http://algo.ftiasch.com/tag/number-theory/ 粘过来有些奇怪 膜拜一下叉姐

今天的话题是解方程x2≡n(modp)，其中p是奇质数。

引理：ap−12≡±1(modp)

证明：由费马小定理，ap−1−1≡(ap−12−1)(ap−12+1)≡0⇒ap−12≡±1

引理：方程有解当且仅当np−12≡1(modp)。

（充分性）xp−1≡1⇒np−12≡(x2)p−12≡1

（必要性）如果方程无解，则可以把1,2,…,(p−1)配成p−12对，每对乘积为n，则np−12≡(p−1)!≡−1(modp)。最后的等号是威尔逊定理。

设a满足ω=a2−n不是二次剩余（即x2≡ω(modp)无解）

定理：x≡(a+ω√)p+12是方程x2≡n(modp)的解。

证明： (a+ω√)p=ap+ωp−12ω√=a−ω√，前面的等号用了二项式定理和(pi)≡0(modp)，后面的等号用了费马小定理和ω不是二次剩余。

(a+ω√)p+1≡(a+ω√)p(a+ω√)≡(a−ω√)(a+ω√)≡a2−ω≡n

具体实现的时候，a的选择可以随机，因为大约有一半非二次剩余，剩下的只有快速幂而已。

根据以上叉姐博客里的结论可以得出也就是通俗的讲随机(0-p-1)选取一个随机数满足x^2=-1(mod p)无解也就是a是素数p的二次非剩余即成立。对二次域(a+sqrt(w))^(p+1)%p即为x0.

x0求出后在通过对和进行欧几里德辗转相除运算，记每次的余数为，当满足条件时停止运算，此时的就是x，这样就得到了x，那么y可以通过。

具体见程序注释。

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<set>
using namespace std;
typedef long long ll;
ll gcd(ll a,ll b)
{
    return b?gcd(b,a%b):a;
}
ll Mulmod(ll a,ll b,ll n)
{
    ll  exp=a%n,res=0;
    while(b)
    {
        if(b&1)
        {
            res += exp;
            if(res>n) res -= n;
        }
        exp<<=1,b>>=1;
        if(exp>n) exp -= n;
    }
    return res;
}
ll exp_mod(ll a,ll b,ll c)
{
    ll k=1;
    while(b)
    {
        if(b&1) k=Mulmod(k,a,c);
        a=Mulmod(a,a,c),b>>=1;
    }
    return k;
}
bool Miller_Rabbin(ll n, ll times)   
{
    if(n==2) return 1;
    if(n<2||!(n&1)) return 0;
    ll a,u=n-1,x,y;
    int t=0;
    while(u%2==0) t++,u/=2;
    srand(100);
    for(int i=0; i<times; i++)
    {
        a=rand()%(n-1)+1,x=exp_mod(a,u,n);
        for(int j=0; j<t; j++)
        {
            y=Mulmod(x,x,n);
            if(y==1&&x!=1&&x!=n-1) return false; //must not
            x=y;
        }
        if(y!=1) return false;
    }
    return true;
}
ll Pollard_Rho(ll n,ll c)
{
    ll x,y,d,i=1,k=2;
    y=x=rand()%(n-1)+1;
    while(1)
    {
        i++;
        x=(Mulmod(x,x,n)+c)%n,d=gcd((x-y+n)%n,n);
        if(d>1&&d<n) return d;
        if(x==y) return n;
        if(i==k) k<<=1,y=x;
    }
}
ll factor[100],tol,col,fac[100][2];
void Find_factor(ll n,ll c)    //大数分解 fac[i][0]存素因子 fac[i][1]存幂
{
    if(n==1) return;
    if(Miller_Rabbin(n,10))
    {
        factor[tol++]=n;
        return;
    }
    ll p=n,k=c;
    while(p>=n) p = Pollard_Rho(p,c--);
    Find_factor(p,k);
    Find_factor(n/p,k);
}
void getfactor() //把素因子转化为fac素组
{
    sort(factor,factor+tol);
    memset(fac,0,sizeof(fac));
    col=0,fac[col][0]=factor[0],fac[col][1]++;
    for(int i=1; i<tol; i++)
        if(factor[i]!=factor[i-1]) fac[++col][0]=factor[i],fac[col][1]++;
        else fac[col][1]++;
    col++;
}
bool judge() //判断是否存在4k+3为奇次幂的
{
    for(int i=0; i<col; i++)
        if(fac[i][0]%4==3&&fac[i][1]&1)
        {
            puts("No Solution");
            return 1;
        }
    return 0;
}
ll sfac[100][2],num;//存的是每个素因子拆分后的x,y
struct qua //二次域的结构体
{
    ll a,b;
};
ll w;//根号下的w
qua mul_quamod(qua a,qua b,ll n)
{
    qua ans;
    ans.a=(Mulmod(a.a,b.a,n)+Mulmod(Mulmod(a.b,b.b,n),w,n))%n;
    ans.b=(Mulmod(a.a,b.b,n)+Mulmod(b.a,a.b,n));
    return ans;
}
qua exp_quamod(qua a,ll b,ll n) //二次域的高次幂取模 原理跟整数的一样
{
    qua ans;
    ans.a=1,ans.b=0;
    while(b)
    {
        if(b&1) ans=mul_quamod(ans,a,n);
        b>>=1,a=mul_quamod(a,a,n);
    }
    return ans;
}
void getsfac(ll p,ll &x,ll &y)//得到4k+1型素因子的x,y
{
    ll a,m,t=(ll)sqrt(p),x0,r;
    while(1)        //随机数求1-sqrt(p)内的一个符合素数p的非二次剩余 这里到根号p是防止a*a超long long
    {
        a=rand()%t,m=a*a+1;
        if((exp_mod(m,(p-1)>>1,p)+1)%p==0) break;
    }
    w=a*a+1; 
    qua qu; //定义二次域求x0
    qu.a=a,qu.b=1;
    qu=exp_quamod(qu,(p+1)>>1,p);
    x0=qu.a,t=p,r=p%x0;     //辗转相处求x
    while(r>=(ll)sqrt(p))
        t=x0,x0=r,r=t%x0;
    x=r,y=(ll)sqrt(p-r*r);  //通过x求y 
}
void solve() //针对每种不同的素因子求得出x^2+y^2=p
{
    num=0;
    for(int i=0; i<col; i++)
        for(int j=0; j<fac[i][1]; j++,num++)
            if(fac[i][0]==2) sfac[num][0]=sfac[num][1]=1;
            else if(fac[i][0]%4==1) getsfac(fac[i][0],sfac[num][0],sfac[num][1]);
            else if(fac[i][0]%4==3) sfac[num][0]=0,sfac[num][1]=fac[i][0],j++;
}
struct finans //定义最后的结果 
{
    ll x,y;
    bool operator<(const finans &o)const
    {
        return x<o.x;
    }
} data;
set<finans>wans[2];
set<finans>::iterator it;
int f1,f2;
void getans() //用了两个set 滚动存放每个素因子通过费马平方和后心的结果 用set的目的是去重
{
    wans[0].clear(),wans[1].clear();
    f1=1,f2=0;
    data.x=sfac[0][0],data.y=sfac[0][1];
    if(data.x>data.y)swap(data.x,data.y);
    wans[f1].insert(data);
    for(int i=1; i<num; i++)
    {
        swap(f1,f2);
        wans[f1].clear();
        for(it=wans[f2].begin(); it!=wans[f2].end(); it++)
        {
            data.x=(ll)fabs(it->x*sfac[i][0]+it->y*sfac[i][1]),
                   data.y=(ll)fabs(it->x*sfac[i][1]-it->y*sfac[i][0]);
            if(data.x>data.y)swap(data.x,data.y);
            wans[f1].insert(data);
            data.x=(ll)fabs(it->x*sfac[i][0]-it->y*sfac[i][1]),
                   data.y=(ll)fabs(it->x*sfac[i][1]+it->y*sfac[i][0]);
            if(data.x>data.y)swap(data.x,data.y);
            wans[f1].insert(data);
        }
    }
}
int main()
{
    ll n;
    while(~scanf("%I64d",&n))
    {
        if(n==1)
        {
            puts("0 1");
            continue;
        }
        tol=0,Find_factor(n,120);
        getfactor();
        if(judge()) continue;
        solve();
        getans();
        for(it=wans[f1].begin(); it!=wans[f1].end(); it++)
            printf("%I64d %I64d\n",it->x,it->y);
    }
    return 0;
}