数论模板
#include<iostream> #include<cstdio> #include<cstring> #include<cmath> #include<algorithm> #define MAXN 1000000 using namespace std; bool not_prime[MAXN]; int prime_number[MAXN]; int nu; int factor[MAXN];//顺手记录一下最小质因子 void prime()//线筛 { for (int i=2;i<=MAXN;i++) { if (!not_prime[i]) { prime_number[++nu]=i; factor[i]=i; } for (int j=1;j<=nu&&i*prime_number[j]<=MAXN;j++) { not_prime[i*prime_number[j]]=1; f[i*prime_number[j]]=prime_number[j]; if (i%prime_number[j]==0) break; } } } int mu[MAXN]={0,1}; void mu()//莫比乌斯函数 { for (int i=2;i<=MAXN;i++) { if (!not_prime[i]) mu[i]=-1; for (int j=1;j<=nu&&i*prime_number[j]<=MAXN;j++) { not_prime[i*prime_number[j]]=1; if (i%prime_number[j]==0) { mu[i*prime_number[j]]=0; break; } else { mu[i*prime_number[j]]=-mu[i]; } } } } int phi[MAXN]; void phi()//欧拉函数 { for (int i=2;i<=MAXN;i++) { if (!not_prime[i]) { phi[i]=i-1; } for (int j=1;j<=nu&&i*prime_number[j]<=MAXN;j++) { not_prime[i*prime_number[j]]=1; if (i%prime_number[j]==0) { phi[i*prime_number[j]]=phi[i]*prime_number[j]; break; } else { phi[i*prime_number[j]]=phi[i]*(prime_number[j]-1); } } } } int inv[MAXN]; //以下连续逆元 void inv1()//费马小定理 { for (int i=1;i<=MAXN;i++) inv[i]=pow(i,phi[p]-1)%p;//p根据题目要求不同来定 } void inv2()//基于线筛的求法 { for (int i=2;i<=MAXN;i++) { int a=p/i; int b=p%i; inv[i]=(p-a)*inv[b]%p; } } void inv3()//基于阶乘的求法 { unsigned long long a[MAXN]={1},rev[MAXN];//这个阶乘略大0-0不想写高精度了unsigned凑活看着吧 for (int i=1;i<=MAXN;i++) a[i]=a[i-1]*i; rev[MAXN]=pow(a[MAXN],phi[p]-1)%p; for (int i=MAXN-1;i>=0;i--) rev[i]=rev[i+1]*(i+1)%p; for (int i=1;i<=MAXN;i++) inv[i]=a[i-1]*rev[i]; } void inv4()//基于a*i+b=p的求法,此做法仅用于p为素数 { inv[1]=1; for (int i=2;i<=MAXN;i++) { int a=p/i; int b=p%i; inv[i]=(p-a)*inv[b]%p; } } int power(int x,int n)//快速幂非递归的 { int ret=1; while (n) { if (n&1) ans*=x; n>>=1; x*=x; } } int powermod(int x,int n,int p)//x^n mod p模幂运算 { if (!n) return 1; int ret=powermod((x*x)%p,n>>1,p); if (n&1) ret=(ret*x)%p; return ret; } int gcd(int a,int b)//欧几里得算法最大公约数 { if (!b) return a; else return (b,a%b); } int lcm(int a,int b)//最小公倍数 { return a*b/gcd(a,b); } int ex_gcd(int a,int b,int &x,int &y)//扩展欧几里得 ax+by=gcd(a,b) { int temp; if (a%b==0) { x=0; y=1; return a; } int ret=ex_gcd(b,a%b,x,y); temp=x; x=y; y=t-a/b*y; return ret; } int CRT(int n,int a[],int mod[])//中国剩余定理 { int temp=1,ret=0,t,y; for (int i=1;i<=n;i++) temp*=mod[i]; for (int i=1;i<=n;i++) { int w=temp/mod[i]; t=ex_gcd(m[i],w,t,y); x=(x+y*w*a[i])%temp; } while (x<=0) x+=temp; return x; }//x=a[1] (mod m[1]) //x=a[2] (mod m[2]) //x=a[3] (mod m[3]) //... int inv(int a,int p)//配合下文食用 { int d,x,y; d=exgcd(a,p,x,y); if (d==1) return (x+p)%p; else return -1; } int C(int n,int m,int p)//普通组合 { if(n<m)return 0; return fac[n]*inv(fac[m]*fac[n-m]%p,p)%p; } int _C(int n,int m,int p)//二逼组合 { if(n<m)return 0; int ans=1,res=1; for(int i=1;i<=n;i++) ans=ans*i%p; for(int i=1;i<=n-m;i++) res=res*i%p; for(int i=1;i<=m;i++) res=res*i%p; res=inv(res,p); ans=(ans*res)%p; return ans; } int Lucas(int n,int m,int p)//文艺组合 Lucas定理求组合数C(n,m)n=a*p+b m=a’*p+b’ C(n,m)=C(a,a’)*C(b,b’) { int ret=1; while (n&&m) { int a=n%p,b=m%p; if (a<b) return 0; ret=(ret*C(a,b,p))%p; n/=p; m/=p; } return ret; } int BSGS(int a,int b,int p)//BSGS算法 a^x=b (mod p) 不想手写hash改成了map { int m=sqrt(p)+.5,temp=inv(power(a,m,p),p),t=1; map<int,int>hash; hash[1]=0; for (int i=1;i<m;i++) { t=t*a%p; hash[t]=i; } for (int i=0;i<=m;i++) { if (hash.count(b)) return i*m+hash[b]; b=b*temp%p; } return -1; } int get_ans(int a,int b,int k)//枚举除法取值,Mobius反演常用 在这之前要维护一个特别的前缀和,内容根据题目而定 { int last=0,ret=0; a/=k;b/=k; for (int i=1;i<=min(a,b);i=last+1) { last=min(a/(a/i),b/(b/i)); ret+=(prev[last]-prev[i-1])*(a/i)*(b/i); } } //容斥原理,结果存在sub中 long long sub = 0; void dfs(int id,int deep,long long sum) { // cout<<sum<<endl; long long temp; for(int i = id; i < primenum; i ++) { temp = sum * prime[i]; // cout<<temp<<endl; if(temp > n) //避免溢出 { return; } if(deep % 2 == 0) { sub -= n / temp; } else { sub += n / temp; } dfs(i + 1,deep + 1,temp); } }