首先是一个结论,若 \(n=2^\alpha\prod p_i^{\beta_i}\prod q_i^{2\gamma_i}\) ,其中 \(p_i\equiv1 \pmod 4,q_i\equiv 3 \pmod 4\) ,那么有解,否则无解。
然后考虑如果 \(n_1=x_1^2+y_1^2,n_2=x_2^2+y_2^2\) ,那么 \(n_1n_2\) 也是可以被表示的,即 \(n_1n_2=(x_1y_1+x_2y_2)^2+(x_1y_2-x_2y_1)^2\) 。
注意到 \(2^\alpha\) 是好分解的,按照 \(\alpha\) 的奇偶性讨论即可,对于 \(\beta_i,2\gamma_i\) 其实可以直接把偶数次方拆出来最后再乘,那么现在就是要分解 \(p\equiv 1\pmod 4\) 的质数了。
考虑如果构造了 \(x^2+y^2=mp\) ,那么根据这个构造一组 \(a^2+b^2=rp\) ,其中 \(r<m\) ,不断重复就可以得到一组解。
考虑构造一组 \((u,v)\) 满足 \(u\equiv x\pmod m,v\equiv y\pmod m\) ,那么有 \(x^2+y^2\equiv u^2+v^2\pmod m\) ,不妨直接令 \(u^2+v^2=mr\) 。
把两个平方和的式子相加可以得到 \((xu+yv)^2+(xv-yu)^2=m^2rp\) ,那么根据同余的性质,可以得到 \((\frac{xu+yv}{m})^2+(\frac{xy-yu}{m})^2=rp\) ,取 \(\mid u,v\mid\leq \frac{m}{2}\) ,可以使得 \(r\leq \frac{m}{2}\) ,这样就能 \(O(\log V)\) 的构造了。
注意到一开始要构造一组 \(x^2+y^2=mp\) ,可以想到令 \(x^2\equiv -1 \pmod p,y=1\) ,那么直接随机 \(a\) ,判断 \(a^{\frac{p-1}{2}}\) 是否为 \(-1\) ,根据二次剩余的结论,这个随机几乎是 \(O(1)\) 的,得到 \(a\) 之后只需要令 \(x\equiv a^{\frac{p-1}{4}}\pmod p\) 即可。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 I;
random_device RD;
mt19937 rnd(RD());
int n;
inline ll add(ll a,ll b,ll p){a+=b;return a>=p?a-p:a;}
inline ll ksc(ll a,ll b,ll p){ll x=(long double)a/p*b;return ((I)a*b-(I)x*p+p)%p;}
inline ll ksm(ll a,ll b,ll p){
ll res=1;
a%=p;
while (b){
if (b&1) res=ksc(res,a,p);
a=ksc(a,a,p);
b>>=1;
}
return res;
}
namespace MillerRabin{
const int P=15;
const int pr[P]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,19260817};
inline bool check(ll x,ll p){
if (x%p==0 || ksm(p%x,x-1,x)!=1) return 0;
ll a=x-1;
while (!(a&1)){
ll t=ksm(p%x,a>>=1,x);
if (t!=1 && t!=x-1) return 0;
if (t==x-1) return 1;
}
return 1;
}
inline bool Isprime(ll x){
if (x>2 && x%2==0) return 0;
if (x<=2) return x==2;
for (int i=0;i<P;++i){
if (x==pr[i]) return 1;
if (!check(x,pr[i])) return 0;
}
return 1;
}
}
namespace PollardRho{
vector<ll> fac;
ll gcd(ll x,ll y){return !y?x:gcd(y,x%y);}
void clear(){fac.clear();}
void unique(){
sort(fac.begin(),fac.end());
auto it=unique(fac.begin(),fac.end());
fac.erase(it,fac.end());
}
ll c;
inline ll f(ll x,ll p){return ((I)x*x+c)%p;}
ll pollard_rho(ll x){
if (x==4) return 2;
while (1){
c=uniform_int_distribution<ll>(1,x-1)(rnd);
ll s=f(f(0,x),x),t=f(0,x);
while (s!=t){
ll val=1;
for (int i=0;i<128 && s!=t;++i){
ll tval=(I)val*abs(s-t)%x;
if (!tval) break;
val=tval;
s=f(f(s,x),x),t=f(t,x);
}
if (val){
ll d=gcd(val,x);
if (d>1) return d;
}
}
}
}
void factor(ll x){
if (x<2) return;
if (MillerRabin::Isprime(x)){
fac.emplace_back(x);
return;
}
ll p=pollard_rho(x);
while ((x%p)==0) x/=p;
factor(x);factor(p);
}
vector<pair<ll,int>> get_factor(ll n){
clear(),factor(n),unique();
vector<pair<ll,int>> res(0);
for (ll p:fac){
int c=0;
while (n%p==0) ++c,n/=p;
res.emplace_back(p,c);
}
return res;
}
}
struct node{ll x,y;node(ll a=0,ll b=0){x=a,y=b;}};
inline node operator + (node a,node b){return node(a.x*b.x+a.y*b.y,abs(a.x*b.y-a.y*b.x));}
ll pow(ll x,ll y){
ll res=1;
while (y){
if (y&1) res*=x;
x*=x;
y>>=1;
}
return res;
}
ll find(ll p){// 4k+1
ll x=uniform_int_distribution<ll>(1,p-1)(rnd);
while (ksm(x,(p-1)/2,p)==1) x=uniform_int_distribution<ll>(1,p-1)(rnd);
return ksm(x,(p-1)/4,p);
}
I myabs(I x){return x<0?-x:x;}
void work(){
ll n;scanf("%lld",&n);
vector<pair<ll,int>> pr=PollardRho::get_factor(n);
node ans(1,0);
for (auto x:pr){
ll p=x.first;int c=x.second;
if (p==2){
if (c&1) ans=ans+node(pow(2,c/2),pow(2,c/2));
else ans=ans+node(pow(2,c/2),0);
}
else if (p%4==3){
if (c&1) return puts("-1"),void();
ans.x*=pow(p,c/2),ans.y*=pow(p,c/2);
}
else{
ans.x*=pow(p,c/2),ans.y*=pow(p,c/2);
if (c&1){
I x=find(p),y=1;
while (x*x+y*y!=p){
I k=(x*x+y*y)/p;
I u=x%k,v=y%k;
u-=(2*u>k)*k,v-=(2*v>k)*k;
tie(x,y)=make_pair((u*y-v*x)/k,(u*x+v*y)/k);
x=myabs(x),y=myabs(y);
}
ans=ans+node(x,y);
}
}
}
printf("%lld %lld\n",ans.x,ans.y);
}
int main(){
int T;scanf("%d",&T);
while (T--) work();
return 0;
}
标签:分成,node,平方,return,int,res,ll,整数,ans
From: https://www.cnblogs.com/jerryjiang/p/17472551.html