[BZOJ4816][SDOI2017]數字表格(莫比烏斯反演)
4816: [Sdoi2017]數字表格
Time Limit: 50 Sec Memory Limit: 128 MB
Submit: 1259 Solved: 625
[Submit][Status][Discuss]Description
Doris剛剛學習了fibonacci數列。用f[i]表示數列的第i項,那麽 f[0]=0 f[1]=1 f[n]=f[n-1]+f[n-2],n>=2 Doris用老師的超級計算機生成了一個n×m的表格,第i行第j列的格子中的數是f[gcd(i,j)],其中gcd(i,j)表示i, j的最大公約數。Doris的表格中共有n×m個數,她想知道這些數的乘積是多少。答案對10^9+7取模。Input
有多組測試數據。
第一個一個數T,表示數據組數。 接下來T行,每行兩個數n,m T<=1000,1<=n,m<=10^6Output
輸出T行,第i行的數是第i組數據的結果Sample Input
3
2 3
4 5
6 7
Sample Output
1
6
960HINT
Source
鳴謝infinityedge上傳
[Submit][Status][Discuss]
不知道為什麽要用fibonacci,感覺既沒有用到矩乘又沒有用到$gcd(f[i],f[j])=f[gcd(i,j)]$的性質。
首先列出連乘式,可以發現很像莫比烏斯反演,先試著推一下式子,從每個數出現的次數入手。
$$Ans(n,m)=\prod_{d}f(d)^{\sum_{i=1}^{n}\sum_{j=1}^m[(i,j)=d]}=\prod_{d=1}^{\min(n,m)}f(d)^{\sum_{p=1}^{\frac{\min(n,m)}{d}}\mu(d)\lfloor \frac{n}{pd}\rfloor \lfloor\frac{m}{pd}\rfloor}$$
到這裏,有一種錯誤的想法(沒錯就是我誤以為可以拿60分的想法):
這樣$$Ans(n,m)=\prod_{d=1}^{\lfloor \frac{\min(n,m)}{d} \rfloor}f(d)^{g(\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor,d)}$$
這個看上去用分塊優化可以做到$$\begin{aligned}O(T\int_1^n\sqrt{\frac{n}{x}}dx)& =O(T\sqrt{n}\int_1^n\sqrt{\frac{1}{x}}dx)\\ & =O(2T\sqrt{n}\sqrt{n})\\ & =O(Tn)\end{aligned}$$但實際上由於分塊區間不連續,復雜度是不對的。
那麽我們繼續化簡剛才的式子:
$$\begin{aligned}Ans(n,m)& =\prod_{d=1}^{\min(n,m)}f(d)^{\sum_{d|T}^{\min(n,m)}\mu(\frac{T}{d})\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T}\rfloor}\\ & =\prod_{T=1}^{\min(n,m)}\prod_{d|T}f(d)^{\mu(\frac{T}{d})\lfloor \frac{n}{T} \rfloor\lfloor\frac{m}{T}\rfloor}\\ & =\prod_{T=1}^{\min(n,m)}g(T)^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}\end{aligned}$$ $g(T)$可以預處理出來,這樣總復雜度就是$O(n\log n+T\sqrt{\min(n,m)})$了。
註意$\mu$在指數上時會有$-1$,這個要求逆元,如果先預處理所有$f[i]$的逆元的話會快一倍。
這道題總體並不難,主要就是將莫比烏斯反演中的一些套路從加法變成乘法了,實現的時候有幾個細節註意一下就好。
(以後再也不手寫LaTeX莫比烏斯反演的題解了)
1 #include<cstdio> 2 #include<algorithm> 3 #define rg register int 4 #define rep(i,l,r) for (rg i=l; i<=r; i++) 5 typedef long long ll; 6 using namespace std; 7 8 const int N=1000100,mod=1000000007; 9 bool b[N]; 10 int n,m,T,ans,tot,p[N],f[N],g[N],G[N],G1[N],miu[N],F[N][3]; 11 12 int ksm(int a,int b){ 13 int res; 14 for (res=1; b; a=1ll*a*a%mod,b>>=1) 15 if (b & 1) res=1ll*res*a%mod; 16 return res; 17 } 18 19 void pre(){ 20 miu[1]=1; 21 for (rg i=2; i<N; i++){ 22 if (!b[i]) p[++tot]=i,miu[i]=-1; 23 for (rg j=1; j<=tot && p[j]*i<N; j++){ 24 rg t=p[j]*i; b[t]=1; 25 if (i%p[j]) miu[t]=-miu[i]; else break; 26 } 27 } 28 for (rg i=1; i<N; i++) g[i]=1,F[i][0]=ksm(f[i],mod-2),F[i][1]=1,F[i][2]=f[i]; 29 for (rg i=1; i<N; i++) 30 for (rg j=i; j<N; j+=i) 31 g[j]=1ll*g[j]*F[i][miu[j/i]+1]%mod; 32 G[0]=1; G[1]=g[1]; for (int i=2; i<N; i++) G[i]=1ll*G[i-1]*g[i]%mod; 33 G1[N-1]=ksm(G[N-1],mod-2); for (int i=N-2; ~i; i--) G1[i]=1ll*G1[i+1]*g[i+1]%mod; 34 } 35 36 int main(){ 37 freopen("product.in","r",stdin); 38 freopen("product.out","w",stdout); 39 f[1]=1; for (rg i=2; i<N; i++) f[i]=(f[i-1]+f[i-2])%mod; 40 pre(); 41 for (scanf("%d",&T); T--; ){ 42 scanf("%d%d",&n,&m); rg lst=0; 43 if (n>m) swap(n,m); ans=1; 44 for (rg i=1; i<=n; i=lst+1){ 45 lst=min(n/(n/i),m/(m/i)); 46 ans=1ll*ans*ksm(1ll*G[lst]*G1[i-1]%mod,1ll*(n/i)*(m/i)%(mod-1))%mod; 47 } 48 printf("%d\n",ans); 49 } 50 return 0; 51 }
[BZOJ4816][SDOI2017]數字表格(莫比烏斯反演)