【数论】[51Nod 1847] 奇怪的数学题

Description

求$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}sgcd(i,j)^m$
其中$sgcd(i,j)$表示i,j的次大公约数

答案对2^32取模
$n\leq 10^9,m\leq 50$
特别的,若$gcd(i,j)=1$，则$sgcd(i,j)=0$

Solution

次大公约数实际上就是最大公约数的最大因子

也就是最大公约数除以最大公约数的最小质因子

根据这个思路，我们可以用洲阁筛或者min_25筛
原式化一下
令$f(d)$表示$d$的最大因子的m次方

$$Ans=\sum\limits_{d=2}^{n}f(d)\sum\limits_{i=1}^{\lfloor{n\over d}\rfloor}\sum\limits_{i=1}^{\lfloor{n\over d}\rfloor}[\gcd(i,j)=1]$$

$$=\sum\limits_{d=2}^{n}f(d)((2\sum\limits_{i=1}^{\lfloor{n\over d}\rfloor}φ(i))-1)$$

后面的可以杜教筛求出

对于前面的，我们想要求$f(d)$的前缀和，设为$S(d)$
这与洲阁筛(min_25筛)的思路很接近

在做min_25筛的时候我们有DP状态$g(t,n)=\sum\limits_{i=2}^{n}i^m[i为质数或i的最大质因子大于第t个质数]$

这里不细讲min_25筛了

因为分块的时候只有$2\sqrt n$个状态是有用的
在做min_25筛的时候，我们可以边DP边将答案累加进S

n以内最小质因子为$P_i$（表示第i个质数）的数的答案就是$g(\lfloor {n\over P_i}\rfloor,i-1)*P_i^m$

那么需要累加进S的就是$g(\lfloor {n\over P_i}\rfloor,i-1)$

注意1的情况要另外处理

由于模数不是质数，一开始预处理自然数幂和要用斯特林数的方法
注意要unsigned long long

Code

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define N 320005
#define M 54
#define mo 4294967296
#define LL unsigned long long
using namespace std;
int n,m,n1,l,pr[N];
LL sp[N],g[N],s1[N],sr[N],st[M][M],id[2*N],id1[2*N],id2[2*N],f[2*N],sm[M],h[N],ct[N],phi[N],si[N],hs[N];
bool bz[N];
LL ksm(LL k,LL n)
{
    LL s=1;
    for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;
    return s;
}
void prp()
{
    s1[1]=1;
    phi[1]=si[1]=1;
    fo(i,2,N-5)
    {
        if(!bz[i]) s1[i]=sr[i]=ksm(i,m),pr[++l]=i,ct[i]=1,phi[i]=i-1;
        for(int j=1;j<=l&&i*pr[j]<=N-5;j++)
        {
            bz[i*pr[j]]=1;
            s1[i*pr[j]]=s1[i]*sr[pr[j]]%mo;
            if(i%pr[j]==0) 
            {
                phi[i*pr[j]]=phi[i]*pr[j]%mo;
                break;
            }
            phi[i*pr[j]]=phi[i]*(pr[j]-1);
        }
        sp[i]=(sp[i-1]+sr[i])%mo;
        ct[i]=(ct[i-1]+ct[i])%mo;
        s1[i]=(s1[i-1]+s1[i])%mo;
        si[i]=(si[i-1]+phi[i])%mo;
    }
}
int num(int k)
{
    if(!k) return 0;
    return(k<=n1)?id1[k]:id2[n/k];
}
LL sum(LL n)
{
    sm[0]=n;
    fo(i,1,m)
    {
        sm[i]=0;
        if(n+1>=i+1)
        {
            sm[i]=1;
            bool pd=0;
            fod(j,n+1,n-i+1) 
            {
                if(j%(i+1)==0&&!pd) pd=1,sm[i]=sm[i]*(LL)(j/(i+1))%mo;
                else sm[i]=sm[i]*(LL)j%mo;
            }
        }
        LL v=(i%2)?1:-1;
        fo(j,0,i-1)
        {
            sm[i]=(sm[i]+v*st[i][j]*sm[j]%mo+mo)%mo;
            v=-v;
        }
    }
    return sm[m];
}
LL get(LL k)
{
    if(k<=N-5) return si[k];
    if(hs[n/k]) return hs[n/k];
    LL d=2,s=k*(k+1)/2%mo;
    while(d<=k)
    {
        LL dm=k/d,d1=k/dm;
        s=(s-get(dm)*(d1-d+1)%mo+mo)%mo;
        d=d1+1;
    }
    return(hs[n/k]=s);
}
int main()
{
    cin>>n>>m;
    n1=sqrt(n);
    prp();
    st[0][0]=1;
    fo(i,1,m)
    {
        fo(j,1,m)
        {
            st[i][j]=((LL)(i-1)*st[i-1][j]%mo+st[i-1][j-1]%mo)%mo;
        }
    }
    LL d=1;
    while(d<=n)
    {
        LL dm=n/d,d1=n/dm;
        id[++id[0]]=dm;
        if(dm<=n1) id1[dm]=id[0];
        else id2[d1]=id[0];
        if(dm<=N-5) g[id[0]]=s1[dm]-1;
        else g[id[0]]=sum(dm)-1;
        h[id[0]]=dm-1;
        d=d1+1;
    }
    fo(j,1,l)
    {
        if(pr[j]>n1) break;
        int cnt=0;
        fo(i,1,id[0])
        {
            if(pr[j]*pr[j]>id[i]) break;
            cnt=i;
            int p=num(id[i]/pr[j]);
            f[i]=(f[i]+g[p]-sp[pr[j]-1]+mo)%mo;
            g[i]=(g[i]-sr[pr[j]]*(g[p]-sp[pr[j]-1]+mo)%mo+mo)%mo;
            h[i]=(h[i]-(h[p]-ct[pr[j]-1]+mo)%mo+mo)%mo;
        }
    }
    fo(i,1,id[0]) f[i]=(f[i]+h[i])%mo;
    LL ans=0;
    d=2;
    while(d<=n)
    {
        LL dm=n/d,d1=n/(n/d);
        ans=(ans+(f[num(d1)]-f[num(d-1)]+mo)%mo*((LL)2*get(dm)%mo-1)%mo)%mo;
        d=d1+1;
    } 
    printf("%lld\n",(ans+mo)%mo);
}