LG P6788 「EZEC-3」四月樱花

Description

在樱花盛开的四月，Muxii 望着满天飘落的樱花，向身旁的 ZZH 问道：

“究竟有多少朵樱花在这个四月飘落？”

ZZH 答道：“樱花飘落的朵数  $s$与时间 $t$ 有如下关系：

$s=prod_{x=1}^t prod_{y|x} frac{y^{d(y)}}{prod_{z|y}(z+1)^2}$

其中 $d(y)$ 表示 $y$ 的约数个数。”

但作为一个文科生萌新，Muxii 显然无法清楚地知道具体的数目，因此他只好继续向 ZZH 询问这个问题的答案。

由于数量可能很大，所以你只需要替 ZZH 告诉 Muxii 他所需要的答案对 $p$ 取模的结果就好了。

Solution

题中所求为

$$ans=left(prod_{x=1}^nprod_{y|x}frac{y^{d(y)}}{prod_{z|y}(z+1)^2} ight)mod p$$

大力操作式子：因为

$$y^{d(y)}=prod_{z|y}y=prod_{z|y}zcdotfrac{y}{z}=prod_{z|y}z^2$$

和

$$sum_{z|y,y|n}=d(frac{n}{z})$$

所以

$$s=prod_{x=1}^nprod_{y|x}frac{y^{d(y)}}{prod_{z|y}(z+1)^2}=prod_{x=1}^nprod_{y|x}prod_{z|y}left(frac{z}{z+1} ight)^2=prod_{z=1}^nleft(frac{z}{z+1} ight)^{2sumd(lfloorfrac{n}{z} floor)}$$

其中$sumd(n)=sum_{i=1}^n d(i)$

由于要求$d$的前缀和，所以有杜教筛

设$f=d=1*1，g=mu$，则$f*g=1*1*mu=1$，可以使用杜教筛，现在需要求出$mu$的前缀和

再次使用杜教筛即可

求原式的值可以用整除分块做

时间复杂度$Theta(n^{frac 23}+sqrt{n}log n)$

#include<iostream>
#include<cstdio>
using namespace std;
long long n,mod,pmod,prime[6000005],mu[6000005],minn[6000005],d[6000005],tot,summu[605],sumd[605],ans=1;
bool isprime[6000005],vst[600];
inline long long read()
{
    long long w=0,f=1;
    char ch=0;
    while(ch<'0'||ch>'9')
    {
        if(ch=='-')
            f=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9')
    {
        w=(w<<1)+(w<<3)+ch-'0';
        ch=getchar();
    }
    return w*f;
}
long long Mu(long long x)
{
    return x<=5000000?mu[x]:summu[n/x];
}
long long D(long long x)
{
    return x<=5000000?d[x]:sumd[n/x];
}
void djs(long long N)
{
    if(N<=5000000||vst[n/N])
    {
        return;
    }
    vst[n/N]=true;
    for(long long l=2;l<=N;)
    {
        long long r=N/(N/l);
        djs(N/l);
        (summu[n/N]+=(pmod-(r-l+1)*Mu(N/l)%pmod)%pmod)%=pmod;
        l=r+1;
    }
    (summu[n/N]+=1)%=pmod;
    for(long long l=2;l<=N;)
    {
        long long r=N/(N/l);
        (sumd[n/N]+=pmod-((Mu(r)-Mu(l-1)+pmod)%pmod)*D(N/l)%pmod)%=pmod;
        l=r+1;
    }
    (sumd[n/N]+=N)%=pmod;
}
long long ksm(long long a,long long p)
{
    long long ret=1;
    while(p)
    {
        if(p&1)
        {
            (ret*=a)%=mod;
        }
        (a*=a)%=mod;
        p>>=1;
    }
    return ret;
}
int main()
{
    n=read();
    mod=read();
    pmod=mod-1;
    mu[1]=d[1]=1;    
    for(long long i=2;i<=5000000;i++)
    {
        if(!isprime[i])
        {
            prime[++tot]=i;
            mu[i]=-1;
            d[i]=2;
            minn[i]=1;
        }
        for(long long j=1;j<=tot&&i*prime[j]<=5000000;j++)
        {
            isprime[i*prime[j]]=true;
            if(!(i%prime[j]))
            {
                minn[i*prime[j]]=minn[i]+1;
                d[i*prime[j]]=d[i]/(minn[i]+1)*(minn[i*prime[j]]+1);
                mu[i*prime[j]]=0;
                break;
            }
            minn[i*prime[j]]=1;
            d[i*prime[j]]=d[i]*d[prime[j]];
            mu[i*prime[j]]=mu[i]*-1;
        }
    }
    for(long long i=2;i<=5000000;i++)
    {
        (d[i]+=d[i-1])%=pmod;
        ((mu[i]+=mu[i-1])+=pmod)%=pmod;
    }
    djs(n);
    for(long long l=1;l<=n;)
    {
        long long r=n/(n/l);
        (ans*=ksm(l*ksm(r+1,mod-2)%mod,D(n/l)))%=mod;
        l=r+1;
    }
    (ans*=ans)%=mod;
    printf("%lld
",ans);
    return 0;
}