并不对劲的bzoj4816:loj2000:p3704[SDOI2017]数字表格

题目大意

有函数(f(x))，(f(0)=0,f(1)=1,f(x)=f(x-1)+f(x-2))
(t)((tleq1000))组询问，每次给定(n,m)((n,mleq10^6))，求：

[prod_{i=1}^{n}prod_{j=1}^{m}f(gcd(i,j)) ]

题解

这个人（点这里）讲得很清楚(color{white}{ ext{shing太强了}})
假设(n)是(n,m)中较小的那个，讲一些在得出答案(=prod_{d=1}^{n}f(d)^{sum_{k=1}^{lfloor frac{n}{d} floor}{mu(k)lfloor frac{n}{dk} floorlfloor frac{m}{dk} floor}})之后的事
设(p=dk)则有答案(=prod_{d=1}^{n}f(d)^{sum_{dmid p}^{n}{mu(frac{p}{d})lfloor frac{n}{p} floorlfloor frac{m}{p} floor}})
把幂次上的(sum)拿下来变成(prod)，答案(=prod_{d=1}^{n}prod_{dmid p}^{n}f(d)^{{mu(frac{p}{d})lfloor frac{n}{p} floorlfloor frac{m}{p} floor}}=prod_{d=1}^{n}(prod_{dmid p}^{n}f(d)^{mu(frac{p}{d})} )^{lfloor frac{n}{p} floorlfloor frac{m}{p} floor})
括号里的部分(g(d)=prod_{dmid p}^{n}f(d)^{mu(frac{p}{d})})与(n,m)无关可以预处理，将每个(f(a)^{mu(b)})乘到(g(a*b))，预处理的时间复杂度为(frac{n}{1}+frac{n}{2}+frac{n}{3}+...+frac{n}{n})，大约是(n*log n)
接下来数论分块计算(prod_{d=1}^{n}g(d)^{lfloor frac{n}{p} floorlfloor frac{m}{p} floor})就行了，这部分的时间复杂度是(Theta(T(sqrt n +sqrt m)))

代码

#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<iomanip>
#include<iostream>
#include<map>
#include<queue>
#include<set>
#include<stack>
#include<vector>
#define rep(i,x,y) for(register int i=(x);i<=(y);++i)
#define dwn(i,x,y) for(register int i=(x);i>=(y);--i)
#define maxn 1000010
#define lim 1000000 
#define LL long long
using namespace std;
int read()
{
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)&&ch!='-')ch=getchar();
	if(ch=='-')f=-1,ch=getchar();
	while(isdigit(ch))x=(x<<1)+(x<<3)+ch-'0',ch=getchar();
	return x*f;
}
void write(int x)
{
	if(x==0){putchar('0'),putchar('
');return;}
	int f=0;char ch[20];
	if(x<0)putchar('-'),x=-x;
	while(x)ch[++f]=x%10+'0',x/=10;
	while(f)putchar(ch[f--]);
	putchar('
');
	return;
}
const LL mod=1e9+7;
int n,m,t,f[maxn],pf[maxn],g[maxn],pg[maxn],p[maxn],no[maxn],cnt,mu[maxn];
int mul(int x,LL y){int res=1;while(y){if(y&1ll)res=(LL)res*(LL)x%mod;x=(LL)x*(LL)x%mod,y>>=1;}return res;}
int main()
{
	no[1]=mu[1]=1;
	rep(i,2,lim)
	{
		if(!no[i])p[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*p[j]<=lim;j++)
		{
			no[i*p[j]]=1;
			if(i%p[j]==0){mu[i*p[j]]=0;break;}
			else mu[i*p[j]]=-mu[i];
		}
	}
	f[1]=pf[1]=g[1]=1;
	rep(i,2,lim)f[i]=(f[i-1]+f[i-2])%mod,pf[i]=mul(f[i],mod-2),g[i]=1;
	rep(i,1,lim)
	{
		for(int j=i;j<=lim;j+=i)
		{
			if(mu[j/i]==1)g[j]=(LL)g[j]*(LL)f[i]%mod;
			else if(mu[j/i]==-1)g[j]=(LL)g[j]*(LL)pf[i]%mod;
		}
	}g[0]=pg[0]=1;
	rep(i,1,lim)g[i]=(LL)g[i-1]*(LL)g[i]%mod,pg[i]=mul(g[i],mod-2);
	t=read();
	while(t--)
	{
		n=read(),m=read();
		if(n>m)swap(n,m);int ans=1;
		for(int l=1,r=0;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans=(LL)ans*mul((LL)g[r]*(LL)pg[l-1]%mod,(LL)(n/l)*(LL)(m/l))%mod;
		}
		write(ans);
	}
	return 0;
}