【扩展欧几里德算法】
〖相关资料〗
《欧几里德算法的证明》 《欧几里德算法与扩展欧几里德算法》 《扩展欧几里德专题》
〖相关知识〗
对于非负整数 $a,b$，求数对 $x,y$ ，满足 $ax+by=gcd(a,b)$ 。
对于不定整数方程 $ax+by=c$ ，若 $cmod gcd(a,b)=0$ ，则该方程存在整数解，否则不存在整数解。
若已知 $ax+by=gcd(a,b)$ 的一组解 $x_{0},y_{0}$ ，则其他整数解满足 $x=x_{0}+frac{b}{gcd(a,b)}cdot t$ ， $y=y_{0}-frac{a}{gcd(a,b)}cdot t$ 。
出现负数时取绝对值。
〖模板代码〗
[扩展欧几里得算法]
1 LL exgcd(LL a,LL b,LL& x,LL& y) 2 { 3 if(!b){x=1;y=0;return a;} 4 LL ans=exgcd(b,a%b,x,y); 5 LL tmp=x;x=y;y=tmp-a/b*y; 6 return ans; 7 }
[扩欧求逆元]
1 LL exgcd(LL a,LL b,LL& x,LL& y) 2 { 3 if(!b){x=1;y=0;return a;} 4 LL ans=exgcd(b,a%b,x,y); 5 LL tmp=x;x=y;y=tmp-a/b*y; 6 return ans; 7 } 8 LL inv(LL t,LL mod) 9 { 10 if(!t)return 0; 11 LL a=t,b=mod,x,y,g=exgcd(a,b,x,y); 12 if(g!=1)return -1; 13 x=(x%b+b)%b;if(!x)x=b; 14 return x; 15 }
〖相关题目〗
1.【bzoj1477】青蛙的约会
题意：两只青蛙在同一条纬度线上，各自朝西跳，直到碰面为止。除非这两只青蛙在同一时间跳到同一点上，不然是永远都不可能碰面的。规定纬度线上东经0度处为原点，由东往西为正方向，单位长度1米，这样我们就得到了一条首尾相接的数轴。设青蛙A的出发点坐标是x，青蛙B的出发点坐标是y。青蛙A一次能跳m米，青蛙B一次能跳n米，两只青蛙跳一次所花费的时间相同。纬度线总长L米。现在要你求出它们跳了几次以后才会碰面。
分析：TJMの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 LL x,y,X,Y,m,n,L,a,b,c,d; 7 LL exgcd(LL a,LL b,LL& x,LL& y) 8 { 9 if(!b){x=1;y=0;return a;} 10 LL ans=exgcd(b,a%b,x,y); 11 LL tmp=x;x=y;y=tmp-a/b*y; 12 return ans; 13 } 14 int main() 15 { 16 scanf("%lld%lld%lld%lld%lld",&X,&Y,&m,&n,&L); 17 a=m-n;b=L;c=Y-X;d=exgcd(a,b,x,y); 18 if(c%d){printf("Impossible");return 0;} 19 x=x*(c/d);b/=d;b=b>0?b:-b; 20 x=((x%b)+b)%b;x=x?x:b; 21 printf("%lld",x); 22 return 0; 23 }
2.【poj2142】The Balance
题意：有两个类型的砝码,质量分别为a,b；现在要求称出质量为d的物品，要用多少a砝码(x)和多少b砝码(y)，使得(x+y)最小。
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 LL A,B,d,g,x1,Y1,x2,y2; 7 LL exgcd(LL a,LL b,LL& x,LL& y) 8 { 9 if(!b){x=1;y=0;return a;} 10 LL ans=exgcd(b,a%b,x,y); 11 LL tmp=x;x=y;y=tmp-a/b*y; 12 return ans; 13 } 14 void work(LL a,LL b,LL& x,LL& y) 15 { 16 g=exgcd(a,b,x,y); 17 x*=d/g;LL t=b/g;t=t>0?t:-t; 18 x=(x%t+t)%t;y=(a*x-d)/b>0?(a*x-d)/b:(-(a*x-d)/b); 19 } 20 int main() 21 { 22 while(1) 23 { 24 scanf("%lld%lld%lld",&A,&B,&d); 25 if(A==0&&B==0&&d==0)break; 26 work(A,B,x1,Y1);work(B,A,y2,x2); 27 if(x1+Y1<x2+y2)printf("%lld %lld ",x1,Y1); 28 else printf("%lld %lld ",x2,y2); 29 } 30 return 0; 31 }
3.【LOJ6257】「CodePlus 2017 12 月赛」可做题2
题意：若一个数列a满足条件an=an−1+an−2，n≥3，而a1,a2为任意实数，则我们称这个数列为广义斐波那契数列。现在请你求出满足条件a1=i，a2为区间[l,r]中的整数，且akmodp=m的广义斐波那契数列有多少个。
分析：矩阵快速幂+exgcd计算周期。需要推公式。
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 LL T,A,L,R,K,P,M,g,x,y,ans; 7 struct node 8 { 9 LL m[3][3]; 10 void clear(){memset(m,0,sizeof(m));} 11 void init(){m[1][1]=m[2][2]=1;m[1][2]=m[2][1]=0;} 12 void base(){m[1][1]=m[1][2]=m[2][1]=1;m[2][2]=0;} 13 }bas,tmp; 14 LL read() 15 { 16 LL x=0,f=1;char c=getchar(); 17 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 18 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 19 return x*f; 20 } 21 node operator * (node a,node b) 22 { 23 node c;c.clear(); 24 for(int i=1;i<=2;i++) 25 for(int j=1;j<=2;j++) 26 for(int k=1;k<=2;k++) 27 c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%P)%P; 28 return c; 29 } 30 node qm(node a,LL b) 31 { 32 node ans;ans.init(); 33 while(b) 34 { 35 if(b&1)ans=ans*a; 36 a=a*a;b>>=1; 37 } 38 return ans; 39 } 40 LL exgcd(LL a,LL b,LL& x,LL& y) 41 { 42 if(!b){x=1;y=0;return a;} 43 LL ans=exgcd(b,a%b,x,y); 44 LL tmp=x;x=y;y=tmp-a/b*y; 45 return ans; 46 } 47 void work() 48 { 49 A=read();L=read()-1;R=read();K=read();P=read();M=read(); 50 bas.base();tmp=qm(bas,K-2); 51 A=A%P;M=((M%P-tmp.m[2][1]*A%P)%P+P)%P; 52 g=exgcd(tmp.m[1][1],P,x,y); 53 if(M%g){printf("0 ");return;} 54 x=(x%P+P)%P;x=x*(M/g)%P;P/=g;x%=P; 55 ans=(R<x?0:(R-x)/P+1); 56 ans-=(L<x?0:(L-x)/P+1); 57 printf("%lld ",ans); 58 } 59 int main() 60 { 61 T=read(); 62 while(T--)work(); 63 return 0; 64 }
4.【poj2891】Strange Way to Express Integers
题意：解同余方程组，模数不一定互质。
分析：zmhの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=1e5+5; 7 LL n,a[N],m[N]; 8 int read() 9 { 10 int x=0,f=1;char c=getchar(); 11 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 12 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 13 return x*f; 14 } 15 LL exgcd(LL a,LL b,LL& x,LL& y) 16 { 17 if(!b){x=1;y=0;return a;} 18 LL ans=exgcd(b,a%b,x,y); 19 LL tmp=x;x=y;y=tmp-a/b*y; 20 return ans; 21 } 22 LL work() 23 { 24 LL M=m[1],A=a[1],x,y,g; 25 for(int i=2;i<=n;i++) 26 { 27 g=exgcd(M,m[i],x,y); 28 if((A-a[i])%g)return -1; 29 x=(A-a[i])/g*x%m[i]; 30 A-=x*M;M=M/g*m[i];A%=M; 31 } 32 return (A+M)%M; 33 } 34 int main() 35 { 36 while(scanf("%lld",&n)==1) 37 { 38 for(int i=1;i<=n;i++)m[i]=read(),a[i]=read(); 39 printf("%lld ",work()); 40 } 41 return 0; 42 }
5.【hdu1573】x问题
题意：求在小于等于N的正整数中有多少个X满足：X mod a[0] = b[0], X mod a[1] = b[1], X mod a[2] = b[2], …, X mod a[i] = b[i], … (0 < a[i] <= 10)。
分析：拦路雨偏似雪花の博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=15; 7 LL T,n,cnt,ans,m[N],a[N]; 8 LL read() 9 { 10 LL x=0,f=1;char c=getchar(); 11 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 12 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 13 return x*f; 14 } 15 LL exgcd(LL a,LL b,LL& x,LL& y) 16 { 17 if(!b){x=1;y=0;return a;} 18 LL ans=exgcd(b,a%b,x,y); 19 LL tmp=x;x=y;y=tmp-a/b*y; 20 return ans; 21 } 22 void work() 23 { 24 n=read();cnt=read(); 25 for(int i=1;i<=cnt;i++)m[i]=read(); 26 for(int i=1;i<=cnt;i++)a[i]=read(); 27 LL M=m[1],A=a[1],x,y,g; 28 for(int i=2;i<=cnt;i++) 29 { 30 g=exgcd(M,m[i],x,y); 31 if((A-a[i])%g){printf("0 ");return;} 32 x=(A-a[i])/g*x%m[i]; 33 A-=x*M;M=M/g*m[i];A%=M; 34 } 35 A=(A+M)%M; 36 if(n<A)printf("0 "); 37 else 38 { 39 ans=(n-A)/M+1; 40 if(A==0)ans--; 41 printf("%lld ",ans); 42 } 43 } 44 int main() 45 { 46 T=read(); 47 while(T--)work(); 48 return 0; 49 }
【中国剩余定理】
〖相关资料〗
〖相关知识〗
假设质数 $m_{1} , m_{2} , ... , m_{n}$ 两两互质，则对于任意整数 $a_{1} , a_{2} , ... , a_{n}$ 同余方程组 $xequiv a_{i} mod m_{i}$ 在模 $M=m_{1}cdot m_{2}cdot ... cdot m_{k}$ 下有唯一解。令 $M_{i}=frac{M}{m_{i}}$ ，$t_{i}=M_{i} ^{-1} pmod {m_{i}}$ ，则 $x=left ( sum _{i=1}^{n} a_{i} t_{i} M_{i} ight )mod M$ 。
〖模板代码〗
1 void exgcd(int a,int b,int& x,int& y) 2 { 3 if(!b){x=1;y=0;return;} 4 exgcd(b,a%b,x,y); 5 int tmp=x;x=y;y=tmp-a/b*y; 6 } 7 int CRT() 8 { 9 int M=1,ans=0,x,y,mi; 10 for(int i=1;i<=n;i++)M*=m[i]; 11 for(int i=1;i<=n;i++) 12 { 13 mi=M/m[i]; 14 exgcd(mi,m[i],x,y); 15 ans=(ans+mi*x*a[i])%M; 16 } 17 if(ans<0)ans+=M; 18 return ans; 19 }
〖相关题目〗
1.【poj1006】Biorhythms
题意：人自出生起就有体力，情感和智力三个生理周期，分别为23，28和33天。一个周期内有一天为峰值，通常这三个周期的峰值不会是同一天。现在给出三个日期，分别对应于体力，情感，智力出现峰值的日期。然后再给出一个起始日期，要求从这一天开始，算出最少再过多少天后三个峰值同时出现。
分析：中国剩余定理裸题
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 using namespace std; 5 int n,d,ans,tim,a[4],m[4]; 6 int read() 7 { 8 int x=0,f=1;char c=getchar(); 9 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 10 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 11 return x*f; 12 } 13 void exgcd(int a,int b,int& x,int& y) 14 { 15 if(!b){x=1;y=0;return;} 16 exgcd(b,a%b,x,y); 17 int tmp=x;x=y;y=tmp-a/b*y; 18 } 19 int CRT() 20 { 21 int M=1,ans=0,x,y,mi; 22 for(int i=1;i<=n;i++)M*=m[i]; 23 for(int i=1;i<=n;i++) 24 { 25 mi=M/m[i]; 26 exgcd(mi,m[i],x,y); 27 ans=(ans+mi*x*a[i])%M; 28 } 29 if(ans<0)ans+=M; 30 return ans; 31 } 32 int main() 33 { 34 n=3;m[1]=23;m[2]=28;m[3]=33; 35 while(1) 36 { 37 a[1]=read();a[2]=read();a[3]=read();d=read(); 38 if(a[1]==-1&&a[2]==-1&&a[3]==-1&&d==-1)break; 39 ans=CRT();while(ans<=d)ans+=21252; 40 printf("Case %d: the next triple peak occurs in %d days. ",++tim,ans-d); 41 } 42 return 0; 43 }
2.【bzoj1951】古代猪文
题意：求G^M mod P，M=∑ i|N C(N,i)，P=999911659。
分析：hzwerの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=4e4+5; 7 const int P=999911659; 8 const int t[4]={2,3,4679,35617}; 9 LL n,g,tmp,m[4],fac[4][N]; 10 LL qm(LL a,LL b,LL mod) 11 { 12 LL ans=1;a%=mod; 13 while(b) 14 { 15 if(b&1)ans=ans*a%mod; 16 a=a*a%mod;b>>=1; 17 } 18 return ans; 19 } 20 LL inv(LL a,LL mod){return qm(a,mod-2,mod);} 21 LL C(LL n,LL m,int x) 22 { 23 if(n<m)return 0; 24 return fac[x][n]*inv(fac[x][m]*fac[x][n-m],t[x]); 25 } 26 LL Lucas(LL n,LL m,LL mod,int x) 27 { 28 if(m==0)return 1; 29 return C(n%mod,m%mod,x)*Lucas(n/mod,m/mod,mod,x)%mod; 30 } 31 LL CRT() 32 { 33 LL M=P-1,ans=0,mi; 34 for(int i=0;i<4;i++) 35 { 36 mi=M/t[i]; 37 ans=(ans+mi*inv(mi,t[i])%M*m[i])%M; 38 } 39 if(ans<0)ans+=M; 40 return ans; 41 } 42 int main() 43 { 44 scanf("%lld%lld",&n,&g); 45 for(int i=0;i<4;i++) 46 { 47 fac[i][0]=1; 48 for(int j=1;j<=t[i];j++) 49 fac[i][j]=fac[i][j-1]*j%t[i]; 50 } 51 if(g==P){printf("0");return 0;}g%=P; 52 for(LL i=1;i*i<=n;i++) 53 if(n%i==0) 54 { 55 tmp=n/i; 56 for(int j=0;j<4;j++) 57 { 58 m[j]=(m[j]+Lucas(n,i,t[j],j))%t[j]; 59 if(tmp!=i)m[j]=(m[j]+Lucas(n,tmp,t[j],j))%t[j]; 60 } 61 } 62 printf("%lld",qm(g,CRT(),P)); 63 return 0; 64 }
3.【bzoj2142】礼物
题意：见原题
分析：Sengxianの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=1e5+5; 7 int mod,n,m,sum,cnt; 8 int w[10],mm[10],a[10]; 9 int fac[N]; 10 int read() 11 { 12 int x=0,f=1;char c=getchar(); 13 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 14 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 15 return x*f; 16 } 17 int power(int a,int b,int mod) 18 { 19 int ans=1; 20 while(b) 21 { 22 if(b&1)ans=1ll*ans*a%mod; 23 a=1ll*a*a%mod;b>>=1; 24 } 25 return ans; 26 } 27 int Fac(int n,int P,int p,LL& e) 28 { 29 if(n<p)return fac[n]; 30 e+=n/p; 31 return 1ll*power(fac[P-1],n/P,P)*fac[n%P]%P*Fac(n/p,P,p,e)%P; 32 } 33 void exgcd(int a,int b,int& x,int& y) 34 { 35 if(!b){x=1;y=0;return;} 36 exgcd(b,a%b,x,y); 37 int tmp=x;x=y;y=tmp-a/b*y; 38 } 39 int inv(int a,int p) 40 { 41 int x,y;exgcd(a,p,x,y); 42 return (x%p+p)%p; 43 } 44 int cal(int P,int p) 45 { 46 fac[0]=fac[1]=1; 47 for(int i=2;i<P;i++)fac[i]=1ll*fac[i-1]*(i%p?i:1)%P; 48 LL e1=0,e2=0,a=Fac(n,P,p,e1),b=1; 49 for(int i=1;i<=m;i++)b=b*Fac(w[i],P,p,e2)%P; 50 return a*inv(b,P)%P*power(p,e1-e2,P)%P; 51 } 52 LL crt() 53 { 54 LL M=1,ans=0; 55 for(int i=1;i<=cnt;i++)M*=mm[i]; 56 for(int i=1;i<=cnt;i++) 57 { 58 LL w=M/mm[i]; 59 ans=(ans+w*inv(w,mm[i])%M*a[i]%M)%M; 60 } 61 return ans; 62 } 63 int main() 64 { 65 mod=read();n=read();m=read(); 66 for(int i=1;i<=m;i++)w[i]=read(),sum+=w[i]; 67 if(sum>n){printf("Impossible");return 0;} 68 if(sum<n)w[++m]=n-sum; 69 for(int i=2;1ll*i*i<=mod;i++) 70 if(mod%i==0) 71 { 72 int p=1;while(mod%i==0)mod/=i,p*=i; 73 mm[++cnt]=p;a[cnt]=cal(p,i); 74 } 75 if(mod>1)mm[++cnt]=mod,a[cnt]=cal(mod,mod); 76 printf("%lld",crt()); 77 return 0; 78 }
4.【bzoj1129】[POI2008]Per
题意：见原题
分析：Sengxianの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=3e5+5; 7 int n,mod,T,cnt,tmp; 8 int s[N],t[N],m[N],a[N],c[N]; 9 LL fac[N],x[N],e[N]; 10 int read() 11 { 12 int x=0,f=1;char c=getchar(); 13 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 14 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 15 return x*f; 16 } 17 int lowbit(int x){return x&(-x);} 18 void add(int x,int v){while(x<=T)t[x]+=v,x+=lowbit(x);} 19 int query(int x){int ans=0;while(x)ans+=t[x],x-=lowbit(x);return ans;} 20 int power(int a,int b,int mod) 21 { 22 int ans=1; 23 while(b) 24 { 25 if(b&1)ans=1ll*ans*a%mod; 26 a=1ll*a*a%mod;b>>=1; 27 } 28 return ans; 29 } 30 void exgcd(int a,int b,int& x,int& y) 31 { 32 if(!b){x=1;y=0;return;} 33 exgcd(b,a%b,x,y); 34 int tmp=x;x=y;y=tmp-a/b*y; 35 } 36 int inv(int a,int p) 37 { 38 int x,y;exgcd(a,p,x,y); 39 return (x%p+p)%p; 40 } 41 int cal(int P,int p) 42 { 43 memset(t,0,sizeof(t)); 44 memset(e,0,sizeof(e)); 45 memset(c,0,sizeof(c)); 46 x[0]=1; 47 for(int i=1;i<=n;i++) 48 { 49 e[i]=e[i-1];tmp=i; 50 while(tmp&&tmp%p==0)e[i]++,tmp/=p; 51 x[i]=x[i-1]*tmp%P; 52 } 53 for(int i=1;i<=n;i++)c[s[i]]++; 54 for(int i=1;i<=T;i++)add(i,c[i]); 55 LL e1=e[n-1],e2=0,a=x[n-1],b=1; 56 for(int i=1;i<=T;i++)e2+=e[c[i]],b=b*x[c[i]]%P; 57 LL w=query(s[1]-1); 58 while(w&&w%p==0)e1++,w/=p; 59 if(!w&&e1<e2)e1=e2; 60 LL ans=a*inv(b,P)%P*power(p,e1-e2,P)%P*w%P; 61 e2-=e[c[s[1]]];b=b*inv(x[c[s[1]]],P)%P; 62 c[s[1]]--;add(s[1],-1); 63 e2+=e[c[s[1]]];b=b*x[c[s[1]]]%P; 64 for(int i=2;i<=n;i++) 65 { 66 e1=e[n-i];a=x[n-i]; 67 LL w=query(s[i]-1); 68 while(w&&w%p==0)e1++,w/=p; 69 if(!w&&e1<e2)e1=e2; 70 w=w*a%P*inv(b,P)%P*power(p,e1-e2,P)%P; 71 ans=(ans+w)%P; 72 e2-=e[c[s[i]]];b=b*inv(x[c[s[i]]],P)%P; 73 c[s[i]]--;add(s[i],-1); 74 e2+=e[c[s[i]]];b=b*x[c[s[i]]]%P; 75 } 76 return ans; 77 } 78 int crt() 79 { 80 int M=1,ans=0; 81 for(int i=1;i<=cnt;i++)M*=m[i]; 82 for(int i=1;i<=cnt;i++) 83 { 84 int w=M/m[i]; 85 ans=(ans+1ll*w*inv(w,m[i])%M*a[i]%M)%M; 86 } 87 return (ans+1)%M; 88 } 89 int main() 90 { 91 n=read();mod=read(); 92 for(int i=1;i<=n;i++)s[i]=t[i]=read(); 93 sort(t+1,t+n+1); 94 T=unique(t+1,t+n+1)-t-1; 95 for(int i=1;i<=n;i++) 96 s[i]=lower_bound(t+1,t+T+1,s[i])-t; 97 for(int i=2;1ll*i*i<=mod;i++) 98 if(mod%i==0) 99 { 100 int p=1;while(mod%i==0)mod/=i,p*=i; 101 m[++cnt]=p;a[cnt]=cal(p,i); 102 } 103 if(mod>1)m[++cnt]=mod,a[cnt]=cal(mod,mod); 104 printf("%d ",crt()); 105 return 0; 106 }
【Lucas定理】
〖相关资料〗
《Lucas定理与大组合数的取模的求法总结》 《关于LUCAS二项式系数同余定理的一些推广》
〖相关知识〗
$Lucas(n,m,p)=inom{n mod p}{m mod p} cdot Lucas(frac{n}{p},frac{m}{p},p)$
〖模板代码〗
[Lucas定理]
1 LL qm(LL a,LL b) 2 { 3 LL ans=1; 4 while(b) 5 { 6 if(b&1)ans=ans*a%mod; 7 a=a*a%mod;b>>=1; 8 } 9 return ans; 10 } 11 LL inv(LL a){return qm(a,mod-2);} 12 LL C(LL n,LL m) 13 { 14 if(m>n)return 0; 15 LL a=1,b=1; 16 while(m) 17 { 18 a=a*n%mod;n--; 19 b=b*m%mod;m--; 20 } 21 return a*inv(b)%mod; 22 } 23 LL Lucas(LL n,LL m) 24 { 25 if(m==0)return 1; 26 return C(n%mod,m%mod)*Lucas(n/mod,m/mod)%mod; 27 }
[扩展Lucas定理]
1 LL qm(LL a,LL b,LL mod) 2 { 3 LL ans=1; 4 while(b) 5 { 6 if(b&1)ans=ans*a%mod; 7 a=a*a%mod;b>>=1; 8 } 9 return ans; 10 } 11 LL exgcd(LL a,LL b,LL& x,LL& y) 12 { 13 if(!b){x=1;y=0;return a;} 14 LL ans=exgcd(b,a%b,x,y); 15 LL tmp=x;x=y;y=tmp-a/b*y; 16 return ans; 17 } 18 LL inv(LL t,LL mod) 19 { 20 if(!t)return 0; 21 LL a=t,b=mod,x,y,g=exgcd(a,b,x,y); 22 if(g!=1)return 0; 23 x=(x%b+b)%b;if(!x)x=b; 24 return x; 25 } 26 LL mul(LL n,LL pi,LL pk) 27 { 28 if(!n)return 1; 29 LL ans=1; 30 if(n/pk) 31 { 32 for(LL i=2;i<=pk;i++) 33 if(i%pi)ans=ans*i%pk; 34 ans=qm(ans,n/pk,pk); 35 } 36 for(LL i=2;i<=n%pk;i++) 37 if(i%pi)ans=ans*i%pk; 38 return ans*mul(n/pi,pi,pk)%pk; 39 } 40 LL C(LL n,LL m,LL mod,LL pi,LL pk) 41 { 42 if(m>n)return 0; 43 LL a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk),k=0,ans; 44 for(LL i=n;i;i/=pi)k+=i/pi; 45 for(LL i=m;i;i/=pi)k-=i/pi; 46 for(LL i=n-m;i;i/=pi)k-=i/pi; 47 ans=a*inv(b,pk)%pk*inv(c,pk)%pk*qm(pi,k,pk)%pk; 48 return ans*(mod/pk)%mod*inv(mod/pk,pk)%mod; 49 } 50 LL F(LL n,LL m,LL mod) 51 { 52 LL ans=0,tmp=mod; 53 for(int i=2;i*i<=mod;i++) 54 if(tmp%i==0) 55 { 56 LL now=1; 57 while(tmp%i==0)now*=i,tmp/=i; 58 ans=(ans+C(n,m,mod,i,now))%mod; 59 } 60 if(tmp>1)ans=(ans+C(n,m,mod,tmp,tmp))%mod; 61 return ans; 62 }
〖相关题目〗
1.【hdu3037】Saving Beans
题意：求n个数的和不超过m的方案数。
分析：tju_virusの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 LL T,n,m,mod; 7 LL read() 8 { 9 LL x=0,f=1;char c=getchar(); 10 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 11 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 12 return x*f; 13 } 14 LL qm(LL a,LL b) 15 { 16 LL ans=1; 17 while(b) 18 { 19 if(b&1)ans=ans*a%mod; 20 a=a*a%mod;b>>=1; 21 } 22 return ans; 23 } 24 LL inv(LL a){return qm(a,mod-2);} 25 LL C(LL n,LL m) 26 { 27 if(m>n)return 0; 28 LL a=1,b=1; 29 while(m) 30 { 31 a=a*n%mod;n--; 32 b=b*m%mod;m--; 33 } 34 return a*inv(b)%mod; 35 } 36 LL Lucas(LL n,LL m) 37 { 38 if(m==0)return 1; 39 return C(n%mod,m%mod)*Lucas(n/mod,m/mod)%mod; 40 } 41 int main() 42 { 43 T=read(); 44 while(T--) 45 { 46 n=read();m=read();mod=read(); 47 printf("%lld ",Lucas(n+m,m)); 48 } 49 return 0; 50 }
2.【2015 ICL, Finals, Div. 1】J. Ceizenpok’s formula
题意：求C（n,k）%m，m不一定为质数。
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 LL n,m,mod,ans; 7 LL qm(LL a,LL b,LL mod) 8 { 9 LL ans=1; 10 while(b) 11 { 12 if(b&1)ans=ans*a%mod; 13 a=a*a%mod;b>>=1; 14 } 15 return ans; 16 } 17 LL exgcd(LL a,LL b,LL& x,LL& y) 18 { 19 if(!b){x=1;y=0;return a;} 20 LL ans=exgcd(b,a%b,x,y); 21 LL tmp=x;x=y;y=tmp-a/b*y; 22 return ans; 23 } 24 LL inv(LL t,LL mod) 25 { 26 if(!t)return 0; 27 LL a=t,b=mod,x,y,g=exgcd(a,b,x,y); 28 if(g!=1)return 0; 29 x=(x%b+b)%b;if(!x)x=b; 30 return x; 31 } 32 LL mul(LL n,LL pi,LL pk) 33 { 34 if(!n)return 1; 35 LL ans=1; 36 if(n/pk) 37 { 38 for(LL i=2;i<=pk;i++) 39 if(i%pi)ans=ans*i%pk; 40 ans=qm(ans,n/pk,pk); 41 } 42 for(LL i=2;i<=n%pk;i++) 43 if(i%pi)ans=ans*i%pk; 44 return ans*mul(n/pi,pi,pk)%pk; 45 } 46 LL C(LL n,LL m,LL mod,LL pi,LL pk) 47 { 48 if(m>n)return 0; 49 LL a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk),k=0,ans; 50 for(LL i=n;i;i/=pi)k+=i/pi; 51 for(LL i=m;i;i/=pi)k-=i/pi; 52 for(LL i=n-m;i;i/=pi)k-=i/pi; 53 ans=a*inv(b,pk)%pk*inv(c,pk)%pk*qm(pi,k,pk)%pk; 54 return ans*(mod/pk)%mod*inv(mod/pk,pk)%mod; 55 } 56 int main() 57 { 58 scanf("%I64d%I64d%I64d",&n,&m,&mod); 59 for(LL tmp=mod,i=2;i<=mod;i++) 60 if(tmp%i==0) 61 { 62 LL now=1; 63 while(tmp%i==0)now*=i,tmp/=i; 64 ans=(ans+C(n,m,mod,i,now))%mod; 65 } 66 printf("%I64d",ans); 67 return 0; 68 }
【BSGS算法】
〖相关资料〗
《Lucas定理&扩展Lucas定理&BSGS算法&扩展BSGS算法》
〖模板代码〗
[BSGS算法]
1 LL bsgs(LL a,LL b) 2 { 3 a%=mod;b%=mod;m.clear(); 4 if(!a&&!b)return 1; 5 if(!a)return -1; 6 int n=ceil(sqrt(mod));LL e=1; 7 m[1]=n+1; 8 for(int i=1;i<n;i++) 9 { 10 e=e*a%mod; 11 if(!m[e])m[e]=i; 12 } 13 e=power(e*a%mod,mod-2); 14 for(int i=0;i<n;i++) 15 { 16 LL now=m[b]; 17 if(now) 18 { 19 if(now==n+1)now=0; 20 return i*n+now; 21 } 22 b=b*e%mod; 23 } 24 return -1; 25 }
[扩展BSGS算法]
1 struct hashtable{LL to,next,w;}e[N]; 2 void add(LL v,LL id) 3 { 4 int hash=v%N; 5 e[++cnt]=(hashtable){v,first[hash],id}; 6 first[hash]=cnt; 7 } 8 LL find(LL v) 9 { 10 int hash=v%N; 11 for(int i=first[hash];i;i=e[i].next) 12 if(e[i].to==v)return e[i].w; 13 return -1; 14 } 15 LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);} 16 LL power(LL a,LL b,LL mod) 17 { 18 LL ans=1; 19 while(b) 20 { 21 if(b&1)ans=ans*a%mod; 22 a=a*a%mod;b>>=1; 23 } 24 return ans; 25 } 26 LL exbsgs() 27 { 28 a%=mod;b%=mod; 29 if(b==1)return 0; 30 LL g=1,d=1,num=0; 31 while((g=gcd(a,mod))!=1) 32 { 33 if(b%g)return -1; 34 num++;b/=g;mod/=g;d=d*(a/g)%mod; 35 if(b==d)return num; 36 } 37 LL m=ceil(sqrt(mod)),tmp=b,j; 38 cnt=0;memset(first,0,sizeof(first)); 39 for(int i=0;i<m;i++)add(tmp,i),tmp=tmp*a%mod; 40 tmp=power(a,m,mod);d=1; 41 for(int i=1;i<=m;i++) 42 { 43 d=d*tmp%mod; 44 if(~(j=find(d)))return i*m-j+num; 45 } 46 return -1; 47 }
〖相关题目〗
1.【bzoj2242】[SDOI2011]计算器
题意：见原题
分析：hzwerの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cstring> 5 #include<cmath> 6 #include<map> 7 #define LL long long 8 using namespace std; 9 LL T,type,a,b,mod,ans; 10 map<LL,LL> m; 11 LL read() 12 { 13 LL x=0,f=1;char c=getchar(); 14 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 15 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 16 return x*f; 17 } 18 LL power(LL a,LL b) 19 { 20 LL ans=1; 21 while(b) 22 { 23 if(b&1)ans=ans*a%mod; 24 a=a*a%mod;b>>=1; 25 } 26 return ans; 27 } 28 LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);} 29 void exgcd(LL a,LL b,LL& x,LL& y) 30 { 31 if(!b){x=1;y=0;return;} 32 exgcd(b,a%b,x,y); 33 LL t=x;x=y;y=t-a/b*y; 34 } 35 LL exg(LL a,LL b) 36 { 37 LL g=gcd(a,mod); 38 if(b%g)return -1; 39 a/=g;b/=g;mod/=g; 40 LL x,y;exgcd(a,mod,x,y); 41 x=x*b%mod;while(x<0)x+=mod; 42 return x; 43 } 44 LL bsgs(LL a,LL b) 45 { 46 a%=mod;b%=mod;m.clear(); 47 if(!a&&!b)return 1; 48 if(!a)return -1; 49 int n=ceil(sqrt(mod));LL e=1; 50 m[1]=n+1; 51 for(int i=1;i<n;i++) 52 { 53 e=e*a%mod; 54 if(!m[e])m[e]=i; 55 } 56 e=power(e*a%mod,mod-2); 57 for(int i=0;i<n;i++) 58 { 59 LL now=m[b]; 60 if(now) 61 { 62 if(now==n+1)now=0; 63 return i*n+now; 64 } 65 b=b*e%mod; 66 } 67 return -1; 68 } 69 int main() 70 { 71 T=read();type=read(); 72 while(T--) 73 { 74 a=read();b=read();mod=read(); 75 if(type==1)ans=power(a,b); 76 else if(type==2)ans=exg(a,b); 77 else ans=bsgs(a,b); 78 if(ans==-1)printf("Orz, I cannot find x! "); 79 else printf("%lld ",ans); 80 } 81 return 0; 82 }
2.【bzoj1467】Pku3243 clever Y
题意：X^Y mod Z=K。给出X、Z、K，计算Y。
分析：BeiYu_oiの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=99991; 8 LL a,b,mod,cnt,ans,first[N]; 9 LL read() 10 { 11 LL x=0,f=1;char c=getchar(); 12 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 13 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 14 return x*f; 15 } 16 struct hashtable{LL to,next,w;}e[N]; 17 void add(LL v,LL id) 18 { 19 int hash=v%N; 20 e[++cnt]=(hashtable){v,first[hash],id}; 21 first[hash]=cnt; 22 } 23 LL find(LL v) 24 { 25 int hash=v%N; 26 for(int i=first[hash];i;i=e[i].next) 27 if(e[i].to==v)return e[i].w; 28 return -1; 29 } 30 LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);} 31 LL power(LL a,LL b,LL mod) 32 { 33 LL ans=1; 34 while(b) 35 { 36 if(b&1)ans=ans*a%mod; 37 a=a*a%mod;b>>=1; 38 } 39 return ans; 40 } 41 LL exbsgs() 42 { 43 a%=mod;b%=mod; 44 if(b==1)return 0; 45 LL g=1,d=1,num=0; 46 while((g=gcd(a,mod))!=1) 47 { 48 if(b%g)return -1; 49 num++;b/=g;mod/=g;d=d*(a/g)%mod; 50 if(b==d)return num; 51 } 52 LL m=ceil(sqrt(mod)),tmp=b,j; 53 cnt=0;memset(first,0,sizeof(first)); 54 for(int i=0;i<m;i++)add(tmp,i),tmp=tmp*a%mod; 55 tmp=power(a,m,mod);d=1; 56 for(int i=1;i<=m;i++) 57 { 58 d=d*tmp%mod; 59 if(~(j=find(d)))return i*m-j+num; 60 } 61 return -1; 62 } 63 int main() 64 { 65 a=read();mod=read();b=read(); 66 while(a||b||mod) 67 { 68 ans=exbsgs(); 69 if(~ans)printf("%lld ",ans); 70 else printf("No Solution "); 71 a=read();mod=read();b=read(); 72 } 73 return 0; 74 }
【欧拉函数】
〖相关资料〗
《小于等于N的所有整数与N关于gcd(i,N)的那些事》（1..n中与n互质的数的和为n*φ(n)/2的证明）
〖相关知识〗
$varphi (1)=1$ ， $varphi (x)=xcdot (1-frac{1}{pri_{1}})cdot ... cdot (1-frac{1}{pri_{n}})$ （其中 $pri_{i}$ 是 $x$ 的所有质因数）。
$n=sum _{d|n}varphi(d)$ ， $sum _{i=1}^{n} i[gcd(i,n)=1]=frac{ncdot varphi(n)+[n=1]}{2}$ ， $sum _{i=1}^{n} gcd(i,n)=sum _{d|n} varphi(d) cdot frac{n}{d}$
〖模板代码〗
1 void getphi() 2 { 3 phi[1]=1; 4 for(int i=2;i<=n;i++) 5 { 6 if(!isp[i]){pri[++tot]=i;phi[i]=i-1;} 7 for(int j=1;j<=tot;j++) 8 { 9 if(i*pri[j]>n)break; 10 isp[i*pri[j]]=true; 11 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 12 else phi[i*pri[j]]=phi[i]*(pri[j]-1); 13 } 14 } 15 }
〖相关题目〗
1.【bzoj2190】[SDOI2008]仪仗队
题意：仪仗队是由学生组成的N * N的方阵，为了保证队伍在行进中整齐划一，C君会跟在仪仗队的左后方，根据其视线所及的学生人数来判断队伍是否整齐。现在，C君希望你告诉他队伍整齐时能看到的学生人数。
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=4e4+5; 7 int n,tot,ans,pri[N],phi[N]; 8 bool isp[N]; 9 void getphi() 10 { 11 phi[1]=1; 12 for(int i=2;i<=n;i++) 13 { 14 if(!isp[i]){pri[++tot]=i;phi[i]=i-1;} 15 for(int j=1;j<=tot;j++) 16 { 17 if(i*pri[j]>n)break; 18 isp[i*pri[j]]=true; 19 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 20 else phi[i*pri[j]]=phi[i]*(pri[j]-1); 21 } 22 } 23 } 24 int main() 25 { 26 scanf("%d",&n);getphi(); 27 for(int i=1;i<n;i++)ans+=phi[i]; 28 printf("%d",ans*2+1); 29 return 0; 30 }
2.【bzoj3643】Phi的反函数
题意：给定 x，求最小的n使得phi(n)=x。
分析：lych_cysの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=5e4+5; 8 int n,m,cnt,pri[N]; 9 LL ans=1ll*2147483647+1; 10 bool f[N]; 11 bool prm(int x) 12 { 13 int y=sqrt(x); 14 for(int i=1;pri[i]<=y;i++) 15 if(x%pri[i]==0)return false; 16 return true; 17 } 18 void dfs(int k,int now,LL sum) 19 { 20 if(sum>=ans)return; 21 if(now==1){ans=min(ans,sum);return;} 22 if(now>m&&prm(now+1))ans=min(ans,sum*(now+1)); 23 for(int i=k+1;pri[i]-1<=m;i++) 24 { 25 if(pri[i]-1>now)break; 26 if(now%(pri[i]-1))continue; 27 LL x=now/(pri[i]-1),y=sum*pri[i];dfs(i,x,y); 28 while(x%pri[i]==0) 29 { 30 x/=pri[i];y*=pri[i]; 31 dfs(i,x,y); 32 } 33 } 34 } 35 int main() 36 { 37 scanf("%d",&n);m=sqrt(n); 38 for(int i=2;i<=50000;i++) 39 { 40 if(!f[i])pri[++cnt]=i; 41 for(int j=1;j<=cnt;j++) 42 { 43 if(i*pri[j]>50000)break; 44 f[i*pri[j]]=true; 45 if(i%pri[j]==0)break; 46 } 47 } 48 dfs(0,n,1); 49 if(ans<=2147483647)printf("%lld",ans); 50 else printf("-1"); 51 return 0; 52 }
3.【bzoj2749】[HAOI2012]外星人
题意：给定一个数，用Πp[i]^q[i](p<=10^5,q<=10^9)的形式表示，问最少需要对这个数字x进行几次x=Φ(x)操作使得x=1。
分析：BearChildの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=1e5; 7 int T,n,cnt,p,q,odd,pri[N+5],f[N+5]; 8 LL ans; 9 int read() 10 { 11 int x=0,f=1;char c=getchar(); 12 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 13 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 14 return x*f; 15 } 16 void init() 17 { 18 f[1]=1; 19 for(int i=2;i<=N;i++) 20 { 21 if(!f[i]){pri[++cnt]=i;f[i]=f[i-1];} 22 for(int j=1;j<=cnt;j++) 23 { 24 if(i*pri[j]>N)break; 25 f[i*pri[j]]=f[i]+f[pri[j]]; 26 if(i%pri[j]==0)break; 27 } 28 } 29 } 30 int main() 31 { 32 init();T=read(); 33 while(T--) 34 { 35 n=read();odd=1;ans=0; 36 for(int i=1;i<=n;i++) 37 { 38 p=read();q=read(); 39 ans+=1ll*f[p]*q; 40 if(!(p&1))odd=0; 41 } 42 printf("%lld ",ans+odd); 43 } 44 return 0; 45 }
4.【bzoj1408】[Noi2002]Robot
题意：求m的所有约数中，满足可以分解成(奇数个不同素数/偶数个不同素数/其他)的所有的phi之和。
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int mod=1e4; 7 int n,p,e,ans1,ans2,t1,t2,m=1; 8 int read() 9 { 10 int x=0,f=1;char c=getchar(); 11 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 12 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 13 return x*f; 14 } 15 int qm(int a,int b) 16 { 17 int ans=1; 18 while(b) 19 { 20 if(b&1)ans=ans*a%mod; 21 a=a*a%mod;b>>=1; 22 } 23 return ans; 24 } 25 int main() 26 { 27 n=read(); 28 for(int i=1;i<=n;i++) 29 { 30 p=read();e=read(); 31 m=m*qm(p,e)%mod; 32 if(p==2)continue; 33 t1=(ans1+ans2*(p-1)%mod)%mod; 34 t2=(ans2+(ans1+1)*(p-1)%mod)%mod; 35 ans1=t1;ans2=t2; 36 } 37 printf("%d %d %d",ans1,ans2,((m-1-ans1-ans2)%mod+mod)%mod); 38 return 0; 39 }
5.【codeforces870F】Paths
题意：给定数字n，建立一个无向图。对于所有1~n之间的数字，当数字gcd(u,v)≠1时将u、v连一条边，边权为1。d(u, v)表示u到v的最短路，求所有d(u, v)的和，其中1 ≤ u < v ≤ n。
分析：Zsnuoの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #define LL long long 5 using namespace std; 6 const int N=1e7+5; 7 int n,m,tot,now,pri[N],p[N],phi[N],num[N],sum[N]; 8 LL one,two,three; 9 int main() 10 { 11 scanf("%d",&n); 12 phi[1]=1; 13 for(int i=2;i<=n;i++) 14 { 15 if(!p[i]){p[i]=pri[++tot]=i;phi[i]=i-1;} 16 for(int j=1;j<=tot;j++) 17 { 18 if(i*pri[j]>n)break; 19 p[i*pri[j]]=pri[j]; 20 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 21 else phi[i*pri[j]]=phi[i]*(pri[j]-1); 22 } 23 } 24 for(int i=2;i<=n;i++)one+=i-1-phi[i]; 25 for(int i=2;i<=n;i++)num[p[i]]++; 26 for(int i=2;i<=n;i++)sum[i]=sum[i-1]+num[i]; 27 for(int i=2;i<=n;i++)two+=1ll*num[i]*sum[n/i]; 28 for(int i=2;i<=n;i++)if(1ll*p[i]*p[i]<=n)two--; 29 two=two/2-one;m=n-1; 30 for(int i=tot;i>=1;i--) 31 if(pri[i]*2>n)m--; 32 else break; 33 three=1ll*m*(m-1)/2-one-two; 34 printf("%I64d ",one+two*2+three*3); 35 return 0; 36 }
【素性测试】
〖相关资料〗
〖模板代码〗
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #include<cstdlib> 6 #define LL long long 7 using namespace std; 8 const int N=105; 9 const LL pri[9]={2,3,5,7,11,13,17,19,23}; 10 LL n,cnt,a[N]; 11 LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);} 12 LL mul(LL x,LL y,LL mo) 13 { 14 LL tmp=(x*y-(LL)((double)x*y/mo+0.1)*mo)%mo; 15 if(tmp<0)tmp+=mo;return tmp; 16 } 17 LL ksm(LL x,LL y,LL mo) 18 { 19 LL ans=1; 20 while(y) 21 { 22 if(y&1)ans=mul(ans,x,mo); 23 x=mul(x,x,mo);y>>=1; 24 } 25 return ans; 26 } 27 bool MR(LL n) 28 { 29 if(n==2)return true; 30 if(n%2==0)return false; 31 int lg=0;LL w=n-1; 32 while(w%2==0)lg++,w>>=1; 33 for(int i=0;i<9;i++) 34 { 35 if(n==pri[i])return true; 36 LL x=ksm(pri[i],w,n); 37 if(x==1||x==n-1)continue; 38 for(int j=1;j<=lg;j++) 39 { 40 LL y=mul(x,x,n); 41 if(x!=1&&x!=n-1&&y==1)return false; 42 x=y; 43 } 44 if(x!=1)return false; 45 } 46 return true; 47 } 48 LL rho(LL n) 49 { 50 LL c=rand()*rand()%(n-1)+1,x1=rand()*rand()%n,x2=x1,k=2,p=1; 51 for(int i=1;p==1;i++) 52 { 53 x1=(mul(x1,x1,n)+c)%n; 54 if(x1==x2)return 1; 55 p=gcd(n,abs(x1-x2)); 56 if(i==k)k<<=1,x2=x1; 57 } 58 return p; 59 } 60 void div(LL n) 61 { 62 if(n==1)return; 63 if(MR(n)){a[++cnt]=n;return;} 64 LL p=1; 65 while(p==1)p=rho(n); 66 div(p);div(n/p); 67 } 68 int main() 69 { 70 scanf("%lld",&n);div(n); 71 sort(a+1,a+cnt+1); 72 cnt=unique(a+1,a+cnt+1)-a-1; 73 return 0; 74 }
〖相关题目〗
1.【bzoj4802】欧拉函数
题意：已知N，求phi(N)。N<=1018。
分析：Galaxiesの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #include<cstdlib> 6 #define LL long long 7 using namespace std; 8 const int N=105; 9 const LL pri[9]={2,3,5,7,11,13,17,19,23}; 10 LL n,cnt,a[N]; 11 LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);} 12 LL mul(LL x,LL y,LL mo) 13 { 14 LL tmp=(x*y-(LL)((double)x*y/mo+0.1)*mo)%mo; 15 if(tmp<0)tmp+=mo;return tmp; 16 } 17 LL ksm(LL x,LL y,LL mo) 18 { 19 LL ans=1; 20 while(y) 21 { 22 if(y&1)ans=mul(ans,x,mo); 23 x=mul(x,x,mo);y>>=1; 24 } 25 return ans; 26 } 27 bool MR(LL n) 28 { 29 if(n==2)return true; 30 if(n%2==0)return false; 31 int lg=0;LL w=n-1; 32 while(w%2==0)lg++,w>>=1; 33 for(int i=0;i<9;i++) 34 { 35 if(n==pri[i])return true; 36 LL x=ksm(pri[i],w,n); 37 if(x==1||x==n-1)continue; 38 for(int j=1;j<=lg;j++) 39 { 40 LL y=mul(x,x,n); 41 if(x!=1&&x!=n-1&&y==1)return false; 42 x=y; 43 } 44 if(x!=1)return false; 45 } 46 return true; 47 } 48 LL rho(LL n) 49 { 50 LL c=rand()*rand()%(n-1)+1,x1=rand()*rand()%n,x2=x1,k=2,p=1; 51 for(int i=1;p==1;i++) 52 { 53 x1=(mul(x1,x1,n)+c)%n; 54 if(x1==x2)return 1; 55 p=gcd(n,abs(x1-x2)); 56 if(i==k)k<<=1,x2=x1; 57 } 58 return p; 59 } 60 void div(LL n) 61 { 62 if(n==1)return; 63 if(MR(n)){a[++cnt]=n;return;} 64 LL p=1; 65 while(p==1)p=rho(n); 66 div(p);div(n/p); 67 } 68 int main() 69 { 70 srand(20010311); 71 scanf("%lld",&n);div(n); 72 sort(a+1,a+cnt+1); 73 cnt=unique(a+1,a+cnt+1)-a-1; 74 for(int i=1;i<=cnt;i++)n=n/a[i]*(a[i]-1); 75 printf("%lld",n); 76 return 0; 77 }
【莫比乌斯反演】
〖相关资料〗
〖相关知识〗
$F(n)=sum _{d|n} f(d)Rightarrow f(n)=sum _{d|n} mu (d)F(frac{n}{d})$
$F(n)=sum _{n|d} f(d)Rightarrow f(n)=sum _{n|d} mu (frac{d}{n})F(d)$
$sum _{d|n} mu(d)=left{egin{matrix} 1(n=1)\ 0(n>1) end{matrix} ight.$ ， $sum _{d|n} frac{mu(d)}{d}=frac{varphi(n)}{n}$
〖模板代码〗
1 void pre() 2 { 3 miu[1]=1; 4 for(int i=2;i<=N;i++) 5 { 6 if(!np[i]){miu[i]=-1;pri[++cnt]=i;} 7 for(int j=1;i*pri[j]<=N;j++) 8 { 9 np[i*pri[j]]=true; 10 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 11 miu[i*pri[j]]=-miu[i]; 12 } 13 } 14 }
〖相关题目〗
1.【bzoj2440】[中山市选2011]完全平方数
题意：求第k个无平方因子数。无平方因子数，即分解之后所有质因数的次数都为1的数。
分析：对莫比乌斯函数的应用
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=50000; 7 int T,n,cnt,pri[N+5],miu[N+5]; 8 bool np[N+5]; 9 void pre() 10 { 11 miu[1]=1; 12 for(int i=2;i<=N;i++) 13 { 14 if(!np[i]){miu[i]=-1;pri[++cnt]=i;} 15 for(int j=1;i*pri[j]<=N;j++) 16 { 17 np[i*pri[j]]=true; 18 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 19 miu[i*pri[j]]=-miu[i]; 20 } 21 } 22 } 23 bool check(LL x) 24 { 25 LL sum=0; 26 for(LL i=1;i*i<=x;i++)sum+=x/(i*i)*miu[i]; 27 return sum>=n; 28 } 29 int main() 30 { 31 pre(); 32 scanf("%d",&T); 33 while(T--) 34 { 35 scanf("%d",&n); 36 LL l=n,r=1644934081,ans=0,mid; 37 while(l<=r) 38 { 39 mid=(l+r)>>1; 40 if(check(mid))ans=mid,r=mid-1; 41 else l=mid+1; 42 } 43 printf("%lld ",ans); 44 } 45 return 0; 46 }
2.【bzoj2301】[HAOI2011]Problem b
题意：对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k。
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=50005; 7 int n,a,b,c,d,k,cnt; 8 int pri[N],miu[N],sum[N]; 9 bool np[N]; 10 int read() 11 { 12 int x=0,f=1;char c=getchar(); 13 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 14 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 15 return x*f; 16 } 17 void pre() 18 { 19 miu[1]=1; 20 for(int i=2;i<=50000;i++) 21 { 22 if(!np[i]){pri[++cnt]=i;miu[i]=-1;} 23 for(int j=1;i*pri[j]<=50000;j++) 24 { 25 np[i*pri[j]]=true; 26 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 27 else miu[i*pri[j]]=-miu[i]; 28 } 29 } 30 for(int i=1;i<=50000;i++)sum[i]=sum[i-1]+miu[i]; 31 } 32 int calc(int n,int m) 33 { 34 if(n>m)swap(n,m); 35 int ans=0,pos; 36 for(int i=1;i<=n;i=pos+1) 37 { 38 pos=min(n/(n/i),m/(m/i)); 39 ans+=(sum[pos]-sum[i-1])*(n/i)*(m/i); 40 } 41 return ans; 42 } 43 int main() 44 { 45 pre();n=read(); 46 while(n--) 47 { 48 a=read()-1;b=read();c=read()-1;d=read();k=read(); 49 a/=k;b/=k;c/=k;d/=k; 50 printf("%d ",calc(a,c)+calc(b,d)-calc(a,d)-calc(b,c)); 51 } 52 return 0; 53 }
3.【bzoj2820】YY的GCD
题意：求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对。
分析：DQSSSの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=1e7+5; 7 const int mx=1e7; 8 int T,n,m,pos,cnt; 9 int pri[N],miu[N]; 10 LL ans,f[N]; 11 bool np[N]; 12 int read() 13 { 14 int x=0,f=1;char c=getchar(); 15 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 16 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 17 return x*f; 18 } 19 void pre() 20 { 21 miu[1]=1; 22 for(int i=2;i<=mx;i++) 23 { 24 if(!np[i]){pri[++cnt]=i;miu[i]=-1;} 25 for(int j=1;i*pri[j]<=mx;j++) 26 { 27 np[i*pri[j]]=true; 28 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 29 else miu[i*pri[j]]=-miu[i]; 30 } 31 } 32 for(int i=1;i<=cnt;i++) 33 { 34 int p=pri[i]; 35 for(int j=1;j*p<=mx;j++)f[j*p]+=miu[j]; 36 } 37 for(int i=1;i<=mx;i++)f[i]+=f[i-1]; 38 } 39 int main() 40 { 41 pre();T=read(); 42 while(T--) 43 { 44 n=read();m=read(); 45 if(n>m)swap(n,m);ans=0; 46 for(int i=1;i<=n;i=pos+1) 47 { 48 pos=min(n/(n/i),m/(m/i)); 49 ans+=(f[pos]-f[i-1])*(n/i)*(m/i); 50 } 51 printf("%lld ",ans); 52 } 53 return 0; 54 }
4.【bzoj1101】[POI2007]Zap
题意：对于给定的整数a,b和d，有多少正整数对x,y，满足x<=a，y<=b，并且gcd(x,y)=d。
分析：hzwerの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 const int N=5e4+5; 7 const int mx=5e4; 8 int T,n,m,d,pos,cnt,ans; 9 int pri[N],miu[N],sum[N]; 10 bool np[N]; 11 int read() 12 { 13 int x=0,f=1;char c=getchar(); 14 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 15 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 16 return x*f; 17 } 18 void pre() 19 { 20 miu[1]=1; 21 for(int i=2;i<=mx;i++) 22 { 23 if(!np[i]){pri[++cnt]=i;miu[i]=-1;} 24 for(int j=1;i*pri[j]<=mx;j++) 25 { 26 np[i*pri[j]]=true; 27 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 28 else miu[i*pri[j]]=-miu[i]; 29 } 30 } 31 for(int i=1;i<=mx;i++)sum[i]=sum[i-1]+miu[i]; 32 } 33 int main() 34 { 35 pre();T=read(); 36 while(T--) 37 { 38 n=read();m=read();d=read(); 39 n/=d;m/=d; 40 if(n>m)swap(n,m);ans=0; 41 for(int i=1;i<=n;i=pos+1) 42 { 43 pos=min(n/(n/i),m/(m/i)); 44 ans+=(sum[pos]-sum[i-1])*(n/i)*(m/i); 45 } 46 printf("%d ",ans); 47 } 48 return 0; 49 }
5.【bzoj3529】[Sdoi2014]数表
题意：有一张n×m的数表，其第i行第j列（1<=i<=n，1<=j<=m）的数值为能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。
分析：PoPoQQQのPPT，见资料
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<iostream> 5 #define LL long long 6 using namespace std; 7 const int N=1e5+5; 8 const int M=2e4+5; 9 int T,cnt,mx,now,num; 10 int pri[N],miu[N],tr[N],ans[M]; 11 bool np[N]; 12 pair<int,int>f[N]; 13 struct node{int n,m,a,id;}q[M]; 14 int read() 15 { 16 int x=0,f=1;char c=getchar(); 17 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 18 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 19 return x*f; 20 } 21 bool operator < (node a,node b){return a.a<b.a;} 22 int lowbit(int x){return x&(-x);} 23 void add(int x,int v){while(x<=mx)tr[x]+=v,x+=lowbit(x);} 24 int query(int x){int ans=0;while(x)ans+=tr[x],x-=lowbit(x);return ans;} 25 void pre() 26 { 27 miu[1]=1; 28 for(int i=2;i<=mx;i++) 29 { 30 if(!np[i]){pri[++cnt]=i;miu[i]=-1;} 31 for(int j=1;i*pri[j]<=mx;j++) 32 { 33 np[i*pri[j]]=true; 34 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 35 else miu[i*pri[j]]=-miu[i]; 36 } 37 } 38 for(int i=1;i<=mx;i++) 39 for(int j=i;j<=mx;j+=i) 40 f[j].first+=i; 41 for(int i=1;i<=mx;i++)f[i].second=i; 42 } 43 void work(int x) 44 { 45 int n=q[x].n,m=q[x].m,id=q[x].id,pos; 46 for(int i=1;i<=n;i=pos+1) 47 { 48 pos=min(n/(n/i),m/(m/i)); 49 ans[id]+=(query(pos)-query(i-1))*(n/i)*(m/i); 50 } 51 } 52 int main() 53 { 54 T=read(); 55 for(int i=1;i<=T;i++) 56 { 57 q[i].n=read();q[i].m=read();q[i].a=read();q[i].id=i; 58 if(q[i].n>q[i].m)swap(q[i].n,q[i].m); 59 mx=max(q[i].n,mx); 60 } 61 pre();sort(q+1,q+T+1);sort(f+1,f+mx+1); 62 for(int i=1;i<=T;i++) 63 { 64 while(now+1<=mx&&f[now+1].first<=q[i].a) 65 { 66 now++;num=f[now].second; 67 for(int j=num;j<=mx;j+=num) 68 add(j,f[now].first*miu[j/num]); 69 } 70 work(i); 71 } 72 for(int i=1;i<=T;i++)printf("%d ",ans[i]&0x7fffffff); 73 return 0; 74 }
6.【bzoj2154】Crash的数字表格
题意：一张N*M的表格，第i行第j列的那个格子里写着数为LCM(i, j)，求表格里所有数的和mod20101009的值。
分析：PoPoQQQのPPT，见资料
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<iostream> 5 #define LL long long 6 using namespace std; 7 const int N=1e7+5; 8 const int mod=20101009; 9 LL n,m,cnt,pri[N],miu[N],s[N]; 10 bool np[N]; 11 void pre() 12 { 13 miu[1]=1; 14 for(int i=2;i<=n;i++) 15 { 16 if(!np[i]){pri[++cnt]=i;miu[i]=-1;} 17 for(int j=1;i*pri[j]<=n;j++) 18 { 19 np[i*pri[j]]=true; 20 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 21 else miu[i*pri[j]]=-miu[i]; 22 } 23 } 24 for(LL i=1;i<=n;i++)s[i]=(s[i-1]+(i*i*miu[i])%mod)%mod; 25 } 26 LL sum(LL x,LL y){return (x*(x+1)/2%mod)*(y*(y+1)/2%mod)%mod;} 27 LL f(LL n,LL m) 28 { 29 if(n>m)swap(n,m); 30 LL ans=0,pos; 31 for(LL i=1;i<=n;i=pos+1) 32 { 33 pos=min(n/(n/i),m/(m/i)); 34 ans=((ans+(s[pos]-s[i-1])*sum(n/i,m/i)%mod)%mod+mod)%mod; 35 } 36 return ans; 37 } 38 int main() 39 { 40 scanf("%lld%lld",&n,&m); 41 if(n>m)swap(n,m);pre(); 42 LL ans=0,pos; 43 for(LL i=1;i<=n;i=pos+1) 44 { 45 pos=min(n/(n/i),m/(m/i)); 46 ans=(ans+(i+pos)*(pos-i+1)/2%mod*f(n/i,m/i)%mod)%mod; 47 } 48 printf("%lld ",ans); 49 return 0; 50 }
7.【bzoj2693】jzptab
题意：一张N*M的表格，第i行第j列的那个格子里写着数为LCM(i, j)，求表格里所有数的和mod1e8+9的值。
分析：PoPoQQQのPPT，见资料
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<iostream> 5 #define LL long long 6 using namespace std; 7 const int N=1e7+5; 8 const int mx=1e7; 9 const int mod=1e8+9; 10 LL T,n,m,cnt,ans,pos; 11 LL pri[N],miu[N],h[N],s[N]; 12 bool np[N]; 13 void pre() 14 { 15 miu[1]=1;h[1]=1; 16 for(LL i=2;i<=mx;i++) 17 { 18 if(!np[i]){pri[++cnt]=i;miu[i]=-1;h[i]=(i-i*i)%mod;} 19 for(LL j=1;i*pri[j]<=mx;j++) 20 { 21 np[i*pri[j]]=true; 22 if(i%pri[j]==0) 23 { 24 miu[i*pri[j]]=0; 25 h[i*pri[j]]=h[i]*pri[j]%mod; 26 break; 27 } 28 miu[i*pri[j]]=-miu[i]; 29 h[i*pri[j]]=h[i]*h[pri[j]]%mod; 30 } 31 } 32 for(LL i=1;i<=mx;i++)s[i]=s[i-1]+h[i]; 33 } 34 LL sum(LL x,LL y){return (x*(x+1)/2%mod)*(y*(y+1)/2%mod)%mod;} 35 int main() 36 { 37 pre();scanf("%lld",&T); 38 while(T--) 39 { 40 scanf("%lld%lld",&n,&m); 41 if(n>m)swap(n,m);ans=0; 42 for(LL i=1;i<=n;i=pos+1) 43 { 44 pos=min(n/(n/i),m/(m/i)); 45 ans=((ans+sum(n/i,m/i)*(s[pos]-s[i-1])%mod)%mod+mod)%mod; 46 } 47 printf("%lld ",ans); 48 } 49 return 0; 50 }
8.【51nod1222】最小公倍数计数
题意：见原题
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<algorithm> 5 #define LL long long 6 using namespace std; 7 const int N=320000; 8 int tot,pri[N+5],miu[N+5]; 9 LL a,b; 10 bool np[N+5]; 11 LL calc(LL n) 12 { 13 if(n==0)return 0; 14 LL upk=sqrt(n),ans=0,tmp,up; 15 for(LL k=1;k<=upk;k++) 16 { 17 if(!miu[k])continue; 18 up=n/(k*k);tmp=0; 19 for(LL d=1;d*d*d<=up;d++) 20 { 21 for(LL i=d+1;i*i*d<=up;i++) 22 tmp+=(up/(i*d)-i)*6+3; 23 tmp+=(up/(d*d)-d)*3; 24 tmp++; 25 } 26 ans+=miu[k]*tmp; 27 } 28 return (ans+n)/2; 29 } 30 int main() 31 { 32 miu[1]=1; 33 for(int i=2;i<=N;i++) 34 { 35 if(!np[i]){miu[i]=-1;pri[++tot]=i;} 36 for(int j=1;i*pri[j]<=N;j++) 37 { 38 np[i*pri[j]]=true; 39 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 40 miu[i*pri[j]]=-miu[i]; 41 } 42 } 43 scanf("%lld%lld",&a,&b); 44 printf("%lld",calc(b)-calc(a-1)); 45 return 0; 46 }
【扩展欧拉定理】
〖相关知识〗
$a^{b}equiv egin{cases}
a^{b\%varphi(p)} ~~~~~~~~~~~~(gcd(a,p)=1)\
a^{b}~~~~~~~~~~~~~~~~~~~~ (gcd(a,p)
eq 1,b<varphi(p))\
a^{b\% varphi(p)+varphi(p)}~~~~(gcd(a,p)
eq 1,bgeq varphi(p))
end{cases} pmod{p}$
〖相关题目〗
1.【codeforces906D】Power Tower
题意：见原题
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<map> 5 #define LL long long 6 using namespace std; 7 const int N=1e5+5; 8 LL n,m,q,L,R,w[N]; 9 map<LL,LL> mp; 10 LL read() 11 { 12 LL x=0,f=1;char c=getchar(); 13 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 14 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 15 return x*f; 16 } 17 LL Mod(LL a,LL b){return a<b?a:a%b+b;} 18 LL phi(LL n) 19 { 20 if(mp.count(n))return mp[n]; 21 LL tmp=n,s=n; 22 for(LL i=2;i*i<=n;i++) 23 { 24 if(n%i==0)s=s/i*(i-1); 25 while(n%i==0)n/=i; 26 } 27 if(n>1)s=s/n*(n-1); 28 mp[tmp]=s;return s; 29 } 30 LL power(LL a,LL b,LL mod) 31 { 32 LL ans=1; 33 while(b) 34 { 35 if(b&1)ans=Mod(ans*a,mod); 36 a=Mod(a*a,mod);b>>=1; 37 } 38 return ans; 39 } 40 LL calc(LL L,LL R,LL mod) 41 { 42 if(L==R||mod==1)return Mod(w[L],mod); 43 return power(w[L],calc(L+1,R,phi(mod)),mod); 44 } 45 int main() 46 { 47 n=read();m=read(); 48 for(int i=1;i<=n;i++)w[i]=read(); 49 q=read(); 50 while(q--) 51 { 52 L=read();R=read(); 53 printf("%lld ",calc(L,R,m)%m); 54 } 55 return 0; 56 }
2.【bzoj3884】上帝与集合的正确用法
题意：求2∞mod P。
分析：无
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define LL long long 5 using namespace std; 6 int T,mod; 7 int read() 8 { 9 int x=0,f=1;char c=getchar(); 10 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 11 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 12 return x*f; 13 } 14 LL phi(LL n) 15 { 16 LL s=n; 17 for(LL i=2;i*i<=n;i++) 18 { 19 if(n%i==0)s=s/i*(i-1); 20 while(n%i==0)n/=i; 21 } 22 if(n>1)s=s/n*(n-1); 23 return s; 24 } 25 LL Mod(LL a,LL b){return a<b?a:a%b+b;} 26 LL power(LL a,LL b,LL mod) 27 { 28 LL ans=1; 29 while(b) 30 { 31 if(b&1)ans=Mod(ans*a,mod); 32 a=Mod(a*a,mod);b>>=1; 33 } 34 return ans; 35 } 36 LL calc(LL mod) 37 { 38 if(mod==1)return 1; 39 return power(2,calc(phi(mod)),mod); 40 } 41 int main() 42 { 43 T=read(); 44 while(T--) 45 { 46 mod=read(); 47 printf("%lld ",calc(mod)%mod); 48 } 49 return 0; 50 }
3.【bzoj4869】[Shoi2017]相逢是问候
题意：见原题
分析：llgycの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #define l(x) x<<1 5 #define r(x) x<<1|1 6 #define LL long long 7 using namespace std; 8 const int N=5e4+5; 9 int n,m,mod,c,k,op,L,R,ans; 10 int p[N],a[N],cnt[N*4],t[N*4]; 11 int read() 12 { 13 int x=0,f=1;char c=getchar(); 14 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 15 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 16 return x*f; 17 } 18 int phi(int n) 19 { 20 int s=n; 21 for(int i=2;i*i<=n;i++) 22 { 23 if(n%i==0)s=s/i*(i-1); 24 while(n%i==0)n/=i; 25 } 26 if(n>1)s=s/n*(n-1); 27 return s; 28 } 29 void up(int x) 30 { 31 t[x]=(t[l(x)]+t[r(x)])%mod; 32 cnt[x]=min(cnt[l(x)],cnt[r(x)]); 33 } 34 void build(int x,int l,int r) 35 { 36 if(l==r){t[x]=a[l];return;} 37 int mid=(l+r)>>1; 38 build(l(x),l,mid); 39 build(r(x),mid+1,r); 40 up(x); 41 } 42 int Mod(LL a,int b){return a<b?a:a%b+b;} 43 int power(int a,int b,int mod) 44 { 45 int ans=1; 46 while(b) 47 { 48 if(b&1)ans=Mod(1ll*ans*a,mod); 49 a=Mod(1ll*a*a,mod);b>>=1; 50 } 51 return ans; 52 } 53 int calc(int d,int x) 54 { 55 int ans=x;ans=Mod(ans,p[d]); 56 while(d){d--;ans=power(c,ans,p[d]);} 57 return ans%mod; 58 } 59 void modify(int x,int l,int r) 60 { 61 if(cnt[x]>=k)return; 62 if(l==r) 63 { 64 cnt[x]++; 65 t[x]=calc(cnt[x],a[l]); 66 return; 67 } 68 int mid=(l+r)>>1; 69 if(L<=mid)modify(l(x),l,mid); 70 if(R>mid)modify(r(x),mid+1,r); 71 up(x); 72 } 73 void query(int x,int l,int r) 74 { 75 if(L<=l&&r<=R){ans=(ans+t[x])%mod;return;} 76 int mid=(l+r)>>1; 77 if(L<=mid)query(l(x),l,mid); 78 if(R>mid)query(r(x),mid+1,r); 79 } 80 int main() 81 { 82 n=read();m=read();mod=read();c=read(); 83 for(int i=1;i<=n;i++)a[i]=read(); 84 p[0]=mod;while(p[k]!=1)k++,p[k]=phi(p[k-1]);p[++k]=1; 85 build(1,1,n); 86 while(m--) 87 { 88 op=read();L=read();R=read(); 89 if(!op)modify(1,1,n); 90 else ans=0,query(1,1,n),printf("%d ",ans); 91 } 92 return 0; 93 }
【高斯消元】
〖模板代码〗
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=105; 8 int n,r; 9 double a[N][N]; 10 bool flag; 11 int read() 12 { 13 int x=0,f=1;char c=getchar(); 14 while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();} 15 while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();} 16 return x*f; 17 } 18 void gauss() 19 { 20 for(int i=1;i<=n;i++) 21 { 22 r=i; 23 for(int j=i+1;j<=n;j++) 24 if(fabs(a[j][i])>fabs(a[r][i]))r=j; 25 if(fabs(a[r][i])<1e-8){flag=true;return;} 26 if(r!=i) 27 for(int j=1;j<=n+1;j++) 28 swap(a[r][j],a[i][j]); 29 for(int k=i+1;k<=n;k++) 30 for(int j=n+1;j>=i;j--)//!!! 31 a[k][j]-=a[k][i]/a[i][i]*a[i][j]; 32 } 33 for(int i=n;i>=1;i--) 34 { 35 for(int j=i+1;j<=n;j++) 36 a[i][n+1]-=a[j][n+1]*a[i][j]; 37 a[i][n+1]/=a[i][i]; 38 } 39 } 40 int main() 41 { 42 n=read(); 43 for(int i=1;i<=n;i++) 44 for(int j=1;j<=n+1;j++) 45 a[i][j]=read(); 46 gauss(); 47 if(flag)printf("No Solution"); 48 else 49 for(int i=1;i<=n;i++) 50 printf("%.2lf ",a[i][n+1]); 51 return 0; 52 }