/**********/ #include<iostream> #include<cstdio> #include<cmath> #include<algorithm> #include<cstring> #include<string> #include<cstdlib> #include<vector> #include<stack> #include<map> using namespace std; typedef long long ll; /***********GCD**********************/ ll gcd(ll a,ll b) { if(b==0) return a; return gcd(b,a%b); } /*************快速乘***************/ ll mult_mod(ll a,ll b,ll mod) { a%=mod; b%=mod; ll res=0; while(b) { if(b&1) { res+=a; res%=mod; } a<<=1; if(a>=mod) a%=mod; b>>=1; } return res; } /*************快速幂**************/ ll pow_mod(ll x,ll n,ll mod) { if(n==1) return x%mod; x%=mod; ll t=x,res=1; while(n) { if(n&1) res=mult_mod(res,t,mod); t=mult_mod(t,t,mod); n>>=1; } return res; } /*****更快的*****快速幂***************/ ll pow_mod(ll x,ll n,ll mod) { ll res=1; ll tmp=x%mod; while(n) { if(n&1) { res*=t; res%=mod; } n>>=1; t*=t; t%=mod; } return res; } /************扩展欧几里德****************/ void extend_gcd(ll a,ll b,ll &x,ll &y) { ll d; //d=gcd(a,b) if(b==0) { x=1,y=0; d=a; } else { d=extend_gcd(b,a%b,y,x); ll xx=x,yy=y; x=yy; y=xx-(a/b)*yy; } } /************素数筛法******************/ ll l,u,prime[N]; int tot; int vis[N],ans[10000005]; void isPrime() { tot=0; memset(vis,0,sizeof(vis)); memset(prime,0,sizeof(prime)); for(ll i=2;i<N;i++) { if(!vis[i]) { prime[tot++]=i; for(ll j=i*i;j<N;j+=i) vis[j]=1; } } } /******更快的*素数筛选*******************/ const int MAXN=10000; int prime[MAXN+1]; void isPrime() { memset(prime,0,sizeof(prime)); //素数存储从下标1开始，prime[0]记录素数个数 for(int i=2;i<=MAXN;i++) { if(!prime[i]) prime[++prime[0]]=i; for(int j=1;j<=prime[0]&&prime[j]<=MAXN/i;j++) { prime[prime[j]*i]=1; if(i%prime[j]==0) break; } } } /*********素数判定Miller_Rabin***********/ bool test(ll a,ll n) //Miller_Rabin算法的核心 { ll x=n-1,t=0,res,last; while((x&1)==0) { x>>=1; t++; } last=pow_mod(a,x,n); for(int i=0;i<t;i++) { res=pow_mod(last,2,n); if(res==1&&last!=1&&last!=n-1) return true; last=res; } if(res!=1) return true; return false; } bool millier_rabin(ll n) { if(n==2) return true; if(n==1||(n&1)==0) return false; for(int i=0;i<times;i++) { ll a=rand()%(n-1)+1; if(test(a,n)) return false; } return true; } /*******质因数分解********************/ ll factor[100][2]; int cnt; //分解质因数:A=(p1^a1)*(p2^a2)*...*(pn^an) //factor[i][0]=pi,factor[i][1]=ai; void getFactor(ll x) { cnt=0; ll t=x; for(int i=0;prime[i]<=t/prime[i];i++) { factor[cnt][1]=0; while(t%prime[i]==0) { factor[cnt][0]=prime[i]; while(t%prime[i]==0) { factor[cnt][1]++; t/=prime[i]; } cnt++; } } if(t!=1) { factor[cnt][0]=t; factor[cnt][1]=1; cnt++; } } /*********大数因数分解Pollard_rho*************/ ll pollard_rho(ll x,ll c) { ll x0,y,i=1,k=2; x0=rand()%x; y=x0; while(1) { i++; x0=(mult_mod(x0,x0,x)+c)%x; ll d=gcd(y-x0,x); if(d>1&&d<x) return d; if(y==x0) break; if(i==k) { y=x0; k+=k; } } return x; } void find_fac(ll n,int c) { if(n==1) return; if(millier_rabin(n)) { factor[tot++]=n; return; } ll p=n; while(p>=n) p=pollard_rho(p,c--); find_fac(p,c); find_fac(n/p,c); } /**********正约数和（二分求法）**************/ //将A进行因式分解可得:A = p1^a1 * p2^a2 * p3^a3 * pn^an //A^B=p1^(a1*B)*p2^(a2*B)*...*pn^(an*B); //A^B的所有约数之和sum=[1+p1+p1^2+...+p1^(a1*B)]*[1+p2+p2^2+...+p2^(a2*B)]*[1+pn+pn^2+...+pn^(an*B)]. //等比数列1+pi+pi^2+pi^3+...+pi^n可以由二分求得(即将需要求解的因式分成部分来求解) //若n为奇数,一共有偶数项,则 //1+p+p^2+p^3+........+p^n=(1+p+p^2+....+p^(n/2))*(1+p^(n/2+1)); //若n为偶数,一共有奇数项,则 //1+p+p^2+p^3+........+p^n=(1+p+p^2+....+p^(n/2-1))*(1+p^(n/2+1))+p^n/2; ll sum(ll p,ll n) { if(p==0) return 0; if(n==0) return 1; if(n&1) return ((1+pow_mod(p,n/2+1))%MOD*sum(p,n/2)%MOD)%MOD; else return ((1+pow_mod(p,n/2+1))%MOD*sum(p,n/2-1)+pow_mod(p,n/2)%MOD)%MOD; } /**********欧拉函数**********************/ int Euler(int n) { int res=n; for(int i=2;i*i<=n;i++) { while(n%i==0) { n/=i; res-=(res/i); while(n%i==0) n/=i; } } if(n>1) res-=(res/n); return res; } /******素数筛法+欧拉打表**********/ ll e[N+5],p[N+5]; bool vis[N+5]; void init() { memset(e,0,sizeof(e)); memset(p,0,sizeof(p)); memset(vis,false,sizeof(vis)); ll i,j; p[0]=1;//记录素数总个数 p[1]=2; for(i=3;i<N;i+=2) { if(!vis[i]) { p[++p[0]]=i; for(j=i*i;j<N;j+=i) vis[j]=true; } } e[1]=1; for(i=1;i<=p[0];i++) e[p[i]]=p[i]-1; //为什么要-1？ for(i=2;i<N;i++) { if(!e[i]) { for(j=1;j<=p[0];j++) { if(i%p[j]==0) { if(i/p[j]%p[j]) e[i]=e[i/p[j]]*e[p[j]]; else e[i]=e[i/p[j]]*p[j]; break; } } } } } /*************欧拉打表*更快版***************/ #define N 1000000 int e[N+5]; void init() { int i,j; memset(e,0,sizeof(e)); for(i=2;i<=N;i++) { if(!e[i]) { for(j=i;j<=N;j+=i) { if(!e[j]) e[j]=j; e[j]=e[j]/i*(i-1); } } } } /**********容斥原理*二进制解法**************/ vector<int> v; ll solve(int l,int r,int n) //[l,r]内与n互素的数字个数 { v.clear(); //筛选素数 for(int i=2;i*i<=n;i++) { if(n%i==0) { v.push_back(i); while(n%i==0) n/=i; } } if(n>1) v.push_back(n); //容斥原理的二进制解法 int len=v.size(); ll res=0; for(int i=1;i<(1<<len);i++) { int cnt=0; ll val=1; for(int j=0;j<len;j++) { if(i&(1<<j)) { cnt++; val*=v[j]; } } if(cnt&1) //若为奇数项进行加法，偶数项进行减法 res+=r/val-(l-1)/val; else res-=r/val-(l-1)/val; } return r-l+1-res; }