51nod1847 奇怪的数学题

Description

http://www.51nod.com/Challenge/Problem.html#problemId=1847

计算下式

\[\sum_{i=1}^n\sum_{j=1}^n sgcd(i,j)^k\mod 2^{32} \]

\(sgcd\) 函数表示第二大的 \(\texttt{common divisor}\)，当 \(gcd(i,j)=1\) 时 \(sgcd(i,j)=0\)，例如 \(sgcd(3,5)=0,sgcd(15,10)=1\)

\(n\le 10^9,k\le 50\)

Solution

写这题是为了复习模板

设 \(f(x)\) 表示 \(x\) 的最大真因子，原式改写成

\[\begin{aligned} \\&= \sum_{d=1}^n f(d)^k\sum_{i=1}^n\sum_{j=1}^n [gcd(i,j)=d] \\&=\sum_{d=1}^n f(d)^k (2\sum_{i=1}^{\lfloor \frac nd\rfloor}\varphi(i)-1) \end{aligned}\]

后半部分是杜教筛模板，外层套整除分块即可得到，那么仍需求出 \(f^k(i)\) 函数前缀和

函数本身并没有什么特别的积性，但是考察 \(\min25\) 筛法的第一个 \(\rm{DP}\) 的过程：依次枚举质因子并从每个 \(g(n,i-1)\) 中去掉 \(i\) 在这其中的贡献

形式化的表达式写作

\(g[n][i]=g[n][i-1]-(pri_i^k)\times(g[\lfloor\frac{n}{pri_i}]\rfloor[i-1]-g[pri_{i-1}][i-1])\)

其中 g[n/pri[i]][i-1]-g[pri[i-1]][i-1] 的含义恰是将那些最大质因子为 pri[i] 的数字去掉最大质因子之后剩下的数字的 \(k\) 次幂值之和，那么可以贡献给 \(f^k(pri_i)\)

另外由于要处理 \(sgcd(pri_i,pri_i)\) 的值，所以在 \(\min25\) 筛法的过程中需要计算范围内的质数数量，这相对 trivial

在 \(\min25\) 筛的初始状态需要赋值 \(g(n,0)=\sum\limits_{i=1}^n i^k\)，因为模数的原因不能使用依赖于逆元的 \(\rm Lagrange\) 插值

一种可能的解决方法是出于指数很小使用第二类斯特林数，即

\[\begin{aligned}\\& \sum_{i=1}^n i^k \\&=\sum_{i=1}^n\sum_{j=0}^k j!\begin{Bmatrix}k\\j\end{Bmatrix}\binom ij \\&=\sum_{j=0}^k j! \begin{Bmatrix}k\\j\end{Bmatrix}\sum_{i=1}^n \binom ij \\&=\sum_{j=0}^k j! \begin{Bmatrix}k\\j\end{Bmatrix}\binom {n+1}{j+1} \\&=\sum_{j=0}^k \begin{Bmatrix}k\\j\end{Bmatrix} \frac{(n+1)^{\underline{j+1}}}{j+1} \end{aligned}\]

由于连续 \(j+1\) 个数中必然有一个 \(j+1\) 的倍数，所以找到之除掉后就没有除法操作了

Code

#define uint unsigned int
const int N=2e6+10;
int pri[N],cnt,fl[N];
map<int,uint> mp;
uint phi[N],Str[100][100],res[N],num[N],f[N],sum[N];
int R[N],tot,n,k,bl,id1[N],id2[N];
inline uint getphi(int n){
    if(n<=2e6) return phi[n];
    if(mp.count(n)) return mp[n];
    uint ans=(n&1)?(uint)(n+1)/2*n:(uint)n/2*(n+1);
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(uint)(r-l+1)*getphi(n/l);
    }
    return mp[n]=ans;
}
inline int getid(int x){
    if(x<=bl) return id1[x];
    else return id2[n/x];
}
inline uint calc(int n){
    int pos=getid(n);
    return res[pos]+num[pos];
}
inline uint S(int n,int k){
    uint ans=0;
    rep(j,0,k){
        uint tmp=1,now=n+1;
        bool fl=0;
        rep(i,0,j){
            if(fl==0&&now%(j+1)==0) tmp*=now/(j+1),fl=1;
            else tmp*=now;
            --now;
        }
        ans+=tmp*Str[k][j];
    }
    return ans;
}
inline uint ksm(uint x,int y){
    uint res=1;
    for(;y;y>>=1,x=x*x) if(y&1) res*=x;
    return res;
}
signed main(){
    n=2e6; phi[1]=1;
    for(int i=2;i<=n;++i){
        if(!fl[i]) phi[i]=i-1,pri[++cnt]=i;
        for(int j=1;j<=cnt&&1ll*i*pri[j]<=n;++j){
            fl[i*pri[j]]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }else{
                phi[i*pri[j]]=phi[i]*phi[pri[j]];
            }
        }
    }
    for(int i=1;i<=n;++i) phi[i]+=phi[i-1];
    n=read(); k=read(); bl=sqrt(n); 
    Str[1][1]=1; 
    rep(i,2,50){
        rep(j,1,i){
            Str[i][j]=Str[i-1][j]*j+Str[i-1][j-1];
        }
    }
    for(int l=1,r;l<=n;l=r+1){
        R[++tot]=n/l; r=n/R[tot];
        f[tot]=S(R[tot],k)-1;
        num[tot]=R[tot]-1;
        // Do not consider ones so subtract them
        if(R[tot]<=bl) id1[R[tot]]=tot; else id2[n/R[tot]]=tot;
    }
    for(int j=1;j<=cnt;++j){
        uint now=ksm(pri[j],k);
        sum[j]=sum[j-1]+now;
        for(int i=1;i<=tot&&1ll*pri[j]*pri[j]<=R[i];++i){
            int pos=getid(R[i]/pri[j]);
            f[i]-=(f[pos]-sum[j-1])*now;
            res[i]+=f[pos]-sum[j-1];
            num[i]-=num[pos]-(j-1); // Numbers except from primes
        }
    }
    uint ans=0;
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans+=(2*getphi(n/l)-1)*(calc(r)-calc(l-1));
    }
    print(ans);
    return 0;
}

代码不附带缺省源，如果需要可以去【遇到困难睡大觉】一文中粘贴并去掉 #define int long long

感觉最近代码越来越蓬松了，都快不可读了