Luogu 4917 天守阁的地板(莫比乌斯反演+线性筛)

  既然已经学傻了,这个题当然是上反演辣。

  对于求积的式子,考虑把[gcd=1]放到指数上。一通套路后可以得到∏D∏d∏i∏j (ijd2)μ(d) (D=1~n,d|D,i,j=1~n/D)。

  冷静分析一下,由μ*1=e,后面一串ij相关的式子仅当D=1时有贡献。这一部分就非常好算了。而d对某个D的贡献,容易发现是d2μ(d)*(n/D)^2。设f(D)=∏dμ(d) (d|D),这个式子是可以线性筛的。(事实上从莫比乌斯函数的性质上看好像也很可以求,然而已经不会了)筛完之后就可以愉快的整除分块了。

  于是我们最后得到了一个不需要莫比乌斯函数的式子。复杂度O(n+t√nlogn)。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
#define P 19260817
#define N 1000010
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int T,n,fac[N],g[N],h[N],prime[N],cnt;
bool flag[N];
int ksm(int a,int k)
{
    if (k<0) k=1ll*(P-2)*(-k)%(P-1);
    k%=(P-1);
    int s=1;
    for (;k;k>>=1,a=1ll*a*a%P) if (k&1) s=1ll*s*a%P;
    return s;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("a.in","r",stdin);
    freopen("a.out","w",stdout);
    const char LL[]="%lld
";
#else
    const char LL[]="%I64d
";
#endif
    T=read();
    fac[0]=1;for (int i=1;i<=1000000;i++) fac[i]=1ll*fac[i-1]*i%P;
    for (int i=1;i<=1000000;i++) fac[i]=ksm(fac[i],i);
    flag[1]=1;g[0]=g[1]=1;
    for (int i=2;i<=N-10;i++)
    {
        if (!flag[i]) prime[++cnt]=i,g[i]=ksm(i,-1);
        for (int j=1;j<=cnt&&prime[j]*i<=N-10;j++)
        {
            flag[prime[j]*i]=1;
            if (i%prime[j]==0) {g[prime[j]*i]=g[i];break; }
            g[prime[j]*i]=1;
        }
    }
    for (int i=1;i<=N-10;i++) g[i]=1ll*g[i]*g[i-1]%P;
    while (T--)
    {
        int n=read(),ans=1ll*fac[n]*fac[n]%P;
        for (int i=1;i<=n;i++)
        {
            int t=n/(n/i);
            ans=1ll*ans*ksm(1ll*g[t]*ksm(g[i-1],-1)%P,2ll*(n/i)*(n/i)%(P-1))%P;
            i=t;
        }
        cout<<ans<<endl;
    }
    return 0;
}
原文地址:https://www.cnblogs.com/Gloid/p/9752158.html