[SDOI2017]数字表格莫比乌斯+数论分块

Posted 2020-11-18 linda-fcj

一句话题意：

求:

$N=min(n,m)$
$prod_{d=1}^{N}prod_{i=1,j=1}^{n,m}f[d]*[gcd(i,j)=d]$

 

把$f[d]$提出来：

 

$=prod{f[d]^{sum_{i=1,j=1}^{n,m}[gcd(i,j)=d]}}$
$=prod{f[d]^{sum_{i=1,j=1}^{lfloor{frac{n}{d}} floor,lfloor{frac{m}{d}} floor}[gcd(i,j)=1]}}$

 

看指数：

 

$sum_{i=1,j=1}^{lfloor{frac{n}{d}} floor,lfloor{frac{m}{d}} floor}[gcd(i,j)=1]$
$=sum_{i=1}^{lfloor{frac{N}{d}} floor}mu(i)*lfloor{frac{n}{i*d}} floor*lfloor{frac{m}{i*d}} floor$

 

所以答案是：

 

  $prod{f[d]^{sum_{i=1}^{lfloor{frac{N}{d}} floor}mu(i)*lfloor{frac{n}{i*d}} floor*lfloor{frac{m}{i*d}} floor}}$

 

令$T=i*d$

 

$ prod_{T=1}^{n}(prod_{d|T}f[d]^{mu(T/d)})^{[n/T][m/T]}$

 

 把括号内部看成一个整体，对$[n/T][m/T]$数论分块，把括号内部的用前缀积预处理（注意预处理前缀积逆元）

预处理$O(sqrt{n})$，回答$O(T*sqrt{n})$，只有筛$mu$是$O(n)$，细节很多

CODE:

 

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxn=1e6+11;
const int mod=1e9+7;
long long n,m;
int TT;
int miu[maxn],p[maxn];
long long F[maxn],f[maxn],g[maxn];
bool mark[maxn];

inline long long Pow(long long x,int b)
{
	long long res=1;
	for(;b;b>>=1)
	{
		if(b&1) res=res*x%mod;
		x=x*x%mod;
	}
	return res;
}

inline void get_miu()
{
	miu[1]=1; g[1]=F[1]=F[0]=f[1]=1;
	for(int i=2;i<=maxn-10;i++)
	{
		f[i]=(f[i-1]+f[i-2])%mod;
		g[i]=Pow(f[i],mod-2);
		
		if(!mark[i])
		{
			p[++p[0]]=i;
			miu[i]=-1;
		}
		for(int j=1;j<=p[0] && p[j]*i<=maxn-10;j++)
		{
			mark[p[j]*i]=1;
			if(i%p[j]==0)
			{
				miu[i*p[j]]=0;
				break;
			}
			else miu[i*p[j]]=-miu[i];
		}
	}
}

inline void pre()
{
	get_miu();
	for(int i=0;i<=maxn;i++) F[i]=1;
    for(int bei=1;bei<=maxn-10;bei++)
    {
        if(miu[bei]==0) continue;
        for(int d=1;bei*d<=maxn-10;d++)
        {
            long long T=bei*d;
            F[T]=F[T]*(miu[bei]==1?f[d]:g[d])%mod;
        }
    }
	for(int i=2;i<=maxn-10;i++) F[i]=F[i]*F[i-1]%mod;
}

int main()
{
	scanf("%d",&TT);
	pre();
	long long ans=1;
	while(TT--)
	{
		scanf("%d%d",&n,&m);
		int N=min(n,m),nxt=0;
		ans=1;
		for(int T=1;T<=N;T=nxt+1)	
		{
			nxt=min(n/(n/T),m/(m/T));
			long long di=F[nxt]%mod*Pow(F[T-1],mod-2)%mod;
			di=Pow(di,(m/T)*(n/T)%(mod-1))%mod;//euler
			ans=(ans*di)%mod;
//			printf("%d**
",di);
		}
		printf("%lld
",ans);
	}
	return 0;
}