1171 . 两个数的平方和 V2
时间限制：1 秒 空间限制：65536 KB 分值: 1280
给出一个整数N，将N表示为2个整数i j的平方和（i <= j)，如果有多种表示，按照i的递增序输出。
例如：N = 130，130 = 3^2 + 11^2 = 7^2 + 9^2 （注：3 11同11 3算1种）
Input
一个数N(1 <= N <= 10^18)
Output
共K行：每行2个数，i j，表示N = i^2 + j^2（0 <= i <= j)。 如果无法分解为2个数的平方和，则输出No Solution
题解：这题非常全面。
首先 如果一个奇素数p为4k+1形式那么有且仅有一组x,y可以表示为x^2+y^2=p。
当然2=1^2+1^2。
那么根据费马平方和定理，如果两个数可以表示成两数平方和的形式那么着两个数的乘积也能表示长平方和的形式，即
。
。
首先如果一个数表示成素因子乘积的形式那么4k+3型的素因子p的幂a一定是偶数。偶数幂可以写成
如果N中含有素因子p的幂为奇数，那么直接可以输出no solution了。
接下来就是如何将4k+1形式的素数分解成两个数x,y的平方和了。至于x,y如何求就要分两步。
(1)先找出同余方程
(2)对
这样就得到了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的二次非剩余即
x0求出后在通过对
具体见程序注释。
#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;
}
ap−12≡±1(modp)
