51nod 1847 奇怪的数学题

Description

给出 N,K ,请计算下面这个式子:
∑Ni=1∑Nj=1sgcd(i,j)k
其中,sgcd(i, j)表示(i, j)的所有公约数中第二大的,特殊地,如果gcd(i, j) = 1, 那么sgcd(i, j) = 0。
考虑到答案太大,请输出答案对2^32取模的结果.
1≤N≤109,1≤K≤50
样例解释:
因为gcd(i, j)=1时sgcd(i,j)=0对答案没有贡献,所以我们只考虑gcd(i,j)>1的情况.
当i是2时,j是2时,sgcd(i,j)=1,它的K次方是1
当i是2时,j是4时,sgcd(i,j)=1,它的K次方是1
当i是3时,j是3时,sgcd(i,j)=1,它的K次方是1
当i是4时,j是2时,sgcd(i,j)=1,它的K次方是1
当i是4时,j是4时,sgcd(i,j)=2,它的K次方是8
当i是5时,j是5时,sgcd(i,j)=1,它的K次方是1

Input

两个数字N,K

Output

⼀⾏⼀个数表⽰答案。

Sample Input

5 3

Sample Output

13

Solution

min_25
是一个神奇的东西。

设f(d)表示d的第二大的约数，也就是d除以d的最小质因数。
显然答案可以表示为

$$ans=\sum_{d=1}^nf(d)*\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=d]$$
将后面套路搞一搞可以得到
$$ans=\sum_{d=1}^nf(d)*2\sum_{i=1}^n\phi(i)-1$$
后面用杜教筛可以轻松解决，那么前面就用min_25解决就行了。

p[i]表示第i个质数
设g(n,i)表示$\sum_{j是质数或j的最小质因子不小于p[i]} f[j]^k$
设h(n,i)表示$\sum_{j是质数或j的最小质因子不小于p[i]} j^k$
转移很套路，如下：
当$p[i]^2>n$时，显然$g(n,i)=g(n,i-1)，h(n,i)=h(n,i-1)$
否则
设sp[i]表示第一个到第i个质数的k次方和

$$g(n,i)=g(n,i+1)+\sum_{p[i]^{e+1}<=n}p[i]^{(e-1)k}(h({n\over p[i]^e},i+1)-sp[i])+p[i]^{ek}$$
$$h(n,i)=h(n,i+1)+\sum_{p[i]^{e+1}<=n}p[i]^{ek}(h({n\over p[i]^e},i+1)-sp[i])+p[i]^{(e+1)k}$$
设p0表示 $\sqrt n$以内质数的个数。
初始时g(n,p0+1)的值显然是n以内质数的个数，h(n,p0+1)的值显然是n以内质数的k次方的和，发现不知道怎么搞。

于是再弄一次min_25，将初始的g和h筛出来。
这次
设g(n,i)表示$\sum_{j是质数或j的最小质因子不小于p[i]} 1$
设h(n,i)表示$\sum_{j是质数或j的最小质因子不小于p[i]} j^k$
转移如下

$$g(n,i)=g(n,i-1)-(g({n\over p[i]},i-1)-i+1)$$
$$h(n,i)=h(n,i-1)-p[i]^{k}(h({n\over p[i]},i-1)-sp[i-1])$$
初始值即g(n,0)和h(n,0)分别为n和自然数幂和。（自然数幂和要用第二类斯特林数）
转移完成后，得到的g和h就是上面第一个min_25筛的初始值。
总得来说，这题只需要两次min_25筛加一次杜教筛加一次自然数幂和还有最后的分块求答案就能做了。

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<map>
#define fo(i,a,b) for(ll i=a;i<=b;i++)
#define fd(i,a,b) for(ll i=a;i>=b;i--)
#define N 1010000
#define ll unsigned long long
using namespace std;
map<ll,ll> bphi;
ll ans,n,k,bz[N],p[N],phi[N],gn,g[N],h[N],g0,sp[N],ss[60][60],p0,m,num[N],d1[N],d2[N];
ll mi(ll a,ll b)
{
    ll c=1;
    for(;b;b/=2,a=a*a) if(b%2==1) c=c*a;
    return c;
}
ll gphi(ll n)
{
    if(n<N) return phi[n];
    if(bphi[n]>0) return bphi[n];
    ll ans;
    if(n%2==0) ans=n/2*(n+1);
    else ans=(n+1)/2*n;
    for(ll l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=ans-(r-l+1)*gphi(n/l);
    }
    bphi[n]=ans;
    return ans;
}
ll getk(ll n)
{
    ll ans=0;
    fo(i,0,k)
    if(n>=i)
    {
        ll C=1;
        int flag=1;
        fo(j,n-i+1,n+1)
        {
            if(j%(i+1)==0&&flag) C=j/(i+1)*C,flag=0;
            else C=C*j;
        }
        ans=ans+ss[k][i]*C;
    }
    else break;
    return ans;
}
int main()
{
    scanf("%llu%llu",&n,&k);
    ss[0][0]=1;
    fo(i,1,k+1) fo(j,1,k+1) ss[i][j]=ss[i-1][j-1]+ss[i-1][j]*j;
    gn=(ll)sqrt(n);
    fo(i,2,N-1)
    {
        if(!bz[i])
        {
            p[++p[0]]=i,phi[i]=i-1;
            sp[p[0]]=sp[p[0]-1]+mi(i,k);
            if(i<=gn) p0=p[0];
        }
        fo(j,1,p[0])
        {
            ll k=i*p[j];
            if(k>=N) break;
            bz[k]=1;
            if(i%p[j]==0)
            {
                phi[k]=phi[i]*p[j];
                break;
            }
            phi[k]=phi[i]*phi[p[j]];
        }
    }
    phi[1]=1;
    fo(i,2,N-1) phi[i]+=phi[i-1];
    for(ll l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        if(n/l<=gn) d1[n/l]=++m;else d2[r]=++m;
        num[m]=n/l;
        g[m]=num[m]-1;h[m]=getk(num[m])-1;
    }
    fo(i,1,p0)
    {
        ll pk=mi(p[i],k);
        ll pf=1ll*p[i]*p[i];
        fo(j,1,m)
        if(pf<=num[j])
        {
            ll op=(num[j]/p[i])<=gn?d1[num[j]/p[i]]:d2[n/(num[j]/p[i])];
            g[j]=g[j]-(g[op]-i+1ll);
            h[j]=h[j]-pk*(h[op]-sp[i-1]);
        }
        else break;
    }
    fd(i,p0,1)
    {
        ll pk=mi(p[i],k);
        ll pf=1ll*p[i]*p[i];
        fo(j,1,m)
        if(pf<=num[j])
        {
            for(ll pmk=1,pm=p[i];pm*p[i]<=num[j];pmk=pmk*pk,pm=pm*p[i])
            {
                ll op=(num[j]/pm)<=gn?d1[num[j]/pm]:d2[n/(num[j]/pm)];
                g[j]=g[j]+pmk*pk+pmk*(h[op]-sp[i]);
                h[j]=h[j]+pmk*pk*pk+pmk*pk*(h[op]-sp[i]);
            }
        }
        else break;
    }
    for(ll l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ll op1=(r<=gn)?d1[r]:d2[n/r];
        ll op2=(l-1)<=gn?d1[l-1]:d2[n/(l-1)];
        ll jy=gphi(n/l);
        jy=jy*2-1;
        ans=ans+(g[op1]-g[op2])*jy;
    }
    printf("%llu\n",ans&((1ll<<32)-1));
}