1. 程式人生 > >$bzoj4816-SDOI2016$ 數字表格 莫比烏斯反演

$bzoj4816-SDOI2016$ 數字表格 莫比烏斯反演

cstring return ons -m 利用 初始 n-1 log out

  • 題面描述

    • \(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$ 數字表格 莫比烏斯反演