BZOJ 4816 SDOI2017 数字表格 莫比乌斯反演

题目链接:http://www.lydsy.com/JudgeOnline/problem.php?id=4816

 

题意概述:

多组询问,给出一个n*m的表格,表格中的格子form[i][j]=f[gcd(i,j)],f[i]表示fibonacci数列的第i项。

求表格中所有数字的乘积对 10^9+7 取模的结果。、

T<=1000,1<=n,m<=10^6.

 

分析:

一眼就可以认出来这是莫比乌斯反演。关键在于式子的推倒。

 说一下推式子的要点,主要是能够灵活地利用反演公式,把式子踢来踢去!

啊这个博客实在是不好打公式......(至少我还没有研究出来),具体过程请看ATP的博客!(https://blog.csdn.net/fromatp/article/details/69944236)

 

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 #include<cstdlib>
 5 #include<algorithm>
 6 #include<cmath>
 7 #include<queue>
 8 #include<set>
 9 #include<map>
10 #include<vector>
11 #include<cctype>
12 using namespace std;
13 const int mo=1000000007;
14 const int maxp=80005;
15 const int maxn=1000005;
16 typedef long long LL;
17 
18 int T,N;
19 int f[maxn],pri[maxp],mu[maxn],pcnt,g[maxn],_g[maxn];
20 bool ntp[maxn];
21 struct data{ int n,m; }Q[1005];
22 
23 void exgcd(LL a,LL b,LL &d,LL &x,LL &y){
24     if(!b) d=a,x=1,y=0;
25     else exgcd(b,a%b,d,y,x),y-=a/b*x;
26 }
27 int inv(int a){
28     LL d,x,y;
29     exgcd(a,mo,d,x,y);
30     return (x%mo+mo)%mo;
31 }
32 int qkpow(int x,LL y)
33 {
34     int re=1,tt=x;
35     for(int i=0;(1ll<<i)<=y;i++){
36         if((1ll<<i)&y) re=1ll*re*tt%mo;
37         tt=1ll*tt*tt%mo;
38     }
39     return re;
40 }
41 void ready()
42 {
43     ntp[0]=ntp[1]=1,mu[1]=1;
44     for(int i=2;i<=N;i++){
45         if(!ntp[i]) pri[++pcnt]=i,mu[i]=-1;
46         for(int j=1;j<=pcnt&&1ll*i*pri[j]<=N;j++){
47             ntp[i*pri[j]]=1;
48             if(i%pri[j]==0){ mu[i*pri[j]]=0; break; }
49             mu[i*pri[j]]=-mu[i];
50         }
51     }
52     f[0]=0,f[1]=1,g[0]=_g[0]=g[1]=_g[1]=1;
53     for(int i=2;i<=N;i++){
54         f[i]=(f[i-1]+f[i-2])%mo;
55         g[i]=_g[i]=1;
56     }
57     for(int i=1;i<=N;i++)
58     for(int j=1;1ll*i*j<=N;j++){
59         if(!mu[j]) continue;
60         if(mu[j]==1) g[i*j]=1ll*g[i*j]*f[i]%mo;
61         else _g[i*j]=1ll*_g[i*j]*f[i]%mo;
62     }
63     for(int i=1;i<=N;i++)
64         g[i]=1ll*g[i]*inv(_g[i])%mo;
65     for(int i=2;i<=N;i++)
66         g[i]=1ll*g[i]*g[i-1]%mo;
67 }
68 int solve(int n,int m)
69 {
70     int re=1,last=0;
71     for(int t=1;t<=n;t=last+1){
72         last=min(n/(n/t),m/(m/t));
73         re=1ll*re*qkpow(1ll*g[last]*inv(g[t-1])%mo,1ll*(n/t)*(m/t))%mo;
74     }
75     return re;
76 }
77 int main()
78 {
79     scanf("%d",&T);
80     for(int i=1;i<=T;i++){
81         scanf("%d%d",&Q[i].n,&Q[i].m);
82         if(Q[i].n>Q[i].m) swap(Q[i].n,Q[i].m);
83         N=max(N,Q[i].n);
84     }
85     ready();
86     for(int i=1;i<=T;i++)
87         printf("%d\n",solve(Q[i].n,Q[i].m));
88     return 0;
89 }

 

posted @ 2018-03-26 18:57  KKKorange  阅读(94)  评论(0编辑  收藏