一个数N(1 <= N <= 10^18)
共K行:每行2个数,i j,表示N = i^2 + j^2(0 <= i <= j)。 如果无法分解为2个数的平方和,则输出No Solution
(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;
}