51nod 1847 奇怪的数学题

51nod 1847 奇怪的数学题

原题链接：
https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1847

定义：

$sgcd(a,b)$为$a$与$b$的次大公约数

$f(a)$为$a$的次大约数 $$sgcd(a,b)=f\big(gcd(a,b)\big)$$

特别的：$f(1)=0$

题目要求计算： $$\sum_{i=1}^n\sum_{j=1}^nf\big(gcd(i,j)\big)^k$$

推导：

$$\sum_{i=1}^n\sum_{j=1}^nf\big(gcd(i,j)\big)\\=\sum_{d=1}^nf(d)^k\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=d]\\=\sum_{d=1}^nf(d)^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[i\perp j]\\=\sum_{d=1}^nf(d)^k\Big(-1+2\sum_{i=1}^{\lfloor\frac {n}{d}\rfloor }\varphi(i)\Big)$$

如果可以分段计算$f$的前缀和与$\varphi$的前缀和。就可以快速计算答案。

对于$\varphi$有： $$\varphi*id_0=id$$

经典杜教筛。

对于计算： $$S_f(n)=\sum_{i=1}^nf(i)$$

定义：$C(t,n)$为，$[1,n]$上不可以被前t个素数整除的数字的$k$次幂之和

对于第$t$个素数$P_t$，$[1,n]$上最小素因子为$P_t$的数字必然有：

. 这些数是$P_t$的倍数

. 这些数不能被前$t-1$个素数整除

所以：最小素因子为$P_t$的数字$k$次幂之和为： $$C(t-1,\Big\lfloor\frac{n}{P_t}\Big\rfloor)P_t^k$$

所以： $$C(t ,n)=C(t-1,n)-C(t-1,\Big\lfloor\frac{n}{P_t}\Big\rfloor)P_t^k\\S_f(n)=\sum_{t\geq 1}C(t-1,\Big\lfloor\frac{n}{P_t}\Big\rfloor)$$

定义： $$S(t,n)=S(t-1,n)+C(t-1,\Big\lfloor\frac{n}{P_t}\Big\rfloor)$$

$$S(0,n)=0$$

所以： $$S_f(n)=S(\pi(n) ,n)$$

显然当 $P_{t+1}>n:$ $$C(t,n)=1$$

则$P_t^2>n$时： $$C(t ,n)=C(t-1,n)-P_t^k$$

那么为们只更新 $$C(t-1,\Big\lfloor\frac{n}{P_t}\Big\rfloor)>1$$

这个情况。多开一个数组记录最早出现$1$的哪个素数标号，预处理合适的前若干个素数$k$次幂之和。

总时间复杂度: $$T(n)=\sum_{1<P\leq n}^n\Big( O\big(\sqrt {\frac{n}{P}}\big)+O\big(\sqrt{P}\big)\Big)=O\Big(\frac{n^{\frac 4 3}}{log\ n}\Big)$$

具体方法类似于：http://blog.csdn.net/ZLH_HHHH/article/details/78020896

下面是AC代码：

#include <algorithm>
#include <string.h>
#include <stdio.h>
#include <cmath>
#define MAXN 40000
#define LIM 1000000
using namespace std;

typedef  unsigned int uint;

uint size;
const uint INF=0x3f3f3f3f;
struct mat
{
    uint A[53][53];

    mat()
    {
        memset(A,0,sizeof A);
    }

    void C()
    {
        A[0][0]=1;
        for(uint i=1;i<size;i++)
        {
            A[i][0]=1;
            for(uint j=1;j<=i;j++)  A[i][j]=A[i-1][j]+A[i-1][j-1];
        }
        for(uint i=0;i<size;i++) A[size][i]=A[size-1][i];
        A[size][size]=1;
    }

    mat operator *(const mat &a)const
    {
        mat b;
        for(uint i=0;i<=size;i++)
            for(uint j=0;j<=size;j++)
                for(uint k=0;k<=size;k++)
                    b.A[i][j]+=A[i][k]*a.A[k][j];
        return b;
    }
}bit[33];

