51nod 1222 最小公倍数计数

51nod 1222 最小公倍数计数

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

次题有毒。。。。

计算 $$\sum_{i=1}^n\sum_{j=1}^i\sum_{t=j}^i[lcm(j, t)=i]$$

周阁筛：

非常经典的周阁筛。

$$F(n)=\sum_{a=1}^n\sum_{b=a}^n[lcm(a, b)=n]\\=\sum_{d|n}\sum_{a=1}^n\sum_{b=1}^n[gcd(a ,b)=d][ab=dn]\\=\sum_{d|n}\sum_{ab=\frac{n}{d}}[a\perp b][a\leq b]\\=\frac{1}{2}\Big(1+\sum_{d|n}\sum_{ab=d}[a\perp b]\Big)$$

显然：

$$\sum_{d|n}\sum_{ab=d}[a\perp b]$$

相当于计算。所有之积为$d$的互素$(a,b)$对

即：

$$\sum_{d|n}\sum_{ab=d}[a\perp b]\\=\sum_{a|n}\sum_{b|n}[a\perp b]\\=\tau(n^2)$$

所以：

$$2F(n)=1+\tau(n^2)$$

周阁筛计算。。。。问题在于。周阁筛很容超时。

我们定义$G(i,j)$为$[1,j]$上。不可以被前$i$个素数整除的答案。

经典的周阁筛更新公式：

$$G(i,j)=G(i+1,j)+\sum_{c\geq 1}F(P_{i}^c)G(i+1,\Big\lfloor\frac{j}{P_i^c}\Big\rfloor)$$

经典周阁筛在$\sum_{c\geq 1}$这部分暴力更新。速度比较慢。

这里再次定义定义：

$$A(i,j)=\sum_{c\geq 1}F(P_{i}^c)G(i+1,\Big\lfloor\frac{j}{P_i^c}\Big\rfloor)\\B(i,j)=\sum_{c\geq 1}G(i+1,\Big\lfloor\frac{j}{P_i^c}\Big\rfloor)$$

显然： $$A(i,j)=A(i,\Big\lfloor\frac{j}{P}\Big\rfloor )+2B(i,\Big\lfloor\frac{j}{P}\Big\rfloor )+3G(i,\Big\lfloor\frac{j}{P}\Big\rfloor )$$

这是因为

$$F(P^{k+1})=2+F(P^k)$$

对贡献数组$A[]$辅助数组$B[]$的计算优化类似周阁筛。计算的技巧自行推到。

比经典周阁筛快一倍。。加上我手残。写的好的话更快。。。。。。。

杜教筛方法：

前面推到：

$$2F(n)=1+\sum_{d|n}\sum_{ab=d}[a\perp b]$$

对于：

$$\sum_{d|n}\sum_{ab=d}[a\perp b]\\=\sum_{d|n}\sum_{a|d}\mu^2(a)$$

令$\lambda$为刘维尔函数。$f=\mu^2*I*I=\mu^2*\tau$

显然： $$\sum_{d|n}\mu^2(d)=\sum_{d|n}\mu(d)\lambda(d)$$

因为刘维尔函数是完全积性函数。 则：

$$\big(u^2*\lambda\big)(n)=\sum_{d|n}\mu(d)\lambda(d)\lambda(\frac{n}{d})\\=\lambda(n)\sum_{d|n}\mu(d)=[n=1]$$

又：

$$\lambda(n)=(-1)^{a_1+a_2+a_3..+a_k}$$

其中：

$$n=\prod_{i=1}^kP_i^{a_i}$$

所以：

$$\big(\lambda*I\big)(n)=\sum_{d|n}\lambda(d)\\=\sum_{x_1=0}^{a_1}\sum_{x_2=0}^{a_2}...\sum_{x_k=0}^{a_k}(-1)^{x_1+x_2+...x_k}\\=(x^0+x^1..x^{a_1})...(x^0+x^1..x^{a_k})\Big|_{x=-1}$$

显然。当且仅当。$a_1,a_2,..a_k$都为偶数时。上式不为0

所以： $$h(n)=\big(\lambda*I\big)(n)=[n为完全平方数]$$

所以： $$h=\lambda*I\\f*\lambda=\tau$$

杜教筛计算之。。。（卡常数）

另一种姿态：

$$f(n)=\sum_{a|n}\sum_{b|n}[a\perp b]$$

令

$$g(n,k)=\sum_{a|n}\sum_{b|n}[gcd(a, b)=k]$$

$$G(n,k)=\sum_{k|d}g(n,d)\\=[k|n]\sum_{a|\frac{n}{k}}\sum_{b|\frac{n}{k}}1=\tau^2(\frac{n}{k})[k|n]$$

反演有：

$$f(n)=g(n,1)=\sum_{d>=1}\mu(d)\tau^2(\frac{n}{d})[d|n]\\=\sum_{d|n}\mu(d)\tau^2(\frac{n}{d})$$

所以 $$f*I=\tau^2$$

这也就是说： $$\tau(n^2)=\mu^2*\tau=\mu*\tau^2\\\tau(n^2)*\lambda=e$$

姿态很多。。。使用： $$f*I=\tau^2$$ 这个姿态还是快很多的。只需要计算一遍。（不过没写。

下面给出周阁筛的实现：
5156ms 卡过

