[BZOJ4816/Luogu3704][Sdoi2017]数字表格(莫比乌斯反演,数论分块)

发布于 2018-01-21  111 次阅读


本文章由SYCstudio或本站其它作者所有,请勿擅自转载

本文链接地址:[BZOJ4816/Luogu3704][Sdoi2017]数字表格(莫比乌斯反演,数论分块)

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取模。

Http

BZOJ
Luogu

Tag

莫比乌斯反演,数论分块

题目大意

求\(\prod_{i=1}^{n} \prod_{j=1}^{m} fab(gcd(i,j))\)

解决思路

依照惯例,提取出\(d=gcd(i,j)\),则
\[Ans=\prod_{d=1}^{n} fab(d)^{\sum_{i=1}^{n} \sum{j=1}^{m} [gcd(i,j)==1]}\]
\(fab(d)\)的指数可以根据莫比乌斯反演变成
\[Ans=\prod_{d=1}^{n} fab(d)^{\sum_{i=1}^{n/d} \mu(i) \lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor }\]
令\(T=id\),把\(T\)提出来,得到
\[ Ans=\prod_{T=1}^{n} [\prod_{i|T} fab(\frac{T}{i})^{\mu(i)}]^{\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor} \]
里面这个东西并不是积性函数,怎么办呢?
考虑到\(n\)只有\(10^6\),可以\(O(n*log_{2}n)\)暴力算,提前处理好斐波那契数列和其逆元即可。
那么设\(H(x)=\prod_{i|x} f(\frac{x}{i})^{\mu(i)}\),求出它的前缀积和前缀积的逆元即可数论分块了

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

#define ll long long
#define mem(Arr,x) memset(Arr,x,sizeof(Arr))

const int maxNum=1000010;
const int Mod=1e9+7;
const int inf=2147483647;

int Mu[maxNum];
int pricnt=0,Prime[maxNum/100];
bool notprime[maxNum];
int Fab[maxNum],InvFab[maxNum];
int H[maxNum],InvH[maxNum];

void Init();//线性筛求mu,以及预处理Fab,H及其逆元
int QPow(ll x,ll cnt,ll Mod);//快速幂
int Inv(ll x,ll Mod);//求逆元

int main()
{
    Init();
    int T;scanf("%d",&T);
    while (T--)
    {
        ll n,m;scanf("%lld%lld",&n,&m);if (n>m) swap(n,m);
        ll ans=1;
        for (ll i=1,last;i<=n;i=last+1)
        {
            last=min(n/(n/i),m/(m/i));
            ans=ans*QPow((ll)H[last]*(ll)InvH[i-1]%Mod,(ll)(n/i)*(ll)(m/i)%(Mod-1),Mod)%Mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}

void Init()
{
    notprime[1]=1;Mu[1]=1;
    for (int i=2;i<maxNum;i++)//线性筛出mu
    {
        if (notprime[i]==0) Mu[i]=-1,Prime[++pricnt]=i;
        for (int j=1;(j<=pricnt)&&((ll)i*(ll)Prime[j]<maxNum);j++)
        {
            notprime[i*Prime[j]]=1;
            if (i%Prime[j]==0) break;
            Mu[i*Prime[j]]=-Mu[i];
        }
    }
    Fab[0]=0;Fab[1]=1;H[0]=InvH[0]=1;
    //递推出fab
    for (int i=2;i<maxNum;i++) Fab[i]=(Fab[i-1]+Fab[i-2])%Mod;
    //求出fab的逆元
    for (int i=1;i<maxNum;i++) InvFab[i]=Inv(Fab[i],Mod),H[i]=1;
    for (int i=1;i<maxNum;i++)//求H
        if (Mu[i]!=0)
            for (int j=i;j<maxNum;j=j+i)
                H[j]=(ll)H[j]*(ll)((Mu[i]==1)?(Fab[j/i]):(InvFab[j/i]))%Mod;
    //求H的逆元
    for (int i=1;i<maxNum;i++) InvH[i]=Inv(H[i],Mod);
    //求前缀积
    for (int i=2;i<maxNum;i++) H[i]=(ll)H[i]*(ll)H[i-1]%Mod,InvH[i]=(ll)InvH[i]*(ll)InvH[i-1]%Mod;
    return;
}

int QPow(ll x,ll cnt,ll Mod)
{
    ll ret=1;
    while (cnt)
    {
        if (cnt&1) ret=ret*x%Mod;
        x=x*x%Mod;
        cnt=cnt>>1;
    }
    return ret%Mod;
}

int Inv(ll x,ll Mod)
{
    return QPow(x,Mod-2,Mod)%Mod;
}

本文章由SYCstudio或本站其它作者所有,请勿擅自转载

本文链接地址:[BZOJ4816/Luogu3704][Sdoi2017]数字表格(莫比乌斯反演,数论分块)


HNCJ OIer