【HDU 6607】Easy Math Problem(杜教筛+min25筛+拉格朗日插值)

Problem Description

One day, Touma Kazusa encountered a easy math problem. Given n and k, she need to calculate the following sum modulo 1e9+7.
${\sum }_{i=1}^n{\sum }_{j=1}^{n}gcd(i,j{)}^{k}lcm(i,j)[gcd(i,j)\in prime]\%(1e9+7)$

However, as a poor student, Kazusa obviously did not, so Touma Kazusa went to ask Kitahara Haruki. But Kitahara Haruki is too busy, in order to prove that he is a skilled man, so he threw this problem to you. Can you answer this easy math problem quickly?

Input

There are multiple test cases.（T=5） The first line of the input contains an integer T, indicating the number of test cases. For each test case:

There are only two positive integers n and k which are separated by spaces.

$1\le n\le 10^{10}$
$1\le k\le 100$

Output

An integer representing your answer.

Sample Input

1

10 2
Sample Output

2829
come from 2019 Multi-University Training Contest 3

其实我觉得这题是不难的，和51nod 1238上其实有差不多的地方，这题我去年的时候还做过，唯一不同的是次数和要求的是质数项上的和，这里搞次求值用的是拉格朗日插值，质数项求和用的min25，其实看到质数求和大概率是想到了要用min25。
写起来也比较麻烦，对于我来说，有点力不从心。
代码：

#include<stdio.h>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<queue>
#include<cmath>
#include<vector>
#include<cmath>
#define ll long long
#define INF 0x7f7f7f7f
using namespace std;
const int mod=1000000007;
const int maxx =5e6;
const int maxn =115;
ll inv4=250000002,inv6=166666668;
ll fac[205],inv[205];
ll n;
bool isp[maxx];
int pr[maxx],cnt;
int phi[maxx];
ll sum[maxx];
inline ll P(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;
}
void init()
{
    phi[1]=1;
    for(int i=2;i<maxx;i++)
    {
        if(!isp[i])pr[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&(ll)i*pr[j]<maxx;j++)
        {
            int p=i*pr[j];
            isp[p]=true;
            if(i%pr[j])phi[p]=phi[i]*(pr[j]-1);
            else
            {
                phi[p]=phi[i]*pr[j];
                break;
            }
        }
    }
    for(int i=1;i<maxx;i++)
    {
        sum[i]=sum[i-1]+(ll)i*i%mod*phi[i]%mod;
        if(sum[i]>=mod)sum[i]-=mod;
    }
}
inline ll cube(ll x)
{
    x%=mod;
    x=x*(x+1)%mod;
    return x*x%mod*inv4%mod;
}
inline ll square(ll x)
{
    x%=mod;
    return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
ll M[maxx];
ll work(ll x)
{
    if(x<maxx)return sum[x];
    if(M[n/x])return M[n/x];
    ll res=cube(x);
    for(ll l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        res-=(square(r)-square(l-1))*work(x/l)%mod;
        res%=mod;
    }
    if(res<0)res+=mod;
    return M[n/x]=res;
}
int k;
ll f[205];
void _init(int tot)
{
    fac[0]=1;
    for(int i=1;i<=tot;i++)
        fac[i]=fac[i-1]*i%mod;
    inv[tot]=P(fac[tot],mod-2);
    inv[0]=1;
    for(int i=tot-1;i>=1;i--)
        inv[i]=inv[i+1]*(i+1)%mod;
    for(int i=1;i<=tot;i++)f[i]=(f[i-1]+P(i,k))%mod;
}
inline ll cal(ll x)
{
    if(x<=k+2)return f[x];
    int tot=k+2;
    ll ans=0;
    ll now=1;
    for(int i=1;i<=tot;i++)now=now*((x-i)%mod)%mod;
    for(int i=1;i<=tot;i++)
    {
        ll inv1=P(x-i,mod-2);
        ll inv2=inv[i-1]*inv[tot-i]%mod;
        if((tot-i)&1)
            inv2=mod-inv2;
        ll temp=now*inv1%mod;
        temp=temp*f[i]%mod*inv2%mod;
        ans+=temp;
        if(ans>=mod)ans-=mod;
    }
    return ans;
}
ll N;
ll w[maxx],tot;
int id1[maxx],id2[maxx];
ll g[maxx];
ll sp[maxx];
void min25()
{
    N=sqrt(n);
    tot=0;
    for(ll i=1,last;i<=n;i=last+1)
    {
        ll now=n/i;
        w[++tot]=now;
        last=n/now;
        if(now<=N)id1[now]=tot;
        else id2[last]=tot;
        g[tot]=cal(w[tot])-1;
        //cout<<g[tot]<<endl;
    }
    int num=1;
    while(pr[num]<=N)
    {
        sp[num]=(sp[num-1]+P(pr[num],k))%mod;
        num++;
    }
    for(int i=1;i<num;i++)
    {
        for(int j=1;j<=tot&&(ll)pr[i]*pr[i]<=w[j];j++)
        {
            ll now=w[j]/pr[i];
            int k=(now<=N?id1[now]:id2[n/now]);
            g[j]-=(g[k]-sp[i-1])*(sp[i]-sp[i-1])%mod;
            g[j]%=mod;
        }
    }
}
int main()
{
    init();
    int t;cin>>t;
    while(t--)
    {
        scanf("%lld%d",&n,&k);
        k++;
        _init(k+2);
        min25();
        //cout<<g[1]<<endl;
        for(int i=1;i<=N;i++)M[i]=0;
        ll ans=0;
        for(ll l=1,r;l<=n;l=r+1)
        {
            r=n/(n/l);
            int R=(r<=N?id1[r]:id2[n/r]);
            int L=(l-1<=N?id1[l-1]:id2[n/(l-1)]);
            ans+=(g[R]-g[L])*work(n/l)%mod;
            ans%=mod;
        }
        if(ans<0)ans+=mod;
        printf("%lld\n",ans);
    }
    return 0;
}