#include <algorithm>
#include <string.h>
#include <stdio.h>
#include <cmath>
#define MAXN 340000
using namespace std;
typedef long long LL;
const int INF=0x3f3f3f3f;
int prim[MAXN];
int pa0[MAXN];
bool vis[MAXN];
int tmp[MAXN];
int C0[MAXN];
LL _C0[MAXN];
int G[MAXN];
LL _G[MAXN];
int A[MAXN];
LL _A[MAXN];
int B[MAXN];
LL _B[MAXN];
int Sf[MAXN];
LL slove(LL n)
{
    memset(vis,0,sizeof vis);
    if(n==1)return 1;
    if(n==0)return 0;
    int m=sqrt(n)+1,V=pow(n,0.25)+1;
    for(int i=1;i<=m;i++)
    {
        C0[i]=i;
        _C0[i]=n/i;
    }
    int k=1;
    while(prim[k]<V)
    {
        int P=prim[k];
        for(int i=1;i<m;i++)
        {
            LL t=(LL)i*P;
            LL u=n/t;
            if(u<m)
                _C0[i]-=C0[u];
            else
                _C0[i]-=_C0[t];
        }
        for(int i=m-1;i;i--)    C0[i]-=C0[i/P];
        k++;
    }
    memset(tmp,0x3f,sizeof tmp);
    int fir,lim=m;
    while(prim[k]<m)
    {
        LL P=prim[k];
        fir=lim;
        lim=(int)(1+n/(P*P));
        for(int i=1;i<lim;i++)
        {
            LL t=P*i;
            LL u=n/t;
            if(u<m)
            {
                if(u<P)
                    _C0[i]-=1;
                else
                    _C0[i]-=pa0[u]-k+2;
            }
            else
            {
                if(tmp[t]<INF)
                    _C0[i]-=_C0[t]-(k-tmp[t]);
                else
                    _C0[i]-=_C0[t];
            }
        }
        for(int i=lim;i<fir;i++)tmp[i]=k;
        k++;
    }
    for(int i=1;i<m;i++) if(tmp[i]!=INF)_C0[i]-=k-tmp[i];
    for(int i=1;i<m;i++) _G[i]=_C0[i]*3-2;
    k--;    V--;
    int X=prim[k];
    memset(tmp,0,sizeof tmp);
    while(prim[k]>V)
    {
        LL P=prim[k];
        lim=(int)(n/(P*P)+1);
        for(int i=lim-1;i;i--)
        {
            LL t=i*P;
            if(t>=lim)
            {
                LL u=n/t;
                if(u<P)
                {
                    _B[i]=0;
                    _A[i]=0;
                }
                else
                {
                    _B[i]=1;
                    _A[i]=5;
                }
            }
            else
            {
                _B[i]=_B[t];
                _A[i]=_A[t]+(_B[t]<<1);
            }
            LL u=n/t;
            if(u<m)
            {
                if(u<P)
                {
                    _B[i]++;
                    _A[i]+=3;
                }
                else
                {
                    _B[i]+=Sf[u]-Sf[prim[k]]+1;
                    _A[i]+=(Sf[u]-Sf[prim[k]]+1)*3;
                }
            }
            else
            {
                if(!vis[t])
                {
                    _B[i]+=(_G[t]+Sf[X]-Sf[prim[k]]);
                    _A[i]+=(_G[t]+Sf[X]-Sf[prim[k]])*3;
                }
                else
                {
                    _B[i]+=_G[t];
                    _A[i]+=_G[t]*3;
                }
            }
        }
        for(int i=1;i<lim;i++)
        {
            _G[i]+=_A[i];
            if(!vis[i])
            {
                _G[i]+=Sf[X]-Sf[prim[k]];
                vis[i]=true;
            }
        }
        k--;
    }

    for(int i=lim;i<m;i++) _G[i]+=Sf[X]-Sf[prim[k]];

    for(int i=1;i<m;i++)
    {
        if(i<=V)
            G[i]=1;
        else
            G[i]=Sf[i]-Sf[prim[k]]+1;
    }
    memset(A,0,sizeof A);
    memset(_A,0,sizeof _A);
    memset(B,0,sizeof _B);
    memset(_B,0,sizeof _B);
    while(k)
    {
        LL P=prim[k];
        for(int i=(int)P;i<m;i++)
        {
            int d=i/P;
            B[i]=B[d]+G[d];
            A[i]=A[d]+(B[d]<<1)+3*G[d];
        }
        for(int i=m-1;i;i--)
        {
            LL t=i*P;
            LL u=n/t;
            if(u<m)
            {
                _B[i]=B[u]+G[u];
                _A[i]=A[u]+(B[u]<<1)+3*G[u];
            }
            else
            {
                _B[i]=_B[t]+_G[t];
                _A[i]=_A[t]+(_B[t]<<1)+3*_G[t];
            }
        }
        for(int i=1;i<m;i++)    _G[i]+=_A[i];
        for(int i=m-1;i;i--)    G[i]+=A[i];
        k--;
    }
    return (_G[1]+n)>>1;
}

int main ()
{
    for(int i=2;i<MAXN;i++)
    {
        pa0[i]=pa0[i-1];
        Sf[i]=Sf[i-1];
        if(vis[i])continue;
        for(int j=i<<1;j<MAXN;j+=i)vis[j]=true;
        prim[++pa0[i]]=i;
        Sf[i]+=3;
    }
    LL a,b;
    scanf("%lld %lld",&a,&b);
    printf("%lld\n",slove(b)-slove(a-1));
    return 0;
}