2019 西安邀请赛 B. Product 杜教筛

最新推荐文章于 2020-09-11 23:06:32 发布

长沙橘子猫

最新推荐文章于 2020-09-11 23:06:32 发布

阅读量1.3k

点赞数

分类专栏：数论----杜教筛

本文链接：https://blog.csdn.net/ccsu_cat/article/details/90579297

版权

数论----杜教筛专栏收录该内容

3 篇文章 0 订阅

订阅专栏

B. Product

我们先求幂ans：

$ans=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}gcd(i,j)[k|gcd(i,j)]=\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)d(gcd(i,j))$

切换枚举次序：

$ans=\sum_{x=1}^{n}xd(x)\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=x]=\sum_{x=1}^{n}xd(x)\sum_{i=1}^{\frac{n}{x}}\sum_{j=1}^{\frac{n}{x}}[gcd(i,j)=1]$

将互质转化为欧拉函数：

$ans=\sum_{x=1}^{n}xd(x)(2\sum_{i=1}^{\frac{n}{x}}\varphi (i)-1)$

对于比较大的前缀和欧拉函数Sn，我们使用杜教筛：

$id=\varphi *I$

$S(n)=\sum_{i=1}^{n}i-\sum_{i=2}^{n}S\left \lfloor \frac{n}{i} \right \rfloor$

对于较大的前缀和 xd(x) ：

$\sum_{x=1}^{n}xd(x)=\sum_{x=1}^{n}\sum_{i|x}x=\sum_{i=1}^{n}\sum_{i|x}^{n}x=\sum_{i=1}^{n}i\sum_{x=1}^{\frac{n}{i}}x$

那么，这题就做完咯。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=5e6+5,N=5e6;
int vis[maxn],phi[maxn],pri[maxn/10],d[maxn],cnt,mod;
unordered_map<int,int>mp,mp2;
void add(int &x,int y)
{
    if(y<0)
        x+=y;
    else
        x=x-mod+y;
    if(x<0)
        x+=mod;
}
void init(int n)
{
    phi[1]=d[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
            pri[++cnt]=i,d[i]=2,phi[i]=i-1;
        for(int j=1;j<=cnt&&pri[j]*i<=n;j++)
        {
            vis[pri[j]*i]=1;
            if(i%pri[j])
                phi[i*pri[j]]=phi[i]*phi[pri[j]],d[i*pri[j]]=d[i]*2;
            else
            {
                phi[i*pri[j]]=phi[i]*pri[j];
                int sum=1,x=i;
                while(x%pri[j]==0)
                    x/=pri[j],sum++;
                d[i*pri[j]]=d[i]/sum*(sum+1);
                break;
            }
        }
    }
    for(int i=1;i<=n;i++)
    {
        add(phi[i],phi[i-1]);
        d[i]=1ll*d[i]*i%mod;
        add(d[i],d[i-1]);
    }
}
int ksm(ll x,int y,int M)
{
    ll res=1;
    while(y)
    {
        if(y&1)
            res=res*x%M;
        x=x*x%M;
        y/=2;
    }
    return res;
}
int Abs(int x)
{
    if(x<0)
        x+=mod;
    return x;
}
int calc(int l,int r)
{
    return 1ll*(l+r)*(r-l+1)/2%mod;
}
int work(int n)
{
    if(n<=N)
        return d[n];
    if(mp2[n])
        return mp2[n];
    int ans=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        add(ans,1ll*calc(l,r)*calc(1,n/l)%mod);
    }
    return mp2[n]=ans;
}
int dfs(int n)
{
    if(n<=N)
        return phi[n];
    if(mp[n])
        return mp[n];
    int ans=calc(1,n),res=0;
    for(int l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        add(res,1ll*(r-l+1)*dfs(n/l)%mod);
    }
    add(ans,-res);
    return mp[n]=ans;
}
int main()
{
    int n,m,p,ans=0;
    cin>>n>>m>>p;
    mod=p-1;
    init(N);
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        add(ans,1ll*(2ll*dfs(n/l)-1)%mod*Abs(work(r)-work(l-1))%mod);
    }
    cout<<ksm(m,ans,p);
}