bzoj4816: [Sdoi2017]数字表格

链接

　http://www.lydsy.com/JudgeOnline/problem.php?id=4816

题解

　　大爷就是大爷…我就是蒻，我在这么蒻下去也许就连D都进不了了

　　这道题当时我做不出来也很正常，我从来没化简过带∏的式子…而且这题细节也不少

　　好了开始写题解

　　ans=∏i=1n∏j=1mf[gcd(i,j)]=∏d=1nf[d]∑ni=1∑mj=1[gcd(i,j)=1]　　=∏d=1nf[d]∑nd|xμ(xd)⌊nx⌋⌊mx⌋

　　对指数的求和与直接做乘积是一样的

ans=∏d=1n∏d|xnf[d]μ(xd)⌊nx⌋⌊mx⌋=∏x=1n∏d|xf[d]μ(xd)⌊nx⌋⌊mx⌋

　　考场上我就卡在这里了，下面继续化简

ans=∏x=1n⎛⎝∏d|xf[d]μ(xd)⎞⎠⌊nx⌋⌊mx⌋

　　这一步太神辣

　　其实不是这一步太神是我对指数运算太不熟练了

　　然后令g(n)=∏x=1n∏d|xf[d]μ(xd)=∏x=1n∏d|xf[d]μ(d)

　　这时

ans=∏x=1ng(x)⌊nx⌋⌊mx⌋

　　分块做到O(nlogn+TN−−√)

代码

//莫比乌斯反演、画柿子
#include <cstdio>
#include <algorithm>
#define maxn  1000005
#define mod 1000000007ll
#define ll long long
using namespace std;
ll N, f[maxn], _f[maxn], g[maxn], _g[maxn];
int mu[maxn], prime[maxn/10];
bool mark[maxn];
inline ll pow(ll a, ll b, ll p)
{
ll t, ans;
for(t=a,ans=1;b;b>>=1,t=t*t%p)if(b&1)ans=ans*t%p;
return ans;
}
void init()
{
ll i, j, t;
mu[1]=1;
for(i=2;i<maxn;i++)
{
if(!mark[i]){prime[++*prime]=i;mu[i]=-1;}
for(j=1;i*prime[j]<maxn;j++)
{
mark[i*prime[j]]=1;
if(i%prime[j]==0){mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=-mu[i];
}
}
f[1]=1;
for(i=2;i<maxn;i++)f[i]=(f[i-1]+f[i-2])%mod;
for(i=1;i<maxn;i++)_f[i]=pow(f[i],mod-2,mod);
for(i=0;i<maxn;i++)g[i]=1;
for(i=1;i<maxn;i++)for(j=i;j<maxn;j+=i)
{
if(mu[j/i]==1)t=f[i];
else if(mu[j/i]==-1)t=_f[i];
else t=1;
g[j]=g[j]*t%mod;
}
for(i=1;i<maxn;i++)g[i]=(g[i]*g[i-1])%mod;
for(i=1;i<maxn;i++)_g[i]=pow(g[i],mod-2,mod);
g[0]=_g[0]=1;
}
void solve(int n, int m)
{
int x, d, last;
ll ans=1;
if(n>m)swap(n,m);
for(x=1;x<=n;x=last+1)
{
last=min(n/(n/x),m/(m/x));
ans=ans*pow(g[last]*_g[x-1]%mod,(ll)(n/x)*(m/x)%(mod-1),mod)%mod;
}
printf("%lld\n",ans);
}
int main()
{
int T, n, m;
init();
for(scanf("%d",&T);T;T--)scanf("%d%d",&n,&m),solve(n,m);
return 0;
}