自然数幂求和+杜教筛

自然数幂求和

1.

高斯消元

（观察）可得肯定是一个多项式

推出来前几项 消元就行了

2.

倍增

（from http://blog.csdn.net/werkeytom_ftd/article/details/50741165#reply）

inline ll f(ll x,int y){
    if(x<N&&ff[x][y])return ff[x][y];
    if((it=F[y].find(x))!=F[y].end())return it->second;
    if(!x)return 0;
    if(x&1){
        if(x<N)return ff[x][y]=(f(x-1,y)+pw(x,y)+mod)%mod;
        return F[y][x]=(f(x-1,y)+pw(x,y)+mod)%mod;
    }
    ll res=0;
    for(ll j=0;j<=y;j++)
        res=(res+pw(x/2,y-j)*C[y][j]%mod*f(x/2,j))%mod;
    if(x<N)return ff[x][y]=(res+f(x/2,y)+mod)%mod;
    return F[y][x]=(res+f(x/2,y)+mod)%mod;
}

//By SiriusRen
#include <map>
#include <cstdio>
using namespace std;
const int N=2500000,mod=1000000007;
typedef long long ll;
ll ans,n;
int C[15][15],k,ff[N][6],tot,phi[N],vis[N],prime[N];
map<ll,ll>mp,F[6];map<ll,ll>::iterator it;
void shai(){
    phi[1]=1;
    for(ll i=2;i<N;i++){
        if(!vis[i])prime[++tot]=i,phi[i]=i-1;
        for(ll j=1;i*prime[j]<N&&j<=tot;j++){
            phi[i*prime[j]]=phi[i]*(prime[j]-1),vis[i*prime[j]]=1;
            if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
        }(phi[i]+=phi[i-1])%=mod;
    }
}
inline ll pw(ll x,ll y){
    ll res=1;x%=mod;
    while(y){
        if(y&1)res=res*x%mod;
        x=x*x%mod,y>>=1;
    }return res;
}
ll get_phi(ll x){
    if(x<N)return phi[x];
    if((it=mp.find(x))!=mp.end())return it->second;
    ll res=((x%mod)*((x+1)%mod))%mod*pw(2,mod-2)%mod;
    for(ll i=2,last;i<=x;i=last+1)
        last=x/(x/i),res=(res-get_phi(x/i)*((last-i+1)%mod))%mod;
    return mp[x]=(res+mod)%mod;
}
inline ll f(ll x,int y){
    if(x<N&&ff[x][y])return ff[x][y];
    if((it=F[y].find(x))!=F[y].end())return it->second;
    if(!x)return 0;
    if(x&1){
        if(x<N)return ff[x][y]=(f(x-1,y)+pw(x,y)+mod)%mod;
        return F[y][x]=(f(x-1,y)+pw(x,y)+mod)%mod;
    }
    ll res=0;
    for(ll j=0;j<=y;j++)
        res=(res+pw(x/2,y-j)*C[y][j]%mod*f(x/2,j))%mod;
    if(x<N)return ff[x][y]=(res+f(x/2,y)+mod)%mod;
    return F[y][x]=(res+f(x/2,y)+mod)%mod;
}
signed main(){
    shai(),scanf("%lld%d",&n,&k);
    for(int i=0;i<=10;i++){
        C[i][0]=C[i][i]=1;
        for(int j=1;j<i;j++)
            C[i][j]=C[i-1][j-1]+C[i-1][j];
    }
    for(ll i=1,last;i<=n;i=last+1)
        last=n/(n/i),ans=(ans+(f(last,k)-f(i-1,k)+mod)*(get_phi(n/i)*2-1))%mod;
    printf("%lld
",ans);
}

O(k2lognlogk)