uint B[MAXN];
uint b[MAXN];
uint Sp[LIM];
uint _Sp[LIM];
uint C[MAXN];
uint c[MAXN];
uint Ik[MAXN];
uint prim[MAXN];
uint S[MAXN];
uint s[MAXN];
uint tmp[MAXN];
uint pai[MAXN];
uint paik[MAXN];
uint f[MAXN];
uint _TMP[55],TMP[55];
uint test[100];
bool vis[MAXN];
uint Pow(uint a,uint b)
{
    uint tmp=1;
    while(b)
    {
        if(b&1)
            tmp=tmp*a;
        a=a*a;
        b>>=1;
    }
    return tmp;
}
uint clat_Sk(uint n)
{
    n--;
    uint k=0;
    for(uint i=0;i<=size;i++)TMP[i]=1;
    while(n)
    {
        if(n&1)
        {
            memset(_TMP,0,sizeof _TMP);
            for(uint i=0;i<=size;i++)
                for(uint j=0;j<=size;j++)
                    _TMP[i]+=bit[k].A[i][j]*TMP[j];
            for(uint i=0;i<=size;i++)TMP[i]=_TMP[i];
        }
        n>>=1;
        k++;
    }
    return TMP[size];
}
uint clat_S(uint n)
{
    if(n&1)return n*((n+1)>>1);
    return (n>>1)*(n+1);
}
void clat_Sp(uint n,uint d)
{
    uint m=sqrt(n)+1,ans=0;
    for(uint L=1;L<m;L++)   ans+=(n/L-n/(L+1))*Sp[L];
    for(uint i=n/m;i>1;i--)
    {
        uint u=n/i;
        if(u<LIM)
            u=Sp[u];
        else
            u=_Sp[d*i];
        ans+=u;
    }
    _Sp[d]=clat_S(n)-ans;
}
void slove(uint n)
{
    if(n<LIM)return;
    for(uint i=n/LIM;i;i--)clat_Sp(n/i,i);
}
uint gcd(uint a,uint b)
{
    if(b==0)return a;
    return gcd(b,a%b);
}
uint DFS(uint t,uint n)
{
    if(t==0)return n>0?clat_Sk(n):0;
    return DFS(t-1,n)-DFS(t-1,n/prim[t])*prim[t];
}
int main ()
{
    uint n,K,m,lim,X,fir;
    scanf("%u %u",&n,&K); size=K+1;
    m=sqrt(n)+1;    lim=pow(n,0.25)+1;  X=m;
    bit[0].C();
    for(uint i=1;i<33;i++) bit[i]=bit[i-1]*bit[i-1];
    for(uint i=1;i<MAXN;i++)Ik[i]=Pow(i,K);
    for(uint i=1;i<=m;i++)
    {
        b[i]=i;
        B[i]=n/i;
        c[i]=c[i-1]+Ik[i];
        C[i]=clat_Sk(n/i);
    }
    for(uint i=2;i<MAXN;i++)
    {
        pai[i]=pai[i-1];
        paik[i]=paik[i-1];
        if(vis[i])continue;
        prim[++pai[i]]=i;
        paik[i]+=Ik[i];
        for(uint j=i<<1;j<MAXN;j+=i)   vis[j]=true;
    }
    uint k=1;
    while(prim[k]<lim)
    {
        for(uint i=1;i<m;i++)
        {
            uint t=prim[k]*i;
            uint u=n/t;
            uint bu;
            if(u<X)
            {
                bu=b[u];
                u=c[u];
            }
            else
            {
                u=C[t];
                bu=B[t];
            }
            B[i]-=bu;
            S[i]+=u;
            C[i]-=u*Ik[prim[k]];
        }
        for(uint i=X-1;i;i--)
        {
            b[i]-=b[i/prim[k]];
            s[i]+=c[i/prim[k]];
            c[i]-=c[i/prim[k]]*Ik[prim[k]];
        }
        k++;
    }
    for(uint i=prim[k-1]+1;i<m;i++)s[i]+=pai[i]-k+1;
    lim=m;
    memset(tmp,0x3f,sizeof tmp);
    while(prim[k]<m)
    {
        fir=lim;
        lim=n/(prim[k]*prim[k])+1;
        for(uint i=1;i<lim;i++)
        {
            uint t=i*prim[k];
            uint u=n/t,bu;
            if(u<X)
            {
                if(u<prim[k])   u=bu=1;
                else
                {
                    bu=pai[u]-k+2;
                    u=paik[u]-paik[prim[k-1]]+1;
                }
            }
            else
            {
                if(tmp[t]!=INF)
                {
                    u=C[t]-paik[prim[k-1]]+paik[prim[tmp[t]-1]];
                    bu=B[t]-k+tmp[t];
                }
                else
                {
                    u=C[t];
                    bu=B[t];
                }
            }
            S[i]+=u;
            B[i]-=bu;
            C[i]-=u*Ik[prim[k]];
        }
        for(uint i=lim;i<fir;i++)tmp[i]=k;
        k++;
    }
    for(uint i=1;i<m;i++)
        if(tmp[i]<INF)
        {
            S[i]+=k-tmp[i]-1;
            B[i]-=k-tmp[i]-1;
        }
    for(uint i=1;i<m;i++) S[i]+=B[i]-1;
    memset(vis,0,sizeof vis);
    for(uint i=2;i<MAXN;i++)
    {
        if(i<100)test[i]=test[i-1]+f[i];
        if(vis[i])continue;
        for(uint j=i,t=1;j<MAXN;j+=i,t++)
        {
            if(vis[j])continue;
            vis[j]=true;
            f[j]=Ik[t];
        }
        if(i<100)test[i]+=f[i];
    }
    for(uint i=1;i<LIM;i++)
    {
        Sp[i]=i-Sp[i];
        for(uint j=i<<1;j<LIM;j+=i)Sp[j]+=Sp[i];
        Sp[i]+=Sp[i-1];
    }
    slove(n);

    uint ans=0;

    for(uint L=1;L<m;L++)
    {
        uint u2;
        if(L+1==m)
            u2=s[n/(L+1)];
        else
            u2=S[L+1];
        ans+=(S[L]-u2)*(Sp[L]*2-1);
    }
    for(uint d=n/m;d>1;d--)
    {
        uint u=n/d;
        if(u<LIM)
            u=Sp[u];
        else
            u=_Sp[d];
        ans+=f[d]*(2*u-1);
    }
    printf("%u\n",ans);
    return 0;
}
/*
 938393691 46
 1106647490
 */