BZOJ 1041
题目描述
给出\(n\),求\(x^2+y^2=n^2,x,y,z\in \mathbb{Z}\)的解数.
复杂度
\(O\left(T_{\mathtt{factorization}}(n)\right)\).在解答中使用Pollard \(\rho\)算法,复杂度\(O\left( n^{1/4}\right)\).
题解
见 潘承洞、潘承彪 《初等数论》第三版 第287页的结论,可简单推出
对于\(n\)的一个形如\(4n+1\)的素因子\(p\),对答案的贡献是\(2c+1\),其中\(c\)是\(n\)中\(p\)的指数.(对于每个\(n^2\)的因子,\(p\)的指数有\(2c+1\)总选择)
对于\(n\)的每个形如\(4n+3\)的素因子\(p_i\),总共有\(t\)个,那么对于每一种每个\(c_i\)(对于\(4n+3\)素数\(p_i\)的次数)的和为奇数的选择方法,有一个\(c_i\)的和为偶数的选择与它一一对应.举一个例子:
对于一个选择\(c_1\dots c_t\),如果\(\sum_{i=1}^tc_i\bmod 2=1\),那么有奇数个\(c_i\)是奇数.如果\(c_1\dots c_{i-1}\)都对应的等于\(max\{c_1\}\dots max\{c_{i-1}\}\)(即\(n\)中\(p_{1\dots i-1}\)的次数\(C_{1\dots i-1}\)的两倍),注意这个数是奇数,而\(c_i\not= 2C_i\),这时可以把\(c_i\)加一形成一个和为偶数的选择,容易看出这个是单射, 而对于\(c_1\dots c_t\)都为0的选择,没有对应的和为奇数的选择与它对应,那么这些素数总贡献就是\(1\).
所以答案就是对于\(4n+1\)素因子求他们分别的\(2c_i+1\)的积乘4.
代码
#include <map>
#include <cmath>
#include <cstdio>
#include <vector>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#define mod 1000000007
#define maxint 2147483648ll
#define nc1
#define rg register
#define rgv
#define rg1 rgv
#define rg2 rgv
#define rg3 rgv
#define rg4 rgv
#define rg5 rg
#define rg6 rg
#define sorc
#define psh(a,l,p) (a)[(l)++]=(p)
typedef long long ll;
namespace basicmath{
inline ll mp(ll a,ll b,ll p){
#ifdef __int128
return (__int128)a*(__int128)b%p;
#else
if(p<3037000000ll) return a*b%p;
#ifdef nc1
if(p<1099511627776ll) return (((a*(b>>20)%p)<<20)+(a*(b&((1<<20)-1))))%p;
#endif
ll d=(ll)(a*(long double)b/p+0.5); ll ret=(a*b-d*p)%p;
if (ret<0) ret+=p;
return ret;
#endif
}
inline ll fp(int a,int b,int p){ a%=p,b%=p-1; int ans=1; for(;b;b>>=1,a=(ll)a*a%p) if(b&1) ans=(ll)ans*a%p; return ans; }
inline ll fp(ll a,ll b,ll p){ a%=p,b%=p-1; ll ans=1ll; for(;b;b>>=1,a=mp(a,a,p)) if(b&1) ans=mp(a,ans,p); return ans; }
template<int max=10000007> struct sieve{
int q[max+5]; int pr[max/4],pl,rxsiz;
inline void generate(rg3 int n=max){
if(n<=rxsiz) return; if(n>max) return;
rxsiz=n;
for(rg1 int i=2;i<=n;++i){
if(!q[i]) pr[pl++]=i,q[i]=i;
for(rg2 int j=0;j<pl;++j){
rg4 int t=i*pr[j];
if(t>n) break;
q[t]=pr[j]; if(q[i]==pr[j]) break;
}
}
}
inline int operator()(int n){
#ifdef sorc
if(n>=pl) return 0;
if(n<0) return 0;
#endif
return pr[n];
}
#ifdef sorc
int sbcc;
inline int operator[](int n){
if(n>rxsiz) return sbcc;
if(n<0) return sbcc;
return q[n];
}
#else
inline int operator[](int n){
return q[n];
}
#endif
};
namespace primeTesting{
const int llrtm=7;
ll q[llrtm]={2, 325, 9375, 28178, 450775, 9780504, 1795265022ll};//for numbers less than 2^64 test these seven numbers are enough
inline bool witness(ll a,rg5 ll n){
a%=n;
if(!a) return 0;
rg5 int t=0; ll u=n-1; for(;~u&1;u>>=1) ++t;
rg5 ll x=fp(a,u,n),xx=0; while(t--){
xx=mp(x,x,n); if(xx==1 && x!=1 && x!=n-1) return 1;
x=xx;
}
return xx!=1;
}
inline bool miller(ll n){
if(!n) return 0;
if(n==2) return 1;
for(int i=0;i<llrtm;++i) if(witness(q[i],n)) return 0;
return 1;
}
template<class T> inline bool miller(ll n,T& sv){//T should have the same API as sieve
if(n<=sv.rxsiz && sv[n]==n) return 1;
return miller(n);
}
//the both miller tests are deterministic
struct pmt{
virtual inline bool operator()(ll n){
return 1;
}
};
struct mr_s100k:public pmt{
sieve<100000> sv;
inline mr_s100k(){
sv.generate(100000);
}
inline bool operator()(ll n){
return miller(n,sv);
}
};
struct mr_0:public pmt{
inline bool operator()(ll n){
return miller(n);
}
};
}
using namespace primeTesting;
namespace gcd{
inline ll euclid(ll a,ll b){return a?euclid(b%a,a):b;}
inline ll stein (rg6 ll a,rg6 ll b){
// if(a<0) a=-a;
rg6 ll ret=1;
while(a){
if((~a&1)&&(~b&1)) ret<<=1,a>>=1,b>>=1; else
if(~a&1) a>>=1; else if(~b&1) b>>=1; else
{ if(a<b) std::swap(a,b); a-=b; }
}
return b*ret;
}
}
namespace exgcd{
inline ll euclid(ll a,ll b,ll& x,ll& y){
if(b){
ll d=euclid(b,a%b,x,y);
ll t=x; x=y; y=t-a/b*y;
return d;
}else{
x=1,y=0;
return a;
}
}
}
namespace qmutils{
inline ll qm1(rg ll a,rg ll q){//if a < 0 a+=q
return a+((-((unsigned long long)a>>63))&q);
}
inline ll qm2(rg ll a,rg ll q){//if a>= q a-=q
return a-((-((unsigned long long)(q+(~a))>>63))&q);
}
}
using namespace qmutils;
/*
namespace ec{
using namespace exgcd;
struct pi{
ll x,y,t;
inline pi(ll x=0,ll y=0,ll t=0):x(x),y(y),t(t){}
inline pi(const pi& p):x(p.x),y(p.y),t(p.t){}
};
inline bool operator==(const pi& a,const pi& b){
return a.x==b.x && a.y==b.y;
}
int sd[20];
std::mt19937_64 tr;
inline void mtp(){
std::random_device kkk;
for(int i=0;i<20;++i) sd[i]=kkk();
std::seed_seq pps(sd,sd+20);
tr.seed(pps);
}
inline ll mtr(ll f,ll t){
std::uniform_int_distribution<ll> kk(f,t);
return kk(tr);
}
#define randx mtr
struct cv{//determines a curve
ll a,b;
ll x,y,p;
inline ll cb(ll x){ return mp(mp(x,x,p),x,p); }
inline pi getpi(){
return pi(x,y,0);
}
inline void random(ll px){
p=px;
x=randx(1,p-1); y=randx(1,p-1);
a=randx(1,p-1);
b=qm1(qm1(mp(y,y,p)-cb(x),p)-mp(x,a,p),p);
}
inline ll slope(pi a,pi b,ll& f){
if(a==b){
ll yy=(3*mp(a.x,a.x,p))%p,xx=qm2(a.y+a.y,p);
ll xp,yp;
f=exgcd::euclid(p,xx,xp,yp);
yp=qm2(qm1(yp,p),p);
return mp(yy,yp,p);
}
ll yy=qm1(a.y-b.y,p),xx=qm1(a.x-b.x,p);
ll xp,yp;
f=exgcd::euclid(p,xx,xp,yp);
yp=qm2(qm1(yp,p),p);
return mp(yy,yp,p);
}
inline pi add (pi a,pi b,ll& f){
if(a.t){
f=1;
return b;
}
if(b.t){
f=1;
return a;
}
if(a.x==b.x && a.y==-b.y){
f=1;
return pi(0,0,1);
}
ll sl=slope(a,b,f);
ll xp=qm1(qm1(mp(sl,sl,p)-a.x,p)-b.x,p);
return pi(
// x: sl^2 - xp - xq
xp ,
// y: yp + sl * (xR - xp)
p-qm2(a.y+mp(sl,qm1(xp-a.x,p),p),p) ,
0
);
}
inline pi mult(pi a,ll b,ll& f){
f=1;
pi t=pi(0,0,1);
for(;b;b>>=1,a=add(a,a,f)){
if(f!=1) return t;
if(b&1) t=add(t,a,f);
if(f!=1) return t;
}
}
inline pi mults(pi a,ll b,ll& f){
f=1;
pi t(a);
for(int i=1;i<b;++i){
t=add(t,a,f);
if(f!=1) return t;
}
return t;
}
};
}
*/
#define dgcd gcd::euclid
namespace factor{
struct factr{
ll q[30],p[30],c[30];
int tl;
virtual inline int operator()(ll a,pmt& q){
return 0;
}
};
struct rho:public factr{
ll sq[70],sql;
private:
ll rhox(ll n,ll c,ll u){
ll i=1,k=2,y=u,x0=u;
while(1){
++i; x0=qm2(mp(x0,x0,n)+c,n);
ll d=dgcd(y-x0,n);
if(d!=1 && d!=n) return d;
if(y==x0) return n;
if(i==k) y=x0,k<<=1;
}
}
void fact_f(ll a,pmt& q){
if(q(a)) psh(sq,sql,a); else
{
ll x=a;
while(x==a){
ll c=rand()%a,d=(rand()&31)+3;
x=rhox(a,d,c);
}
fact_f(x,q),fact_f(a/x,q);
}
}
inline void norm(){
std::sort(sq,sq+sql);
tl=0;
for(rg int i=0,j=1;i<sql;i=(j++)){
q[tl]=sq[i],p[tl]=1,c[tl]=sq[i];
while(sq[j]==sq[i] && j<sql) c[tl]*=sq[j++],++p[tl];
++tl;
}
}
public:
inline int operator()(ll a,pmt& q){
sql=0;
fact_f(a,q);
norm();
return tl;
}
};
template<int nx=1000000> struct rhos:public factr{
ll sq[70],sql;
sieve<nx> sv;
private:
ll rhox(ll n,ll c,ll u){
ll i=1,k=2,y=u,x0=u;
while(1){
++i; x0=qm2(mp(x0,x0,n)+c,n);
ll d=dgcd(y-x0,n);
if(d<0) d=-d;
if(d!=1 && d!=n) return d;
if(y==x0) return n;
if(i==k) y=x0,k<<=1;
}
}
void fact_f(ll a,pmt& q){
if(a<=nx){//directly factor a
while(a!=1){
psh(sq,sql,sv[a]);
a/=sv[a];
}
return;
}
if(q(a)) psh(sq,sql,a); else
{
ll x=a;
while(x==a){
ll c=rand()%a,d=(rand()&31)+3;
x=rhox(a,d,c);
}
fact_f(x,q),fact_f(a/x,q);
}
}
inline void norm(){
std::sort(sq,sq+sql);
tl=0;
for(rg int i=0,j=1;i<sql;i=(j++)){
q[tl]=sq[i],p[tl]=1,c[tl]=sq[i];
while(sq[j]==sq[i] && j<sql) c[tl]*=sq[j++],++p[tl];
++tl;
}
}
public:
rhos(){
sv.generate(nx);
}
inline int operator()(ll a,pmt& q){
sql=0;
fact_f(a,q);
// for(int i=0;i<sql;++i) printf("%lld ",sq[i]);
norm();
return tl;
}
};
}
using namespace factor;
struct factorsGetter{
//get all factors
std::vector<ll> fac;
void get(ll q,int le,factr& p){
if(le==p.tl) fac.push_back(q); else{
get(q,le+1,p);
for(int i=1;i<=p.p[le];++i) get((q*=p.q[le]),le+1,p);
}
}
int operator()(factr& p){
fac.clear();
get(1,0,p);
return fac.size();
}
};
}
using namespace basicmath;
rhos<> rhoer;
mr_0 pp;
#define maxf 3
const int pxf[10][2]={{200,1},{2000,1},{10000,2}};
int main(){
ll a;
srand(32345);
scanf("%lld",&a);
int kk=rhoer(a,pp);
ll ans=1;
for(int i=0;i<kk;++i){
// printf("%lld^%lld=%lld\n",rhoer.q[i],rhoer.p[i],rhoer.c[i]);
// fflush(stdout);
if((rhoer.q[i]&3)==1){
ans*=rhoer.p[i]<<1|1;
}
}
// if(a>1 && a&3==1) ans*=3;
printf("%lld\n",ans<<2);
return 0;
}
