bzoj4816: [Sdoi2017]数字表格

链接

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

题解

　　大爷就是大爷…我就是蒻，我在这么蒻下去也许就连$D$都进不了了
　　这道题当时我做不出来也很正常，我从来没化简过带$\prod$的式子…而且这题细节也不少
　　好了开始写题解
　　

$$ans=\prod_{i=1}^n\prod_{j=1}^mf[gcd(i,j)]\\ =\prod_{d=1}^nf[d]^{\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=1]}\\　　=\prod_{d=1}^nf[d]^{\sum_{d|x}^n\mu(\frac{x}{d})\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{x}\rfloor}$$
　　对指数的求和与直接做乘积是一样的
$$ans=\prod_{d=1}^n\prod_{d|x}^nf[d]^{\mu(\frac{x}{d})\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{x}\rfloor}\\ =\prod_{x=1}^n\prod_{d|x}f[d]^{\mu(\frac{x}{d})\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{x}\rfloor}$$
　　考场上我就卡在这里了，下面继续化简
$$ans=\prod_{x=1}^n\left (\prod_{d|x}f[d]^{\mu(\frac{x}{d})}\right )^{\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{x}\rfloor}$$
　　这一步太神辣
　　 其实不是这一步太神是我对指数运算太不熟练了
　　然后令 $$g(n)=\prod_{x=1}^n\prod_{d|x}f[d]^{\mu(\frac{x}{d})}$$
　　这时
$$ans=\prod_{x=1}^ng(x)^{\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{x}\rfloor}$$
　　分块做到 $O(nlogn+T\sqrt N)$

代码

//莫比乌斯反演、画柿子
#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;
}