反演曾经一直是我不敢搞的一个大坑……
又重新学习了一下反演，并且做了一些习题……
大概基础什么的……我就介绍一点常用的
正经反演的式子有这样两种
$$f(n)=sum _{d|n}F(d) mu(frac{n}{d})$$
以及
$$f(n)=sum _{n|d} F(d)mu(frac{d}{n})$$
很多题都要用关于狄利克雷卷积的一个式子
$mu otimes 1 = e$，也即$sum _{d|n} mu(d)=[n==1]$
证明的话……我们可以用二项式定理来证，网上有很多证明，我就不证了
基本上，所有题目的转换都会利用这个式子
我们一般会搞一个$[???==1]$，然后套用上面那个关于$summu$的式子
…………板子题就不放了，每道题都挺水的
一般就是枚举$gcd$，然后判$[gcd==d]$的时候才累加，
然后除掉一个$d$变成$[gcd==1]$，然后套式子.
写一些比较有意思的题目好了。
这题反演的式子还是可以推的……甚至说比较板子
比较好的思想是我们把询问离线按a排序，下标i按照F值排序
然后用树状数组维护答案并且更新
数学与数据结构的优秀结合……
贴代码：
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 using namespace std; 5 typedef long long LL; 6 const int N=100010; 7 const unsigned int inf=0x7fffffff; 8 bool vis[N+10];int tot; 9 unsigned int prime[N],sumd[N+10]; 10 unsigned int tmp1[N+10],tmp2[N+10],mu[N+10],p[N+10]; 11 struct node{unsigned int ans,id,n,m,a;}q[20010]; 12 inline bool mt1(const node &a,const node &b){return a.a<b.a;} 13 inline bool mt2(const int &a,const int &b){return sumd[a]<sumd[b];} 14 inline bool mt3(const node &a,const node &b){return a.id<b.id;} 15 inline void intn() 16 { 17 int k;mu[1]=sumd[1]=1; 18 for(int i=2;i<=N;i++) 19 { 20 if(!vis[i]) 21 { 22 prime[++tot]=tmp2[i]=i; 23 sumd[i]=tmp1[i]=i+1,mu[i]=-1; 24 } 25 for(int j=1;j<=tot&&(k=i*prime[j])<=N;j++) 26 { 27 vis[k]=1; 28 if(i%prime[j]==0) 29 { 30 tmp2[k]=tmp2[i]*prime[j],tmp1[k]=tmp1[i]+tmp2[k]; 31 sumd[k]=sumd[i]/tmp1[i]*tmp1[k]; 32 mu[k]=0;break; 33 } 34 sumd[k]=sumd[i]*sumd[prime[j]]; 35 tmp1[k]=1+prime[j];tmp2[k]=prime[j]; 36 mu[k]=-mu[i]; 37 } 38 } 39 for(int i=1;i<=N;i++)p[i]=i; 40 sort(p+1,p+N+1,mt2); 41 } 42 unsigned int bit[N+10]; 43 inline int lowbit(int a){return a&-a;} 44 inline void add(int pos,unsigned int val){while(pos<=N)bit[pos]+=val,pos+=lowbit(pos);} 45 inline unsigned int query(int pos) 46 {unsigned int ret=0;while(pos)ret+=bit[pos],pos-=lowbit(pos);return ret;} 47 inline unsigned int calc(unsigned int n,unsigned int m) 48 { 49 if(n>m)swap(n,m); 50 unsigned int ret=0,last; 51 for(unsigned int i=1;i<=n;i=last+1) 52 { 53 last=min(n/(n/i),m/(m/i)); 54 ret+=(n/i)*(m/i)*(query(last)-query(i-1)); 55 } 56 return ret; 57 } 58 int main() 59 { 60 intn();int top=1,t;scanf("%d",&t); 61 for(int i=1;i<=t;i++) 62 q[i].id=i,scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a); 63 sort(q+1,q+t+1,mt1); 64 for(int i=1;i<=t;i++) 65 { 66 while(top<=N&&sumd[p[top]]<=q[i].a) 67 { 68 for(int k=1;k*p[top]<=N;k++) 69 add(k*p[top],mu[k]*sumd[p[top]]); 70 top++; 71 } 72 q[i].ans=calc(q[i].n,q[i].m); 73 } 74 sort(q+1,q+t+1,mt3); 75 for(int i=1;i<=t;i++) 76 printf("%d ",q[i].ans&inf); 77 }
我们知道，题给的$f$是一个积性函数
然后我们枚举一下$gcd$，则有
$ans=sum _{i=1} ^{n} sum _{j=1} ^{m} f((i,j))$
$=sum _{d=1} ^{n} f(d)sum _{i=1} ^{n} sum _{j=1} ^{m} [(i,j)==d]$
$=sum _{d=1} ^{n} f(d)sum _{i=1} ^{left lfloor frac {n}{d} ight floor } sum _{j=1} ^{left lfloor frac {m}{d} ight floor } [(i,j)==1]$
$=sum _{d=1} ^{n} f(d) sum _{g=1} ^{ left lfloor frac {n}{d} ight floor } mu(g) left lfloor frac {n}{dg} ight floor left lfloor frac {m}{dg} ight floor$
然后我们转而枚举$T=dg$
则原式有
$sum _{T=1} ^{n} left lfloor frac {n}{T} ight floor left lfloor frac {m}{T} ight floor sum _{d|T}f(d)mu(frac{T}{d})$
那么我们如果能处理内层函数$sum _{d|T}f(d)mu(frac{T}{d})$的前缀和，就能$O(sqrt{n})$解决这题了
然后……我们来考察一下内层函数怎么求
我们当然可以$O(nlogn)$暴力预处理
但是不知道能不能过……由于$f$和$mu$都是积性函数，因此它们的狄利克雷卷积也是积性函数
我们设$g(T)=sum _{d|T}f(d)mu(frac{T}{d})$
我们设$T=prod p_{i} ^ { a_{i} }$
那么$frac{T}{d}$中每个$p_{i}$都只有0个或者1个（多了$mu$就是0了）
又由于$f(d)$与最大的$a_{i}$有关，我们考虑按照a_{i}讨论
我们设$a_{i}$最大的因子集合为A，不最大的为B
如果存在$i eq j$，有$a_{i} eq a_{j}$
不管A集合如何选择，B集合的每个因子都有选和不选2种选择，从而使得$mu$值一正一负，最后值为0
即，存在$i eq j,a_{i} eq a_{j}$时，$g(T)=0$
反之，如果对于任意$i eq j$，有$a_{i}=a_{j}$
那么我们所有的因子都属于A集合，每个因子都是有选和不选2种选择，从而使得$mu$值一正一负。
但是对于另外一个f函数来说，如果$mu$中出现了全部因子，那么$f=a-1$，否则$f=a$
因此此时，$g(T)=-1*(-1)^{k}=(-1)^{k+1}$,k为T的质因子个数
有了这些……我们就可以写线性筛的式子了
线性筛挺简单的……自己根据关系推一下就行
这题反演的式子好说，关键是构造函数线筛的分析。
这个分析很妙……结合了函数的特点进行了讨论。
放个代码：
1 #include <cstdio> 2 #include <cstring> 3 using namespace std; 4 #define N 10000010 5 #define K 10000000 6 #define LL long long 7 char B[1<<15],*S=B,*T=B; 8 #define getc (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++) 9 inline int read() 10 { 11 int x=0;register char c=getc; 12 while(c<'0'||c>'9')c=getc; 13 while(c>='0'&&c<='9')x=10*x+(c^48),c=getc; 14 return x; 15 } 16 int prime[K/10],tot,vis[N]; 17 int mincnt[N],g[N],ppow[N]; 18 inline void intn() 19 { 20 register int i,j,tmp; 21 for(i=2;i<=K;++i) 22 { 23 if(!vis[i])prime[++tot]=i,ppow[i]=mincnt[i]=g[i]=1; 24 for(j=1;j<=tot&&(tmp=i*prime[j])<=K;++j) 25 { 26 vis[tmp]=1; 27 if(i%prime[j]==0) 28 { 29 mincnt[tmp]=mincnt[i]+1,ppow[tmp]=ppow[i]; 30 if(1==ppow[tmp])g[tmp]=1; 31 else g[tmp]=(mincnt[ppow[tmp]]==mincnt[tmp])?-g[ppow[tmp]]:0; 32 break; 33 } 34 mincnt[tmp]=1,ppow[tmp]=i, 35 g[tmp]=(mincnt[i]==1)?-g[i]:0; 36 } 37 } 38 for(i=1;i<=K;++i)g[i]+=g[i-1]; 39 } 40 inline int min(int a,int b){return a<b?a:b;} 41 int main() 42 { 43 // freopen("Ark.in","r",stdin); 44 register int t,n,m,i,j; 45 t=read();LL ans=0;intn(); 46 while(t--) 47 { 48 n=read(),m=read(); 49 if(n>m)n^=m,m^=n,n^=m; 50 for(i=1,ans=0;i<=n;i=j+1) 51 j=min(n/(n/i),m/(m/i)),ans+=(n/i)*1ll*(m/i)*(g[j]-g[i-1]); 52 printf("%lld ",ans); 53 } 54 }
一开始依然使用枚举$gcd$的套路
过程不写了，结果大概是：
$prod_{d=1}^{n} fibo(d) ^ { sum_{g=1} ^ { left lfloor frac{n}{d} ight floor } mu(g) left lfloor frac{n}{dg} ight floorleft lfloor frac{m}{dg} ight floor }$
然后我们枚举$T=dg$，有
$prod_{T=1}^{n} ( prod_{d|t} fibo(d) ^ { mu( frac{T}{d} ) } ) ^{left lfloor frac{n}{dg} ight floor left lfloor frac{m}{dg} ight floor} $
然后我们设$g(T)=prod_{d|t} fibo(d) ^ { mu( frac{T}{d} ) } $
只要预处理一下这个g函数的前缀积以及其逆元，我们就能除法分块求这个问题了。
然后呢……我们可以用组合数时候那种线性求逆元来搞，就是用费马小定理求出最后一项逆元，然后在用乘除关系线性回推
但是直接求也不是很好求……我们处理一下斐波那契数列的值及其逆元，然后用$nln_{n}$的枚举倍数处理出来每个$g(T)$
然后再用一遍线性求逆元……
这代码只能说，短但是步骤很多……
这个线性求逆元的想法很好！
贴个代码：
1 #include <cstdio> 2 #include <cstring> 3 using namespace std; 4 #define N 1000010 5 #define K 1000000 6 #define LL long long 7 #define mod 1000000007 8 char B[1<<15],*S=B,*T=B; 9 #define getc (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++) 10 inline int read() 11 { 12 int x=0;register char c=getc; 13 while(c<'0'||c>'9')c=getc; 14 while(c>='0'&&c<='9')x=10*x+(c^48),c=getc; 15 return x; 16 } 17 inline LL quick_mod(LL di,LL mi) 18 { 19 LL ans=1; 20 for(di%=mod;mi;mi>>=1,di=di*di%mod) 21 if(mi&1)ans=ans*di%mod; 22 return ans; 23 } 24 LL f[N],g[N],invf[N],invg[N],h[N],x[N],invx[N]; 25 int prime[N/10],tot,vis[N],mu[N]; 26 inline void intn() 27 { 28 register int i,j,tmp; 29 for(f[0]=0,f[1]=g[1]=1,i=2;i<=K;++i) 30 f[i]=(f[i-1]+f[i-2])%mod,g[i]=g[i-1]*f[i]%mod; 31 for(invg[K]=quick_mod(g[K],mod-2),i=K;i;--i) 32 h[i]=1,invg[i-1]=invg[i]*f[i]%mod,invf[i]=invg[i]*g[i-1]%mod; 33 for(mu[1]=1,i=2;i<=K;++i) 34 { 35 if(!vis[i])prime[++tot]=i,mu[i]=-1; 36 for(j=1;j<=tot&&(tmp=prime[j]*i)<=K;++j) 37 { 38 vis[tmp]=1; 39 if(i%prime[j]==0){mu[tmp]=0;break;} 40 mu[tmp]=-mu[i]; 41 } 42 } 43 for(i=2;i<=K;++i) 44 for(j=1;(tmp=i*j)<=K;++j) 45 if(mu[j]==1)h[tmp]=h[tmp]*f[i]%mod; 46 else if(mu[j]==-1)h[tmp]=h[tmp]*invf[i]%mod; 47 for(x[1]=h[1]=1,i=2;i<=K;++i)x[i]=x[i-1]*h[i]%mod; 48 invx[K]=quick_mod(x[K],mod-2); 49 for(i=K;i;--i)invx[i-1]=invx[i]*h[i]%mod; 50 } 51 inline int min(int a,int b){return a<b?a:b;} 52 int main() 53 { 54 register int i,j,t,n,m; 55 t=read();intn();LL ans,tmp; 56 while(t--) 57 { 58 n=read(),m=read(); 59 if(n>m)n^=m,m^=n,n^=m; 60 for(ans=i=1;i<=n;i=j+1) 61 { 62 j=min(n/(n/i),m/(m/i)), 63 ans=ans*quick_mod( x[j]*invx[i-1]%mod , (n/i)*1ll*(m/i)%(mod-1) )%mod; 64 } 65 printf("%lld ",ans); 66 } 67 }
……让人无法形容的恐惧
我只能说恐惧……恐惧……
其实，反演的部分还是可以推式子的
但是最后的那一下……简直算是神来之笔
谁会去猜测$sum_{i=1}^{n} i^{k}$可以写成一个$k+1$次的多项式啊！
然后我们拿到这个多项式之后才能继续推式子……
或者……我们可以用更加令人懵逼的伯努利数
感觉到了自己的愚蠢
附上蒟弱的代码
1 #include <cstdio> 2 #include <cstring> 3 using namespace std; 4 #define LL long long 5 #define N 1010 6 #define D 110 7 #define mod 1000000007 8 char cB[1<<15],*S=cB,*T=cB; 9 #define getc (S==T&&(T=(S=cB)+fread(cB,1,1<<15,stdin),S==T)?0:*S++) 10 inline int read() 11 { 12 int x=0;register char c=getc; 13 while(c<'0'||c>'9')c=getc; 14 while(c>='0'&&c<='9')x=10*x+(c^48),c=getc; 15 return x; 16 } 17 inline int quick_mod(int di,int mi) 18 { 19 if(mi<0)return quick_mod(quick_mod(di,mod-2),-mi); 20 int ret=1; 21 for(di%=mod;mi;mi>>=1,di=(LL)di*di%mod) 22 if(mi&1)ret=(LL)ret*di%mod; 23 return ret; 24 } 25 int d,w,p[N],A,a[N]; 26 int fac[D],inv[D],invf[D],B[D]; 27 #define C(i,j) ((LL)fac[i]*invf[j]%mod*invf[i-j]%mod) 28 int main() 29 { 30 // freopen("Ark.in","r",stdin); 31 register int i,j,k,ans,sum,tot=1; 32 d=read(),w=read(); 33 for(i=1;i<=w;++i) 34 p[i]=read(),a[i]=read(),tot=(LL)tot*quick_mod(p[i],a[i])%mod; 35 for(fac[0]=fac[1]=1,i=2;i<=d+2;++i)fac[i]=(LL)fac[i-1]*i%mod; 36 for(inv[0]=inv[1]=1,i=2;i<=d+2;++i)inv[i]=(LL)(mod-mod/i)*inv[mod%i]%mod; 37 for(invf[0]=invf[1]=1,i=2;i<=d+2;++i)invf[i]=(LL)invf[i-1]*inv[i]%mod; 38 for(B[0]=1,i=1;i<=d+2;++i) 39 { 40 for(B[i]=0,j=0;j<i;++j) 41 B[i]=(B[i]+(LL)C(i+1,j)*B[j]%mod)%mod; 42 B[i]=((LL)-B[i]*inv[i+1]%mod+mod)%mod; 43 } 44 for(ans=0,i=1;i<=d+1;++i) 45 { 46 A=((LL)( ((d+1-i)&1)?-1:1 )*C(d+1,i)*B[d+1-i]%mod*inv[d+1]%mod+mod)%mod; 47 if(A==0)continue; 48 // printf("%d ",A); 49 for(sum=quick_mod(tot,i),j=1;j<=w;++j) 50 sum=(LL)sum*( mod+ 1ll- quick_mod(p[j],d-i) ) %mod; 51 // printf("sum=%d ",sum); 52 ans=(ans+(LL)A*sum%mod)%mod; 53 } 54 printf("%d ",(ans%mod+mod)%mod); 55 }
我觉得……我2天打不完式子
我还是……粘dalao的题解过来吧
但是有几个要点可以总结
1.$left lfloor frac{nk}{m} ight floor = frac{nk-nk\%m}{m} = frac{nk}{m}-frac{nk\%m}{m}$
也就是说我们可以把m的整数倍提出来
2.这道题大量的运用了把问题拆分再分段处理的思想……这样的确可以简化问题，这种不行就拆分的转换是很巧妙的
代码：
1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 using namespace std; 5 #define mod 998244353 6 #define N 500010 7 #define K 500000 8 #define LL long long 9 #define db double 10 int n,m,x;double tmp; 11 //把整除拆成 12 inline int Sum(int x){return ((LL)x*(x+1)>>1)%mod;} 13 int prime[N/10],tot,vis[N],mu[N]; 14 inline void intn(int m) 15 { 16 register int i,j,tmp; 17 for(mu[1]=1,i=2;i<=m;++i) 18 { 19 if(!vis[i])prime[++tot]=i,mu[i]=-1; 20 for(j=1;j<=tot&&(tmp=i*prime[j])<=m;++j) 21 { 22 vis[tmp]=1; 23 if(i%prime[j]==0){mu[tmp]=0;break;} 24 mu[tmp]=-mu[i]; 25 } 26 } 27 for(i=2;i<=m;++i)mu[i]+=mu[i-1]; 28 } 29 inline int min(int a,int b){return a<b?a:b;} 30 inline int calc(int a,int b) 31 { 32 register int i,j,ret=0; 33 if(a>b)a^=b,b^=a,a^=b; 34 for(i=1;i<=a;i=j+1) 35 j=min(a/(a/i),b/(b/i)),ret=(ret+(LL)(mu[j]-mu[i-1])*(a/i)*(b/i)%mod)%mod; 36 return ret; 37 } 38 int main() 39 { 40 // freopen("Ark.in","r",stdin); 41 register int i,j,k; 42 scanf("%d%d%lf",&n,&m,&tmp);x=tmp; 43 int ans=((LL)Sum(n)*Sum(m)-(LL)Sum(n)*m-(LL)Sum(m)*n)%mod; 44 if(n>m)n^=m,m^=n,n^=m;intn(n); 45 for(i=1;i<=n;++i) 46 ans=(ans+(LL)(2*i*(x/i)+i)*calc(n/i,m/i)%mod)%mod; 47 ans=(ans+mod)%mod,printf("%lld ",(LL)ans*499122177%mod); 48 }