BZOJ 4816 [SDOI2017]数字表格 (莫比乌斯反演)

题目大意：略

洛谷传送门

第一次不看题解做出来反演题，我好菜啊

最初推的式子里是把矩阵快速幂然后欧拉降幂，然后发现弄出来的矩阵不能直接相乘..

但这依然改变不了反演玩的是套路的事实

以下默认$n \leq m$

题目让我们求$\prod\limits_{i=1}^{n} \prod\limits_{j=1}^{m}f[gcd(i,j)]$

显然直接统计$f[gcd(i,j)]$很困难，我们每次指定一个值k，求有多少个$gcd(i,j)==k$，那么$k$的贡献就是$f[k]^{ans}$

变形$\prod\limits_{k=1}^{n} f[k]^{\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m}[gcd(i,j)==k]}$

我们只需要处理出$\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m}[gcd(i,j)==k]$的数量即可，这是一道经典的反演题

再利用常规套路，$\varepsilon =1*\mu$变形，$\sum\limits_{i=1}^{\left \lfloor \frac{n}{k} \right \rfloor} \sum\limits_{j=1}^{\left \lfloor \frac{m}{k} \right \rfloor}\sum\limits_{d|gcd(i,j)}\mu(d)$

指数最后变成了$\sum\limits_{d=1}^{\left \lfloor \frac{n}{k} \right \rfloor}\mu(d) \left \lfloor \frac{n}{kd} \right \rfloor \left \lfloor \frac{m}{kd} \right \rfloor$

原式$\prod\limits_{k=1}^{n} f[k]^{\sum\limits_{d=1}^{\left \lfloor \frac{n}{k} \right \rfloor}\mu(d) \left \lfloor \frac{n}{kd} \right \rfloor \left \lfloor \frac{m}{kd} \right \rfloor}$

还是常规套路提取出$Q=kd$，变形$\prod\limits_{Q=1}^{n}  \prod\limits_{d|Q}f[d]^{\mu(\frac{Q}{d}) \left \lfloor \frac{n}{Q} \right \rfloor \left \lfloor \frac{m}{Q} \right \rfloor}$

令$g(Q)=\prod\limits_{d|Q} f[d]^{\mu(\frac{Q}{d})}$，显然它的前$n$项可以在$O(nlogn)$内处理出

每次询问都需要整除分块，还需要快速幂。总时间$O(T\sqrt{n}logn+nlogn)$

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <algorithm>
 4 #define N1 1010000
 5 #define ll long long
 6 #define dd double
 7 #define maxn 1000000
 8 using namespace std;
 9 
10 const ll mod2=1000000007;
11 int T,n,m;
12 int use[N1],pr[N1],mu[N1],cnt;
13 int qpow(ll x,ll y,const ll &mod)
14 {
15     int ans=1;
16     while(y){
17         if(y&1) ans=ans*x%mod;
18         x=x*x%mod; y>>=1;
19     }return ans;
20 }
21 void exgcd(ll a,ll b,ll &x,ll &y)
22 {
23     if(!b){ x=1; y=0; return;}
24     else{ exgcd(b,a%b,x,y); ll t=x; x=y; y=t-a/b*y; }
25 }
26 ll f[N1],g[N1],sg[N1];
27 void init()
28 {
29     int i,j; ll x,y,inv; mu[1]=1;
30     for(i=2;i<=maxn;i++)
31     {
32         if(!use[i]){ pr[++cnt]=i; mu[i]=-1; }
33         for(j=1;j<=cnt&&i*pr[j]<=maxn;j++)
34         {
35             use[i*pr[j]]=1;
36             if(i%pr[j]){ mu[i*pr[j]]=-mu[i]; }
37             else{ mu[i*pr[j]]=0; break; }
38         }
39     }
40     for(j=1;j<=maxn;j++) g[j]=1;
41     for(f[1]=1,i=2;i<=maxn;i++)
42     {
43         f[i]=(f[i-1]+f[i-2])%mod2;
44         exgcd(f[i],mod2,inv,y); inv=(inv%mod2+mod2)%mod2;
45         for(j=1;j*i<=maxn;j++)
46             if(mu[j]) g[i*j]=1ll*g[i*j]*(mu[j]>0?f[i]:inv)%mod2;
47     }
48     for(sg[0]=1,i=1;i<=maxn;i++) sg[i]=1ll*sg[i-1]*g[i]%mod2;
49 }
50 
51 int main()
52 {
53     scanf("%d",&T);
54     init();
55     int i,j,la,ans=0; ll inv,y;
56     while(T--)
57     {
58         scanf("%d%d",&n,&m); if(n>m) swap(n,m); 
59         for(i=1,ans=1;i<=n;i=la+1)
60         {
61             la=min(n/(n/i),m/(m/i));
62             exgcd(sg[i-1],mod2,inv,y); inv=(inv%mod2+mod2)%mod2;
63             ans=1ll*qpow(1ll*sg[la]*inv%mod2,1ll*(n/i)*(m/i),mod2)*ans%mod2;
64         }
65         printf("%d\n",ans);
66     }
67     return 0;
68 }

 

转载于:https://www.cnblogs.com/guapisolo/p/10222553.html