51nod 1575 Gcd and Lcm

51nod 1575 Gcd and Lcm

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

这个题目。完美的体现了洲阁筛更为通用。。。。

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

令 $F(n)$有:

$$F(n)=\sum_{i=1}^n\sum_{j=1}^nlcm(gcd(i,n),gcd(j,n))\\=\sum_{a|n}\sum_{b|n}\sum_{i=1}^n\sum_{j=1}^n\frac{ab}{gcd(a,b)}[gcd(i,n)=a][gcd(j,n)=b]\\=\sum_{a|n}\sum_{b|n}\frac{ab}{gcd(a,b)}\sum_{i=1}^{\frac{n}{a}}\sum_{j=1}^{\frac n b}[i\perp n][j\perp n]\\=\sum_{a|n}\sum_{b|n}\frac{ab}{gcd(a,b)}\varphi(\frac{n}{a})\varphi(\frac n b)\\=\sum_{d|n}\frac 1 d\sum_{a|n}\sum_{b|n}ab\varphi(\frac{n}{a})\varphi(\frac n b)[gcd(a,b)=d]\\=\sum_{d|n}d\sum_{a\big|\frac n d}\sum_{b\big|\frac n d}ab\varphi(\frac{n}{ad})\varphi(\frac{n}{bd})[a\perp b]\\=\sum_{d|n}\frac n d\sum_{a|d}\sum_{b|d}ab\varphi(\frac d a)\varphi(\frac d b)[a\perp b]$$

令： $$g(n)=g(n,1)\\g(n,k)=\sum_{a|n}\sum_{b|n}ab\varphi(\frac n a)\varphi(\frac n b)[gcd(a,b)=k]$$

令： $$G(n,k)=[k|n]\sum_{k|d}g(n,d)\\=[k|n]k^2\sum_{a\big|\frac n k}\sum_{b\big|\frac n k}ab\varphi(\frac{n}{ka})\varphi(\frac{n}{kb})$$

反演有： $$g(n)=g(n,1)=\sum_{d\geq 1}\mu(d)d^2\sum_{a\big|\frac n d}\sum_{b\big|\frac n d}ab\varphi(\frac{n}{da})\varphi(\frac{n}{db})[d|n]\\=\sum_{d|n}\mu(\frac{n}{d})\big(\frac{n}{d}\big)^2\sum_{a|d}\sum_{b|d}ab\varphi(\frac d a)\varphi(\frac d b)\\=\sum_{d|n}\mu(\frac{n}{d})\big(\frac{n}{d}\big)^2\Big(\sum_{a|d}a\varphi(\frac d a)\Big)^2$$

因为积性函数的狄利克雷积依然为积性函数。所以$g(n)$为积性函数。

因为$F=N*g$

所以$F$为积性函数： $$F(P^k)=\sum_{d|P^k}\frac {P^k} d\sum_{a|d}\sum_{b|d}ab\varphi(\frac d a)\varphi(\frac d b)[a\perp b]\\=\sum_{x=0}^kP^{k-x}\sum_{i=0}^x\sum_{j=0}^xP^{i+j}\varphi(P^{x-i})\varphi(P^{x-j})[P^i\perp P^j]\\=\sum_{x=0}^kP^{k-x}\Big(2\varphi(P^x)\sum_{j=0}^{x-1}P^j\varphi(P^{x-j})-\varphi^2(P^x)+2P^x\varphi(P^x)\Big)\\=\sum_{x=0}^kP^{k-x}\Big(2x\varphi^2(P^x)-\varphi^2(P^x)+2P^x\varphi(P^x)\Big)\\=P^k+\sum_{x=1}^kP^{k-x}\Big((2x-1)(P^x-P^{x-1})^2+2P^x(P^x-P^{x-1})\Big)\\=P^k+(P^k-P^{k-1})\sum_{x=1}^k\Big((2x-1)(P^x-P^{x-1})+2P^x\Big)$$

其中： $$\sum_{x=1}^k(2x-1)(P^x-P^{x-1})\\=\sum_{x=1}^k(2x-1)P^x-\sum_{x=1}^k(2x-1)P^{x-1}\\=\sum_{x=1}^k(2x-1)P^x-\sum_{x=0}^{k-1}(2x+1)P^x\\=-2\sum_{x=1}^{k-1}P^x-1+(2k-1)P^k$$

