杜教筛 与 数论函数（狄雷克卷积）
为了改变数论只会GCD的尴尬局面，我们来开一波数论：
数论函数：
- 数论函数是定义域在正整数的函数。
- 积性函数： f(ab)=f(a)f(b),gcd(a,b)=1 ，完全积性函数： f(ab)=f(a)f(b) 。
- 常见积性函数： φ(n) ，μ(n) （莫比乌斯函数）， d(n) （因子个数）， σ(n) （因子和）。
- 单位函数 ： e(n)=[n=1] 。
- 常见完全积性函数： Idk(n)=n^k ， 1(n)=Id0(n) ， Id(n)=Id1(n) 。
我们 有以下令人窒息的操作：
(F*G)（x）=∑d|n F（d）*G（n/d）
这种操作我们称之为狄雷克卷积。满足交换律，结合律，分配律以及 f∗e=f 。
我们有: d(n)=(1∗1)(n)
σ(n)=(Id∗1)(n)
e=μ∗1
都可以由积性函数的性质再套一个乘法分配率可以证出来。
我们有两个重要的数论函数 φ(n) ，μ(n) 。
φ(n)，1到n里面与N互质的数的个数。有重要的等式：φ*1=id
μ(n)， 我们可以理解为1的逆元。 我们发现狄雷克卷积可以构成一个群，生成元是e ，1* μ=e，我们又知道*1又叫求该函数的狄雷克前缀和，那么再卷上μ就是原来的函数了。
我们设F=f*1，则 f=F* μ
那么我们就可以做一些简单的题目了：
bzoj 3884 .这道题我们要一个结论：
X^i=X^（i>φ（p）?i mod φ（p）+φ（p）：i） （mod p）
特殊的 ，我们有 X^i=X^(i%(p-1)) (mod p) p为质数。
我们还要知道φ（φ（p））<n/2,那么我们不妨迭代指数，最多迭代log p次指数为1，我们停止迭代。
#include<bits/stdc++.h> using namespace std; int sed(int &p){ int ll=0; while (p%2==0) ll++,p>>=1; return ll; } int fai(int x){ int i,re=x; for(i=2;i*i<=x;i++) if(x%i==0) { re/=i;re*=i-1; while(x%i==0) x/=i; } if(x^1) re/=x,re*=x-1; return re; } long long QP(long long x,int y,int p){ long long re=1; while(y) { if(y&1) (re*=x)%=p; (x*=x)%=p; y>>=1; } return re; } long long dos(int p){ if (p==1) return 1; int k=sed(p), fi=fai(p); long long re=dos(fi); (re+=fi-k%fi)%=fi; re=QP(2,re,p)%p; return (re<<k); } int main () { int t,p; scanf("%d",&t); while (t--) { scanf("%d",&p); printf("%lld\n",dos(p)%p); } }
bzoj 2186
这某非就是传说中的gcd？
#include<bits/stdc++.h> #define inff 10000000 #define N 10000090 int n,m,p,tot,tops,t,T; int ans[N];int fa[N]; bool f[N];char c; int b,inf; using namespace std; inline char nc(){ static char buf[100000],*p1=buf,*p2=buf; return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++; } inline void read(int &x) { c=nc(); b=1; for (; !(c>='0' && c<='9'); c=nc()) if (c=='-') b=-1; for (x=0; c>='0' && c<='9'; x=x*10+c-'0',c=nc()); x*=b; } inline void exgcd(int a,int b,int &x,int &y) { if (!b){x=1;y=0;return;} exgcd(b,a%b,x,y); int t=x;x=y;y=t-a/b*x; } inline int getinv(int a) { int x=0,y=0; exgcd(a,p,x,y); return (x%p+p)%p; } void getfine(){ ans[1]=1;t=0; for (int i=2;i<=inf;i++) { if (!f[i]) ans[i]=(long long)1ll*ans[i-1]*(i-1)%p*getinv(i)%p,fa[++t]=i; else ans[i]=ans[i-1]; for(int j=1;j<=t;j++) { if((long long)fa[j]*i>inf)break; f[fa[j]*i]=1; if(i%fa[j]==0)break; } } } inline void write(int x) { if(x<0) putchar('-'),x=-x; if(x>9) write(x/10); putchar(x%10+'0'); } int main (){ //scanf("%d%d",&T,&p); read(T); read(p); inf=min(p,inff); getfine();fa[0]=1; for (int i=1;i<=inf;i++) fa[i]=(long long)(1ll*fa[i-1]*i)%p; while (T--){ //scanf("%d%d",&n,&m); read(n); read(m); int top=1ll*fa[n]%p*ans[m]%p; //printf("%d\n",top); write(top); putchar('\n'); } }
bzoj 2440
二分答案，我们转为一个计数问题，我们统计有多少符合条件的数小于X？
由容斥原理得 ans=∑ x/(i*i)*μ(i)
#include<bits/stdc++.h> #define N 50000 #define ll long long using namespace std; inline ll read() { ll x=0,f=1;char ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();} while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();} return x*f; } int T; ll ans,K; int tot; int mu[50005],pri[50005]; bool mark[50005]; void getmu() { mu[1]=1; for(int i=2;i<=N;i++) { if(!mark[i])pri[++tot]=i,mu[i]=-1; for(int j=1;j<=tot&&pri[j]*i<=N;j++) { mark[i*pri[j]]=1; if(i%pri[j]==0){mu[i*pri[j]]=0;break;} else mu[i*pri[j]]=-mu[i]; } } } ll cal(int x) { ll sum=0; int t=sqrt(x); for(int i=1;i<=t;i++) sum+=x/(i*i)*mu[i]; return sum; } int main() { getmu(); T=read(); while(T--) { K=read(); ll l=K,r=1644934081; while(l<=r) { ll mid=(l+r)>>1; if(cal(mid)>=K)ans=mid,r=mid-1; else l=mid+1; } printf("%lld\n",ans); } return 0; }
给出T组N和M，依次求出
的值
我们考虑反演：∑i∑j (gcd=1) =∑i∑j ∑d |gcd(i,j) μ(i)=∑ d μ(d) (n/d)*(m/d)
我们发现其只有O(N^0.5)个整点，分块处理。
#include<bits/stdc++.h> #define N 1000000 bool mark[N+5]; int prime[N+5]; int num; long long sum[N+5]; int euler[N+5]; inline int Min(int a,int b){return a<b?a:b;} void Euler() { int i,j; num = 0; memset(euler,0,sizeof(euler)); memset(prime,0,sizeof(prime)); memset(mark,0,sizeof(mark)); euler[1] = 1; // multiply function sum[1]=1; for(i=2;i<=N;i++) { if(!mark[i]) { prime[num++] = i; euler[i] = -1; } for(j=0;j<num;j++) { if(i*prime[j]>N){break;} mark[i*prime[j]] = 1; if(i%prime[j]==0) { euler[i*prime[j]] = 0; break; } else { euler[i*prime[j]] = -euler[i]; } } sum[i]=sum[i-1]+euler[i]; } } int main() { int i; int M1,M2; Euler(); int t; scanf("%d",&t); while (t--){ scanf("%d%d",&M1,&M2); long long ans = 0; int min = Min(M1,M2); /*for(i=1;i<=min;i++) { sum += euler[i]*(M1/i)*(M2/i); }*/ for (int i=1,last;i<=min;i=last+1){ last=Min(M1/(M1/i),M2/(M2/i)); ans+=(sum[last]-sum[i-1])*(M1/i)*(M2/i); } printf("%lld\n",ans); } }
给出T组N和M，求
的值
我们考虑反演：∑i∑j (gcd=1) =∑i∑j ∑d |gcd(i,j)φ (i)=∑ dφ(d) (n/d)*(m/d)
我们发现其只有O(N^0.5)个整点，分块处理。
#include<bits/stdc++.h> #define N 1000000 bool mark[N+5]; int prime[N+5]; int num; long long sum[N+5]; int euler[N+5]; inline int Min(int a,int b){return a<b?a:b;} void Euler() { int i,j; num = 0; memset(euler,0,sizeof(euler)); memset(prime,0,sizeof(prime)); memset(mark,0,sizeof(mark)); euler[1] = 1; // multiply function sum[1]=1; for(i=2;i<=N;i++) { if(!mark[i]) { prime[num++] = i; euler[i] = i-1; } for(j=0;j<num;j++) { if(i*prime[j]>N){break;} mark[i*prime[j]] = 1; if(i%prime[j]==0) { euler[i*prime[j]] = euler[i]*prime[j]; break; } else { euler[i*prime[j]] = euler[i]*euler[prime[j]]; } } sum[i]=sum[i-1]+euler[i]; } } int main() { int i; int M1,M2; Euler(); int t; scanf("%d",&t); while (t--){ scanf("%d%d",&M1,&M2); long long ans = 0; int min = Min(M1,M2); /*for(i=1;i<=min;i++) { sum += euler[i]*(M1/i)*(M2/i); }*/ for (int i=1,last;i<=min;i=last+1){ last=Min(M1/(M1/i),M2/(M2/i)); ans+=(long long)1ll*(sum[last]-sum[i-1])*(M1/i)*(M2/i); } printf("%lld\n",ans); } }
那么我们就可以上杜教筛了。
举个例子：这里 我们要在线性以下的时间搞定φ的前缀和。
(⊙o⊙)… 我放弃了，不会打数学公式。
#include<bits/stdc++.h> #define N 4011007 #define LL long long #define mo 1000000007 #define fg #define int LL using namespace std; LL pim[N>>3],fi[N],tot,a[N],b[N],inv2=500000004,s=N-1007,n,SI,fis[N]; inline LL Mo(LL x){x%=mo;if (x<0) x+=mo;return x;} LL tr(LL x) {if (x<=s) return a[x]; return b[n/x];} void put(LL x,LL y) {if (x<=s) a[x]=y; else b[n/x]=y;} LL F(LL x) { if (x<SI) return fis[x]; if (tr(x)) return tr(x); LL ans=Mo(x)*Mo(x+1)%mo*inv2%mo; for (LL i=2,last;i<=x;i=last+1) { last=x/(x/i); ans=Mo(ans-F(x/i)*Mo(last-i+1)%mo); } put(x,ans);return ans; } void Pin(int MA){ fis[1]=1; for (int i=2;i<MA;i++) { if (!fi[i]) pim[++tot]=i,fi[i]=i-1; for (int j=1;j<=tot&&pim[j]*i<MA;j++) if(i%pim[j]) fi[i*pim[j]]=fi[i]*(pim[j]-1); else {fi[i*pim[j]]=fi[i]*pim[j]; break;} fis[i]=Mo(fis[i-1]+fi[i]); } } void Init() { SI=min((int)powl(n,0.666),(LL)N-12); Pin(SI); } signed main () { scanf("%lld",&n); Init(); printf("%lld",F(n)); }
再丢几道题目：
#include<bits/stdc++.h> #define LL long long #define N 4000007 using namespace std; map<LL,int> ma; LL a,n; int LS,fis[N],pim[N>>3],tot,fi[N],u[N]; bool p[N]; LL F(LL x){ if (x<LS) return fis[x]; if (ma.count(x)) return ma[x]; LL ans=1; for (LL i=2,last;i<=x;i=last+1) { last=x/(x/i); ans=ans-(last-i+1)*F(x/i); } ma[x]=ans; return ans; } void LST(int LS){ fis[1]=1; for (int i=2;i<LS;i++) { if (!p[i]) u[i]=-1,pim[++tot]=i; for (int j=1;j<=tot&&i*pim[j]<LS;j++) { p[i*pim[j]]=1; if (i%pim[j]) u[i*pim[j]]=-u[i]; else{ u[i*pim[j]]=0; break;} } fis[i]=(fis[i-1]+u[i]); } } void init() { LS=min((int)pow(n,0.6667),N-47); LST(LS); } int main () { scanf("%lld %lld",&a,&n); init(); printf("%lld",F(n)-F(a-1)); }
#include<bits/stdc++.h> #define mo 1000000007 #define N 5000007 #define LL long long using namespace std; LL n,fis[N],pim[N>>3],fi[N],IS,inv2=500000004,tot,ans,a[N],b[N],s=N-1007; inline LL Mo(LL x){x%=mo; if (x<0) x+=mo; return x;} //map<LL,LL> Ma; LL tr(LL x) {if (x<=s) return a[x]; return b[n/x];} void put(LL x,LL y) {if (x<=s) a[x]=y; else b[n/x]=y;} LL fiz(LL x){ if (x<=IS) return fis[x]; if (tr(x)) return tr(x); // printf("%d\n",x); LL ans=x%mo*((x+1)%mo)%mo*inv2%mo; for (LL i=2,last;i<=x;i=last+1) { last=x/(x/i); ans=(ans-fiz(x/i)* ((last-i+1)%mo)%mo)%mo; } put(x,ans); return ans; } void P_pit(LL ma_l){ fis[1]=1; for (int i=2;i<=ma_l;i++) { if (!fi[i]) pim[++tot]=i,fi[i]=i-1; for (int j=1;j<=tot&&i*pim[j]<=ma_l;j++) if (i%pim[j]) fi[i*pim[j]]=fi[i]*(pim[j]-1); else {fi[i*pim[j]]=fi[i]*pim[j];break;} fis[i]=fi[i]+fis[i-1]; } } void init() { IS=min(N-17,(int)powl(n,0.6666)); P_pit(IS); fiz(n); } int main () { // freopen("a.in","r",stdin); scanf("%lld",&n); init(); for (LL i=1,last;i<=n;i=last+1) { last=n/(n/i); ans=(ans+Mo(n/i)*Mo(n/i)*Mo(fiz(last)-fiz(i-1))%mo); ans=Mo(ans); } printf("%lld\n",ans); return 0; }
来一发技巧。对于单点杜教筛的话，我们可以这样存取值
LL tr(LL x) {if (x<=s) return a[x]; return b[n/x];}
void put(LL x,LL y) {if (x<=s) a[x]=y; else b[n/x]=y;}
那么我们就不用map 或哈希了。
以后再学洲阁筛吧，挖坑就跑。
