51nod1847 奇怪的数学题

Description

给出 $N,K$ ,请计算下面这个式子:
$sum _{i=1}^n sum_{j=1}^n sgcd(i,j)^k$
其中,$sgcd(i, j)$表示$(i, j)$的所有公约数中第二大的,特殊地,如果$gcd(i, j) = 1$, 那么$sgcd(i, j) = 0$。
考虑到答案太大,请输出答案对$2^{32}$取模的结果.
$1 leq N leq 10^9,1 leq K leq 50$

Solution

设$f(x)$为$x$的次大的因数，那么原式可化为

egin{align}
& sum _{i=1}^n sum_{j=1}^n sgcd(i,j)^k\
= & sum _{i=1}^n sum _{j=1}^n f((i,j))^k\
= & sum _{d=1}^n f^k(d) sum _{i=1}^{lfloor frac{n}{d} floor }sum _{j=1}^{lfloor frac{n}{d} floor}[(i,j)=1]\
= & sum _{d=1}^n f^k(d) left ( 2sum _{i=1}^{lfloor frac {n}{d} floor }varphi(i)-1 ight )
end{align}

再定义$G_n=sum _{i=1}^n f(i)$，定义Min_25中的$g$函数为符合要求的$x^k$之和

发现在做$g$转移时：

$$g_{i,j}=g_{i-1,j}-p_j^k(g_{lfloor frac{i}{p_j} floor ,j-1 }-g_{p_{j-1},j-1})$$

括号内的内容即为该数本身除以它的最小质因数，就是题中的次大因数，所以在转移$g$的同时记录括号中的数的和，这样可以得到所有整除点的$G$值

$g$函数初始化的时候要求自然数幂之和，使用第二类斯特林数和普通幂的转化

欧拉函数可以用杜教筛得到整除点的值

#include<unordered_map>
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
int tot,cnt,id1[100005],id2[100005];
unsigned int phi[1000005],S[60][60],pw[1000005],w[200005],g1[200005],g0[200005],sp[1000005],G[1000005],ans,pre,now;
long long n,K,prime[1000005],sq;
bool vst[1000005];
unordered_map<int,unsigned int>mp;
inline int read(){
    int f=1,w=0;
    char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')w=(w<<1)+(w<<3)+ch-'0',ch=getchar();
    return f*w;
}
unsigned int ksm(unsigned int a,unsigned int p){
    unsigned int ret=1;
    while(p){
        if(p&1)ret*=a;
        a*=a,p>>=1;
    }
    return ret;
}
unsigned int calc(long long x){
    unsigned int ret=0,temp;
    for(unsigned int i=1;i<=K;i++){
        int t=(x-i+1)%(i+1);
        temp=1;
        for(long long j=x-i+1;j<=x+1;j++,++t){
            if(t>=i+1)t-=i+1;
            if(!t)temp*=(unsigned int)j/(i+1);
            else temp*=(unsigned int)j;
        }
        ret+=S[K][i]*temp;
    }
    return ret;
}
unsigned int djs(long long x){
    if(x<=1000000)return phi[x];
    if(mp.find(x)!=mp.end())return mp[x];
    unsigned int ret=0;
    if(x&1)ret=(x+1)/2*x;
    else ret=x/2*(x+1);
    for(unsigned int i=2;i<=x;){
        unsigned int j=x/(x/i);
        ret-=djs(x/i)*(j-i+1),i=j+1;
    }
    return mp[x]=ret;
}
int main(){
    n=read(),K=read(),sq=sqrt(n),phi[1]=1;
    for(int i=2;i<=1000000;i++){
        if(!vst[i])prime[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot&&i*prime[j]<=1000000;j++){
            vst[i*prime[j]]=true;
            if(!(i%prime[j])){phi[i*prime[j]]=phi[i]*prime[j];break;}
            else phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    }
    for(int i=1;i<=1000000;i++)phi[i]+=phi[i-1];
    for(int i=1;i<=K;i++){
        S[i][1]=1;
        for(unsigned int j=2;j<=i;j++)S[i][j]=S[i-1][j]*j+S[i-1][j-1];
    }
    for(int i=1;i<=tot;i++)pw[i]=ksm(prime[i],K),sp[i]=sp[i-1]+pw[i];
    for(int i=1;i<=n;){
        int j=n/(n/i);
        w[++cnt]=n/i,g1[cnt]=calc(w[cnt])-1,g0[cnt]=w[cnt]-1,n/i<=sq?id1[n/i]=cnt:id2[n/(n/i)]=cnt,i=j+1;
    }
    for(int i=1;i<=tot;i++)for(int j=1;j<=cnt&&1ll*prime[i]*prime[i]<=w[j];j++){
        int id=w[j]/prime[i]<=sq?id1[w[j]/prime[i]]:id2[n/(w[j]/prime[i])];
        g1[j]-=pw[i]*(g1[id]-sp[i-1]),g0[j]-=g0[id]-(i-1),G[j]+=g1[id]-sp[i-1];
//        cout<<g0[j]<<" ";
    }
    for(int i=2;i<=n;){
        int j=n/(n/i),id=j<=sq?id1[j]:id2[n/j];
        now=G[id]+g0[id],ans+=(djs(n/i)*2-1)*(now-pre),pre=now,i=j+1;
    }
    printf("%u",ans);
    return 0;
}//3656808373