51Nod 1847 奇怪的数学题

奇怪的数学题

传送门

Description

Solution

考虑枚举$gcd$。

$$Ans=\sum_{d=1}^n(\frac{d}{minp(d)})^k\sum_{j=1}^n\sum_{k=1}^n[(j,k)=d]$$
$$Ans=\sum_{d=1}^n(\frac{d}{minp(d)})^k\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}[(j,k)=1]$$
$$Ans=\sum_{d=1}^n(\frac{d}{minp(d)})^k(2\sum_{i=1}^{\lfloor {n\over d }\rfloor}\varphi( i)-1)$$

$minp(d)$表示$d$的最小质因子。
一个很直观的思想就是分块，对于$(2\sum_{i=1}^{\lfloor {n\over d }\rfloor}\varphi( i)-1)$这一部分就是经典的杜教筛问题。
关于前面的部分，可以用$Min\_25$筛处理，但注意到$f(d)=(\frac{d}{minp(d)})$并非是个积性函数，其实只需维护两个函数即可，即维护

$S(n,i)=\sum_{i=2}^nf(d)^k[i的最小质因子不小于p_i或i是质数]$
$S'(n,i)=\sum_{i=2}^nd^k[i的最小质因子不小于p_i或i是质数]$
$g，g'$也是同理。

于是我们便能得到转移式

$$S(n,i)=\left\{ \begin{array}{} S(n,i+1)+\sum_{e\geq 1且\ p_i^{e+1}\ \leq n}\ \ p_i^{(e-1)k}[S'(\lfloor\frac{n}{p_i^e}\rfloor,i+1)-g'(p_i,i)]+p_i^{ek}& &{p_i^2\leq n} \\ S(n,i+1)& &{p_i^2\geq n} \end{array}\right.$$
$$S'(n,i)=\left\{ \begin{array}{} S'(n,i+1)+\sum_{e\geq 1且\ p_i^{e+1}\ \leq n}\ \ p_i^{ek}[S'(\lfloor\frac{n}{p_i^e}\rfloor,i+1)-g(p_i,i)]+p_i^{(e+1)k}& &{p_i^2\leq n} \\ S'(n,i+1)& &{p_i^2\geq n} \end{array}\right.$$

还有一个细节，由于没有逆元求自然数幂和要用第二类斯特林数，取模时最好用自然溢出。

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)

using namespace std;
typedef unsigned long long ll;
const ll N=12e5,M=1e6,mo=((ll)1<<32)-1,mm=mo+1,U=100;

ll n;
int p[N],k,oo,phi[N];
ll g[N],h[N],zh[N],w[N],f1[N],f2[N];
bool bz[N];
int id1[M],id2[M];
ll xs[U];
ll f[U][U];

inline ll ksm(ll o,ll t)
{
    ll y=1;
    for(;t;t>>=1,o=o*o)
    if(t&1)y=y*o;
    return y;
}

inline void prepare()
{
    f[1][1]=1; f[1][0]=0;
    fo(i,2,k){
        f[i][0]=0;
        fo(l,1,i)f[i][l]=(f[i-1][l-1]+f[i-1][l]*l);
    }
    fo(i,0,k)xs[i]=f[k][i];
}

inline ll qq(ll o)
{
    ll ans=0; ++o;
    fo(i,1,k){
        int ok=0; ll lj=1;
        fo(l,0,i)if((o-l)%(i+1)==0&&ok==0)lj=(lj*((o-l)/(i+1))),ok=1;
        else lj=lj*(o-l);
        ans=(ans+lj*xs[i]);
    }
    return ans;
}

ll ask(ll o)
{
    if(o<=M)return f1[o];
    ll k=n/o;
    if(f2[k])return f2[k];
    f2[k]=(o*(o+1)>>1);
    ll l=2,r;
    while(l<=o){
        r=o/(o/l);
        f2[k]=f2[k]-(r-l+1)*ask(o/l);
        l=r+1;
    }
    return f2[k];
}

int main()
{
    cin>>n>>k;
    prepare();
    fo(i,2,M){
        if(!bz[i])p[++oo]=i,phi[i]=i-1,zh[oo]=(zh[oo-1]+ksm(i,k));
        fo(j,1,oo)if(i*p[j]<=M){
            bz[i*p[j]]=true;
            phi[i*p[j]]=phi[i]*(p[j]-1);
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
        }else break;
    }
    phi[1]=1;
    fo(i,1,M)f1[i]=(f1[i-1]+phi[i]);
    ll P=sqrt(n); int u=0,m=0;
    fo(i,1,oo)if(p[i]<=P)u=i;else break;
    ll l=1,r;
    while(l<=n){
        ll len=n/l; r=n/len;
        if(len<=P)id1[len]=++m;else id2[r]=++m;
        g[m]=len-1; h[m]=qq(len)-1; w[m]=len;
        l=r+1;
    }
    fo(i,1,u){
        register ll ww=ksm(p[i],k),U=(ll)p[i]*p[i]; register int op;
        fo(j,1,m)if(U<=w[j])
        {
            op=(w[j]/p[i]<=P)?id1[w[j]/p[i]]:id2[n/(w[j]/p[i])];
            g[j]=g[j]-(g[op]-i+1);
            h[j]=h[j]-(h[op]-zh[i-1])*ww;
        }else break;
    }
    fd(i,u,1){
        register ll ww=ksm(p[i],k),U=(ll)p[i]*p[i];
        fo(j,1,m)if(U<=w[j]){
            register ll uu=w[j]/p[i];
            register int ok;
            for(register ll o=p[i],o3=ww,u=w[j]/o;o<=uu;o=o*p[i],o3=o3*ww,u=u/p[i])
            {
                ok=(u<=P)?id1[u]:id2[n/u];
                h[j]=h[j]+(h[ok]-zh[i])*o3;
                h[j]=h[j]+o3*ww;
            }
            for(register ll o1=1,o2=p[i],u=w[j]/p[i];o2<=uu;o1=o1*ww,o2=o2*p[i],u=u/p[i])
            {
                ok=(u<=P)?id1[u]:id2[n/u];
                g[j]=g[j]+(h[ok]-zh[i])*o1;
                g[j]=g[j]+o1*ww;
            }
        }else break;
    }
    l=1; ll ans=0;
    while(l<=n){
        r=n/(n/l);
        int o1,o2;
        o2=(r<=P)?id1[r]:id2[n/r];
        o1=l-1<=P?id1[l-1]:id2[n/(l-1)];
        ll cz=g[o2]-g[o1];
        ll ok=2*ask(n/r)-1;
        ans=ans+cz*ok;
        l=r+1;
    }
    ans=ans&mo;
    cout<<ans;
}