所以： $$P^k+(P^k-P^{k-1})\sum_{x=1}^k\Big((2x-1)(P^x-P^{x-1})+2P^x\Big)\\=P^k+(P^k-P^{k-1})(-2\sum_{x=1}^{k-1}P^x-1+(2k-1)P^k+2\sum_{x=1}^kP^x)\\=P^k+(P^k-P^{k-1})(2P^k+(2k-1)P^k-1)\\=2(k+1)(P^{2k}-P^{2k-1})+P^{k-1}$$

所以： $$F(P^k)=(2k+1)(P^{2k}-P^{2k-1})+P^{k-1}$$

特别的： $$F(1)=1$$

洲阁筛计算之

#include <algorithm>
#include <string.h>
#include <stdio.h>
#include <cmath>
#define MAXN 1000000
#define SQR 40000
using namespace std;
typedef unsigned int uint;
const uint INF=0x3f3f3f3f;

uint C0[SQR];
uint _C0[SQR];
uint C1[SQR];
uint _C1[SQR];
uint C2[SQR];
uint _C2[SQR];
uint prim[SQR];
uint pa0[SQR];
uint pa1[SQR];
uint pa2[SQR];
uint F[SQR];
uint _F[SQR];
uint tmp[SQR];
uint Sf[SQR];
uint pf[SQR][32];
bool vis[SQR];
uint S1(uint n)
{
    if(n&1)return n*((n+1)>>1);
    return (n>>1)*(n+1);
}
uint S2(uint n)
{
    uint n1=n+1,n2=2*n+1;
    if(n&1)n1>>=1;
    else n>>=1;

    if(n%3==0)n/=3;
    else
    {
        if(n2%3==0)
            n2/=3;
        else
            n1/=3;
    }
    return n*n1*n2;
}
uint Pow(uint a,uint b)
{
    uint tmp=1;
    while(b)
    {
        if(b&1)
            tmp*=a;
        a*=a;
        b>>=1;
    }
    return tmp;
}
uint f(uint P,uint c)
{
    return (2*c+1)*(Pow(P,c<<1)-Pow(P,(c<<1)-1))+Pow(P,c-1);
}
uint gcd(uint a,uint b)
{
    if(b==0)return a;
    return gcd(b,a%b);
}
uint FF(uint a,uint b,uint c)
{
    uint A=gcd(a,b);
    uint B=gcd(a,c);
    return A*B/gcd(A,B);
}
int main ()
{
    for(uint i=2;i<SQR;i++)
    {
        pa0[i]=pa0[i-1];
        pa1[i]=pa1[i-1];
        pa2[i]=pa2[i-1];
        Sf[i]=Sf[i-1];
        if(vis[i])continue;
        prim[++pa0[i]]=i;
        pa1[i]+=i;
        pa2[i]+=i*i;
        Sf[i]+=f(i,1);
        for(uint j=1;j<31;j++) pf[i][j]=f(i,j);
        for(uint j=i<<1;j<SQR;j+=i)vis[j]=true;
    }
    uint n,T;
    scanf("%u",&T);
    while(T--)
    {
        scanf("%u",&n);
        uint m=sqrt(n)+1,V=pow(n,0.25)+1;
        for(uint i=1;i<m;i++)
        {
            C0[i]=i;
            _C0[i]=n/i;
            C1[i]=C1[i-1]+i;
            _C1[i]=S1(n/i);
            C2[i]=C2[i-1]+i*i;
            _C2[i]=S2(n/i);
        }
        uint k=1;
        while(prim[k]<V)
        {
            uint P=prim[k];
            for(uint i=1;i<m;i++)
            {
                uint t=i*prim[k];
                uint u=n/t;
                if(u<m)
                {
                    _C0[i]-=C0[u];
                    _C1[i]-=C1[u]*P;
                    _C2[i]-=C2[u]*P*P;
                }
                else
                {
                    _C0[i]-=_C0[t];
                    _C1[i]-=_C1[t]*P;
                    _C2[i]-=_C2[t]*P*P;
                }
            }
            for(uint i=m-1;i;i--)
            {
                C0[i]-=C0[i/P];
                C1[i]-=C1[i/P]*P;
                C2[i]-=C2[i/P]*P*P;
            }
            k++;
        }
        memset(tmp,0x3f,sizeof tmp);
        uint fir,lim=m;
        while(prim[k]<m)
        {
            uint P=prim[k];
            fir=lim;
            lim=1+n/(prim[k]*prim[k]);
            for(uint i=1;i<lim;i++)
            {
                uint t=P*i;
                uint u=n/t;
                if(u<m)
                {
                   if(u<prim[k])
                   {
                       _C0[i]-=1;
                       _C1[i]-=P;
                       _C2[i]-=P*P;
                   }
                   else
                   {
                       _C0[i]-=pa0[u]-k+2;
                       _C1[i]-=(pa1[u]-pa1[prim[k-1]]+1)*P;
                       _C2[i]-=(pa2[u]-pa2[prim[k-1]]+1)*P*P;
                   }
                }
                else
                {
                    if(tmp[t]<INF)
                    {
                        _C0[i]-=_C0[t]-(k-tmp[t]);
                        _C1[i]-=P*(_C1[t]-(pa1[prim[k-1]]-pa1[prim[tmp[t]-1]]));
                        _C2[i]-=P*P*(_C2[t]-(pa2[prim[k-1]]-pa2[prim[tmp[t]-1]]));
                    }
                    else
                    {
                        _C0[i]-=_C0[t];
                        _C1[i]-=_C1[t]*P;
                        _C2[i]-=_C2[t]*P*P;
                    }
                }
            }
            for(uint i=lim;i<fir;i++)tmp[i]=k;
            k++;
        }
        for(uint i=1;i<m;i++)
        {
            if(tmp[i]==INF)continue;
            _C0[i]-=k-tmp[i];
            _C1[i]-=pa1[prim[k-1]]-pa1[prim[tmp[i]-1]];
            _C2[i]-=pa2[prim[k-1]]-pa2[prim[tmp[i]-1]];
        }
        for(uint i=1;i<m;i++)
            _F[i]=3*(_C2[i]-_C1[i])+_C0[i];
        k--;
        V--;
        uint X=prim[k];
        memset(vis,0,sizeof vis);
        while(prim[k]>V)
        {
            uint P=prim[k];
            uint bf;
            lim=n/(P*P)+1;
            for(uint i=1;i<lim;i++)
            {
                uint t=i,c=0;
                while(t<=n/P)
                {
                    c++;
                    t*=P;
                    bf=pf[P][c];
                    uint u=n/t;
                    if(u<m)
                    {
                        if(u<prim[k])
                            _F[i]+=bf;
                        else
                            _F[i]+=(Sf[u]-Sf[prim[k]]+1)*bf;
                    }
                    else
                    {
                        if(!vis[t])
                            _F[i]+=bf*(_F[t]+Sf[X]-Sf[prim[k]]);
                        else
                            _F[i]+=bf*_F[t];
                    }
                }
                if(!vis[i])
                {
                    _F[i]+=Sf[X]-Sf[prim[k]];
                    vis[i]=true;
                }
            }
            k--;
        }

        for(uint i=lim;i<m;i++)_F[i]+=Sf[X]-Sf[prim[k]];
        for(uint i=1;i<m;i++)
        {
            if(i<=V)
                F[i]=1;
            else
                F[i]=Sf[i]-Sf[prim[k]]+1;
        }
        while(k)
        {
            uint P=prim[k];
            for(uint i=1;i<m;i++)
            {
                uint d=i,c=0;
                while(d<=n/P)
                {
                    d*=P;
                    c++;
                    uint u=n/d;
                    if(u<m)
                        _F[i]+=F[u]*pf[P][c];
                    else
                        _F[i]+=_F[d]*pf[P][c];
                }
            }

            for(uint i=m-1;i;i--)
            {
                uint d=1,c=0;
                while(d*P<=i)
                {
                    d*=P;
                    c++;
                    F[i]+=F[i/d]*pf[P][c];
                }
            }
            k--;
        }
        printf("%u\n",_F[1]);
    }
    return 0;
}