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;
}
posted @ 2015-12-03 15:58  zball  阅读(130)  评论(0编辑  收藏