一道数学题

今天做了一道有趣的数学题。

给定n、k， 求$\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)^k$。 

我们来推一波公式，可以发现原式等于$\sum_{i=1}^ni^k\sum_{i|p}\mu(p/i)(n/p)^2$ 。(/为除号，下取整)这个莫比乌斯反演一下就好了。

设p=ig，那么原式为$\sum_{i=1}^ni^k\sum_g\mu(g)(n/i/g)^2$ 。

然后我们可以枚举n/i的值，这样的值只有根号个。那么这一段$i^k$ 之和可以插值，现在只要解决后面一段就可以了。

设$p(n)=\sum_{i=1}^n\mu(i)(n/i)^2$ ，对于每个n/i，后面一段就是p(n/i)。

根号算显然太慢了，我们可以发现$p(n)=-1+\sum_{i=1}^n2\varphi(i)$ 。

考虑p(n)相当于用莫比乌斯反演在算gcd(i,j)=1的对数，枚举较大的一个用欧拉函数算*2之后扣掉(1,1)即可。

这样我们只要枚举n/i+杜教筛即可。这样的复杂度是$O(n^{\frac{2}{3}})$ 的。

因为杜教筛的时候已经算出了所有n/i的前缀和了。

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <string>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>
#include <sstream>
#include <stack>
#include <iomanip>
using namespace std;
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
typedef double ld;
typedef vector<int> vi;
#define fi first
#define se second
#define fe first
#define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);}
#define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);}
#define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);}
#define es(x,e) (int e=fst[x];e;e=nxt[e])
#define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e])
#define VIZ {printf("digraph G{\n"); for(int i=1;i<=n;i++) for es(i,e) printf("%d->%d;\n",i,vb[e]); puts("}");}
#define VIZ2 {printf("graph G{\n"); for(int i=1;i<=n;i++) for es(i,e) if(vb[e]>=i)printf("%d--%d;\n",i,vb[e]); puts("}");}
#define SZ 666666
ll n,MOD=1e9+7; int k;
ll qp(ll a,ll b)
{
    ll ans=1;a%=MOD;
    while(b)
    {
        if(b&1) ans=ans*a%MOD;
        a=a*a%MOD; b>>=1;
    }
    return ans;
}
//拉格朗日大法好
int r[13],__[25],*ny=__+7;
ll f(ll s)
{
    if(s<=k+1) return r[s];
    s%=MOD; 
    ll tt=0;
    for(int i=0;i<=k+1;i++)
    {
        ll t=r[i];
        for(int j=0;j<=k+1;j++)
        {
            if(i==j) continue;
            t=t*(s-j)%MOD*ny[i-j]%MOD;
        }
        tt=(tt+t)%MOD;
    }
    return tt;
}
#define MM 5000000
int ps[MM+5],pn=0,mu[MM+5],eul[MM+5];
bool np[MM+5];
ll qzheul[MM+5];
void shai()
{
    np[1]=mu[1]=eul[1]=1;
    for(int i=2;i<=MM;i++)
    {
        if(!np[i])
        {mu[i]=-1; ps[++pn]=i; eul[i]=i-1;}
        for(int j=1;j<=pn&&i*ps[j]<=MM;j++)
        {
            np[i*ps[j]]=1;
            if(i%ps[j])
            {
                mu[i*ps[j]]=-mu[i];
                eul[i*ps[j]]=eul[i]*(ps[j]-1);
            }
            else
            {
                mu[i*ps[j]]=0;
                eul[i*ps[j]]=eul[i]*ps[j];
                break;
            }
        }
    }
    for(int i=1;i<=MM;i++)
    qzheul[i]=(qzheul[i-1]+eul[i])%MOD;
}
#define G 233333
ll p1[G],p2[G],rv2=qp(2,MOD-2);
ll &gm(ll x)
{
    if(x<G) return p1[x];
    return p2[n/x];
}
ll eulsum(ll x)
{
    if(x<=MM) return qzheul[x];
    ll &ps=gm(x);
    if(ps!=-1) return ps;
    ll ans,lst;
    if(x%MOD!=0)
        ans=(x%MOD)*((x%MOD)+1)%MOD*rv2%MOD;
    else
        ans=0;
    for(ll p=2;p<=x;p=lst+1)
    {
        lst=x/(x/p);
        ans-=((lst-p+1)%MOD)*eulsum(x/p)%MOD;
        ans%=MOD;
    }
    ans=(ans%MOD+MOD)%MOD;
    return ps=ans;
}
int main()
{
    memset(p1,-1,sizeof(p1));
    memset(p2,-1,sizeof(p2));
    shai();
    cin>>n>>k;
    for(int i=-k-1;i<=k+1;i++)
        ny[i]=qp(i,MOD-2);
    for(int i=1;i<=k+1;i++)
    {
        ll ans=1;
        for(int p=1;p<=k;p++)
            ans=ans*i%MOD;
        r[i]=(r[i-1]+ans)%MOD;
    }
    ll ls=0,ans=0;
    for(ll i=1;i<=n;i=n/(n/i)+1)
    {
        ll cur=f(n/(n/i)),ca=cur-ls;ca%=MOD;
        ll g=n/i,p=-1+eulsum(g)*2; p%=MOD;
        ans=ans+ca*p%MOD;
        ans%=MOD;
        ls=cur;
    }
    ans=(ans%MOD+MOD)%MOD;
    printf("%d\n",int(ans));
}