$bzoj4816-SDOI2016$ 數字表格 莫比烏斯反演
題面描述
\(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取模。
輸入格式
有多組測試數據。
第一個一個數T,表示數據組數。
接下來T行,每行兩個數n,m
T<=1000,1<=n,m<=10^6
輸出格式
- 輸出T行,第i行的數是第i組數據的結果
題解
根據題意,我們能列出初始計算公式\(ans=\prod_{i=1}^n\prod_{j=1}^mf_{gcd}(i,j)\),這個式子在\(n,m\leq10^6?\)時很明顯不可做。
\[ \prod_{i=1}^n\prod_{j=1}^mf_{gcd(i,j)}\ (n<m)\\prod_{d=1}^nf_{d}^{\sum_{x=1}^{\frac{n}{d}}\sum_{y=1}^{\frac{m}{d}}[gcd(x,y)=1]}\\prod_{d=1}^nf_{d}^{\sum_{x=1}^{\frac{n}{d}}\sum_{y=1}^{\frac{m}{d}}\sum_{d'|gcd(x,y)}\mu(d')}\\prod_{d=1}^nf_{d}^{\sum_{d'=1}^{\frac{n}{d}}\lfloor\frac{n}{d*d'}\rfloor\lfloor\frac{m}{d*d'}\rfloor\mu(d')}\\prod_{p=1}^n\prod_{k|p}f_k^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor\mu(\frac{p}{k})}\\prod_{p=1}^n\prod_{k|p}(f_k^{\mu(\frac{p}{k})})^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor}\\prod_{p=1}^npre_p^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor} \]
得到這個式子後,利用\(\lfloor\frac{n}{p}\rfloor\),\(\lfloor\frac{m}{p}\rfloor\)共\(o(n)\)個取值,分塊+前綴積處理即可,時間復雜度\(o(10^6log.mod+Tnlog. mod)\),其中log mod為求逆元的快速冪復雜度。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; typedef long long ll; const int mod=1e9+7; const int MAXN=1e6+6; int T,n,m; int pri[MAXN],mu[MAXN]; bool use[MAXN]; int fb[MAXN],f[MAXN][3]; int pre[MAXN]; int fm[MAXN],fr[MAXN]; int mod_pow(int x,int n){ int ret=1; while (n){ if (n&1) ret=(ll)ret*x%mod; x=(ll)x*x%mod; n>>=1; } return ret; } int main(){ use[1]=0; mu[1]=1; for (int i=2;i<=1e6;i++){ if (!use[i]) pri[++pri[0]]=i,mu[i]=-1; for (int j=1;j<=pri[0]&&i*pri[j]<=1e6;j++){ use[i*pri[j]]=1; if (i%pri[j]!=0) mu[i*pri[j]]=-mu[i]; else break; } } fb[1]=1; fb[2]=1; for (int i=3;i<=1e6;i++) fb[i]=(ll)(fb[i-1]+fb[i-2])%mod; for (int i=1;i<=1e6;i++){ f[i][0]=mod_pow(fb[i],mod-2); f[i][1]=1; f[i][2]=fb[i]; } for (int i=1;i<=1e6;i++) pre[i]=1; for (int i=1;i<=1e6;i++){ for (int j=i;j<=1e6;j+=i){ pre[j]=(ll)pre[j]*f[i][mu[j/i]+1]%mod; } } // for (int i=1;i<=10;i++) printf("%d\n",pre[i]); // cout<<pre[(int)1e6]<<endl; fm[0]=1; for (int i=1;i<=1e6;i++) fm[i]=(ll)fm[i-1]*pre[i]%mod; fr[(int)1e6]=mod_pow(fm[(int)1e6],mod-2); // cout<<fr[(int)1e6]<<endl; for (int i=(int)1e6-1;i>=0;i--) fr[i]=(ll)fr[i+1]*pre[i+1]%mod; // for (int i=1;i<=10;i++) printf("%d %d\n",fm[i],fr[i]); scanf("%d",&T); while (T--){ scanf("%d%d",&n,&m); if (n>m) swap(n,m); int j=0,ans=1; for (int i=1;i<=n;i=j+1){ j=min(n/(n/i),m/(m/i)); ans=(ll)ans*mod_pow((ll)fr[i-1]*fm[j]%mod,(ll)(n/i)*(m/i)%(mod-1))%mod; // cout<<i<<" "<<j<<" "<<ans<<endl; } printf("%d\n",ans); } return 0; }
$bzoj4816-SDOI2016$ 數字表格 莫比烏斯反演