51nod1847 奇怪的数学题

Description

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

计算下式

\sum_{i = 1}^{n} \sum_{j = 1}^{n} s g c d (i, j)^{k} \mod 2^{32}

$\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} [g c d (i, j) = d] \\ = \sum_{d = 1}^{n} f (d)^{k} (2 \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} φ (i) - 1) \end{aligned}

$\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{matrix} k \\ j \end{matrix}} (\binom{i}{j}) \\ = \sum_{j = 0}^{k} j! {\begin{matrix} k \\ j \end{matrix}} \sum_{i = 1}^{n} (\binom{i}{j}) \\ = \sum_{j = 0}^{k} j! {\begin{matrix} k \\ j \end{matrix}} (\binom{n + 1}{j + 1}) \\ = \sum_{j = 0}^{k} {\begin{matrix} k \\ j \end{matrix}} \frac{(n + 1)^{\underline{j + 1}}}{j + 1} \end{aligned}

$\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;
}