【GDOI2018模拟7.12】C

Description

求

$$\sum_{i=1}^n\sum_{j=1}^n gcd(i,j)^k$$

Input

两个整数n,k

Output

答案需要对10^9+7取模。

Sample Input

2 2

Sample Output

7

Solution

双sigma+gcd=反演

$$ans=\sum_{i=1}^n\sum_{j=1}^n gcd(i,j)^k$$
枚举gcd=d
方便起见，除法默认下取整
$$ans=\sum_{d=1}^nd^k\sum_{i=1}^{n/d}\sum_{j=1}^{n/d} [gcd(i,j)=1]$$
后面的gcd就是互质的数的数量，也就是phi，因为数对有序，所以
$$ans=\sum_{d=1}^nd^k(\sum_{i=1}^{n/d}2*\phi(i))-1$$
即
$$ans=\sum_{d=1}^nd^k*2*(\sum_{i=1}^{n/d}\phi(i))-1$$
前面用自然数幂和，因为这题k只有五，所以可以直接用公式
后面用杜教筛

Code

#include<cstdio>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define N 5010000
#define NM 3012041
#define M 5001000
#define ll long long
#define mo 1000000007
#define min(a,b) ((a)<(b)?(a):(b))
using namespace std;
ll n,ans=0;
int k,bz1[NM],bz[N],hs[NM],phi[N],p[360000];
ll calc(ll x)
{
    x=x%mo;
    if(k==1) return (((x*(x+1))%mo)*500000004ll)%mo;
    else if(k==2) return (((x*(x+1))%mo*((x+x+1)%mo))%mo*166666668ll)%mo;
         else if(k==3) return (((x*x)%mo*(((x+1)*(x+1))%mo))%mo*250000002ll)%mo;
              else if(k==4) return ((((x*(x+1))%mo*((2ll*x+1)%mo)%mo)*(((3ll*x*x)%mo+(3ll*x)%mo-1+mo)%mo))%mo*233333335ll)%mo;
                   else return ((((x*x)%mo*(((x+1)*(x+1))%mo))%mo*(((2ll*x*x)%mo+(2ll*x)%mo-1+mo)%mo))%mo*83333334ll)%mo;
}
int hash(ll x)
{
    ll y=x%NM;
    while(bz1[y]!=0&&bz1[y]!=x) y=(y+1)%NM;
    return y;
}
ll dg(ll x)
{
    if(x<=M) return phi[x];
    int y=hash(x);
    if(bz1[y]) return hs[y];
    bz1[y]=x;
    ll ans=0;
    for(ll i=2,i1=2;i1<=x;i=i1+1,i1=min(x,x/(x/i)))
    {
        ans=(ans+(dg(x/i)*(i1-i+1))%mo)%mo;
        if(i1==x) break;
    }
    x%=mo;
    ans=(((x*(x+1))%mo*500000004ll)%mo-ans+mo)%mo;
    hs[y]=ans;return ans;
}
int main()
{
    phi[1]=1;
    fo(i,2,M)
    {
        if(!bz[i]) p[++p[0]]=i,phi[i]=i-1;
        fo(j,1,p[0])
        {
            int k=i*p[j];
            if(k>M) break;
            bz[k]=1;
            if(i%p[j]==0) {phi[k]=phi[i]*p[j];break;}
            phi[k]=phi[i]*(p[j]-1);
        }
    }
    fo(i,2,M) phi[i]=(phi[i-1]+phi[i])%mo;
    scanf("%lld%d",&n,&k);
    ll ls=0;
    for(ll d=1,d1=1;d1<=n;d=d1+1,d1=min(n,n/(n/d)))
    {
        ll an=dg(n/d),jy=calc(d1);
        an=(an+an-1+mo)%mo;
        ans=(ans+((jy-ls+mo)%mo*an)%mo)%mo;
        ls=jy;
        if(d1==n) break;
    }
    printf("%lld",ans);
}