• bzoj4816: [Sdoi2017]数字表格


    Description

    Doris刚刚学习了(fibonacci)数列。用(f[i])表示数列的第(i)项,那么

    [egin{align}f[0]&=0\f[1]&=1\f[n]&=f[n-1]+f[n-2],n>=2end{align} ]

    Doris用老师的超级计算机生成了一个(n imes m)的表格,第(i)行第(j)列的格子中的数是(f[mathrm{gcd}(i,j)]),其中(mathrm{gcd}(i,j))表示(i),(j)的最大公约数。Doris的表格中共有(n imes m)个数,她想知道这些数的乘积是多少。答案对(10^9+7)取模。

    Input

    有多组测试数据。第一个一个数(T),表示数据组数。

    接下来(T)行,每行两个数(n,m)

    (T leq 1000),(1 leq n),(m leq 10^6)

    Output

    输出(T)行,第(i)行的数是第(i)组数据的结果

    Sample Input

    3
    2 3
    4 5
    6 7
    

    Sample Output

    1
    6
    960
    

    题解

    要求的东西是

    [ans=prod_{i=1}^nprod_{j=1}^mf[mathrm{gcd}(i,j)] ]

    假设(n<m),根据常见套路,枚举约数

    [ans=prod_{d=1}^nf(d)^{sum_{i=1}^nsum_{j=1}^m[mathrm{gcd}(i,j)==d]} ]

    把指数拿出来化简

    [egin{align} &sum_{i=1}^nsum_{j=1}^m[mathrm{gcd}(i,j)==d]\ =&sum_{i=1}^{left lfloor frac{n}{d} ight floor}sum_{j=1}^{left lfloor frac{m}{d} ight floor}[mathrm{gcd}(i,j)==1]\ =&sum_{i=1}^{left lfloor frac{n}{d} ight floor}sum_{j=1}^{left lfloor frac{m}{d} ight floor}sum_{k|mathrm{gcd}(i,j)}mu(k)\ =&sum_{k=1}^{left lfloor frac{n}{d} ight floor}mu(k)left lfloor frac{n}{kd} ight floorleft lfloor frac{m}{kd} ight floor end{align} ]

    (T=kd),代入上式

    [sum_{k=1}^{left lfloor frac{n}{d} ight floor}mu(k)left lfloor frac{n}{kd} ight floorleft lfloor frac{m}{kd} ight floor=sum_{T=1}^nmu(frac{T}{d})left lfloor frac{n}{T} ight floorleft lfloor frac{m}{T} ight floor\ ]

    (T)放到外面枚举

    [ans=prod_{T=1}^nprod_{d|T}f(d)^{mu(frac{T}{d})left lfloor frac{n}{T} ight floorleft lfloor frac{m}{T} ight floor} ]

    令函数(g(x)=prod_{d|x}f(d)^{mu(frac{x}{d})})

    [ans=prod_{T=1}^ng(T)^{left lfloor frac{n}{T} ight floorleft lfloor frac{m}{T} ight floor} ]

    只要我们求出了(g(x))的前缀积以及对应的逆元就可以整除分块了。

    只要预处理出莫比乌斯函数和(f)以及(f)的逆元,(g(x))就可以通过枚举每个数的筛它的倍数求出,复杂度是调和级数。

    问题就是如何求所有(f)以及对应的逆元(invf)

    这里我们需要一个辅助的东西(mulf),表示(f)的前缀积,即(mulf(x)=prod_{i=1}^x f(i)).

    假设我们知道了(mulf(i+1))的逆元(invm(i+1)),那么有

    [egin{align} invm(i+1) &equiv invm(i) imes invf(i+1) (mathrm{mod} P)\ invm(i)&equiv invm(i+1) imes f(i+1)(mathrm{mod} P) end{align} ]

    (invf(i)equiv invm(i)*mulf(i-1)(mathrm{mod} P)),这样我们就可以求出(f)以及(invf)了。

    用同样的方法,我们也可以求出(g(x))的前缀积以及对应的逆元。

    代码

    #include<bits/stdc++.h>
    #define MAXN 1000010
    #define P 1000000007
    #define LL long long
    namespace IO{
    	char buf[1<<15],*fs,*ft;
    	inline char gc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
    	inline int qr(){
    		int x=0,rev=0,ch=gc();
    		while(ch<'0'||ch>'9'){if(ch=='-')rev=1;ch=gc();}
    		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=gc();}
    		return rev?-x:x;}
    }using namespace IO;
    using namespace std;
    LL f[MAXN],invm[MAXN],invf[MAXN],g[MAXN],invg[MAXN],mulf[MAXN],mulg[MAXN],ans;
    int p[MAXN],cnt,u[MAXN],T,N,M;
    bool np[MAXN]; 
    inline LL qp(LL x,LL y){
    	LL ret=1;
    	while(y){
    		if(y&1)ret=ret*x%P;
    		x=x*x%P;y>>=1; 
    	}
    	return ret;
    }
    int main(){
    	#ifndef ONLINE_JUDGE
    	freopen("product.in","r",stdin);
    	freopen("product.out","w",stdout);
    	#endif
    	f[0]=0;
    	f[1]=u[1]=np[1]=g[1]=mulf[1]=mulg[1]=invm[0]=1;
    	for(int i=2;i<=1000000;i++){
    		f[i]=(f[i-1]+f[i-2])%P;
    		mulf[i]=mulf[i-1]*f[i]%P;
    		g[i]=1;
    		if(!np[i])p[++cnt]=i,u[i]=-1;
    		for(int j=1;j<=cnt;j++){
    			int t=p[j]*i;
    			if(t>1000000)break;
    			np[t]=1;
    			if(i%p[j]==0){u[t]=0;break;}
    			u[t]=-u[i];
    		}
    	}
    	invm[1000000]=qp(mulf[1000000],P-2);
    	invf[1000000]=qp(f[1000000],P-2);
    	for(int i=999999;i>=1;i--){
    		invm[i]=invm[i+1]*f[i+1]%P;
    		invf[i]=mulf[i-1]*invm[i]%P;
    	} 
    	for(int i=2;i<=1000000;i++){
    		for(int j=1,k;(k=j*i)<=1000000;j++){
    			if(u[j]==1)g[k]=g[k]*f[i]%P;
    			else if(u[j]==-1)g[k]=g[k]*invf[i]%P;
    		}
    		mulg[i]=mulg[i-1]*g[i]%P;
    	}
    	invg[1000000]=qp(mulg[1000000],P-2);
    	for(int i=999999;i>=0;i--)invg[i]=invg[i+1]*g[i+1]%P;
    	T=qr();
    	while(T--){
    		N=qr();M=qr();ans=1;
    		if(N>M)swap(N,M);
    		for(int i=1,j;i<=N;i=j+1){
    			j=min(N/(N/i),M/(M/i));
    			ans=ans*qp(mulg[j]*invg[i-1]%P,1LL*(N/i)*(M/i)%(P-1))%P;
    		}
    		printf("%lld
    ",ans);
    	}
    	return 0;
    }
    
  • 相关阅读:
    第十二周进度表
    第一个冲刺周期-第十天
    第一个冲刺周期-第九天
    团队作业—第二阶段06
    团队作业—第二阶段05
    团队作业—第二阶段04
    团队作业—第二阶段03
    团队作业—第二阶段02
    团队作业—第二阶段01
    对于风行小组第一阶段冲刺成果的概括
  • 原文地址:https://www.cnblogs.com/lrj998244353/p/8858100.html
Copyright © 2020-2023  润新知