这题需要用拉格朗日插值法求前n项的 ik 的和。
然后用min25筛一下质数的贡献。
有点卡常，能预处理的都预处理。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<string>
#include<queue>
#define ll long long
#define llu unsigned ll
using namespace std;
const int maxn=1001000;
ll mod;
ll pre[maxn],suf[maxn],fac[maxn],inv[maxn],a[maxn];
ll y[maxn];
ll inv1;

ll p[maxn],id1[maxn],id2[maxn],pc[maxn];
ll g[maxn],sum[maxn],wi[maxn];
ll cnt,sq,n,k;
bool ha[maxn];

ll mypow(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}

void init(ll k)
{
    int m=k+2;
    inv1=mypow(mod-1,mod-2);
    fac[0]=1,y[0]=0;
    for(int i=1;i<=m;i++)
        y[i]=(y[i-1]+mypow(i,k))%mod,fac[i]=fac[i-1]*i%mod;
    inv[m]=mypow(fac[m],mod-2);
    for(int i=m-1;i>=0;i--)
        inv[i]=inv[i+1]*(i+1)%mod;
}

ll lg(ll n,ll k)
{
    int m=k+2;
    n%=mod;
    pre[0]=1;
    for(int i=1;i<=m;i++)
        pre[i]=pre[i-1]*(n-i+mod)%mod;

    suf[m+1]=1;
    suf[m]=suf[m+1]*(n-m+mod)%mod;
    for(int i=m-1;i>=0;i--)
        suf[i]=suf[i+1]*(n-i+mod)%mod;

    ll ans=0;
    for(int i=1;i<=m;i++)
        ans=(ans+y[i]*pre[i-1]%mod*suf[i+1]%mod*inv[i-1]%mod*inv[m-i]%mod*((m-i)&1?inv1:1)%mod)%mod;
    return ans;
}


void Prime(ll n)
{
   ha[1]=true;
   for(int i=2;i<=n;i++)
   {
       if(!ha[i])
       {
           p[++cnt]=i;
           sum[cnt]=(sum[cnt-1]+mypow(i,k))%mod;
           pc[cnt]=mypow(i,k);
       }
       for(int j=1;j<=cnt&&p[j]*i<=n;j++)
       {
           ha[i*p[j]]=true;
           if(i%p[j]==0) break;
       }
   }
}

int getid(ll m)
{
    if(m<=sq) return id1[m];
    else return id2[n/m];
}

void min25(ll n,ll k)
{
    sq=sqrt(n*1.0);
    Prime(sq);
    init(k);

    int tot=0;
    for(ll l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        wi[++tot]=n/l;
        if(wi[tot]<=sq) id1[wi[tot]]=tot;
        else id2[n/wi[tot]]=tot;
        g[tot]=(lg(wi[tot],k)-1+mod)%mod;
    }

    for(int j=1;j<=cnt;j++)
    {
        for(int i=1;i<=tot&&p[j]*p[j]<=wi[i];i++)
        {
            int cm=getid(wi[i]/p[j]);
            g[i]=((g[i]-pc[j]*(g[cm]-sum[j-1])%mod)%mod+mod)%mod;
        }
    }
}

int main(void)
{
    scanf("%lld%lld%lld",&n,&k,&mod);
    min25(n,k);
    ll ans=0;
    //这里其实还可以分块优化。
    for(ll i=1;i<=sq;i++)
        ans=(ans+i*i%mod*g[getid(n/i)])%mod;
    printf("%lld\n",ans);
    return 0;
}