51nod 1222 最小公倍数计数

版权声明:这里要写版权声明么。。。 https://blog.csdn.net/ZLH_HHHH/article/details/78406585

51nod 1222 最小公倍数计数

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

次题有毒。。。。

计算
i=1nj=1it=ji[lcm(j,t)=i]

周阁筛:

非常经典的周阁筛。

F(n)=a=1nb=an[lcm(a,b)=n]=d|na=1nb=1n[gcd(a,b)=d][ab=dn]=d|nab=nd[ab][ab]=12(1+d|nab=d[ab])

显然:

d|nab=d[ab]

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

即:

d|nab=d[ab]=a|nb|n[ab]=τ(n2)

所以:

2F(n)=1+τ(n2)

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

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

经典的周阁筛更新公式:

G(i,j)=G(i+1,j)+c1F(Pci)G(i+1,jPci)

经典周阁筛在c1这部分暴力更新。速度比较慢。

这里再次定义定义:

A(i,j)=c1F(Pci)G(i+1,jPci)B(i,j)=c1G(i+1,jPci)

显然:
A(i,j)=A(i,jP)+2B(i,jP)+3G(i,jP)

这是因为

F(Pk+1)=2+F(Pk)

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

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

杜教筛方法:

前面推到:

2F(n)=1+d|nab=d[ab]

对于:

d|nab=d[ab]=d|na|dμ2(a)

λ为刘维尔函数。f=μ2II=μ2τ

显然:
d|nμ2(d)=d|nμ(d)λ(d)

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

(u2λ)(n)=d|nμ(d)λ(d)λ(nd)=λ(n)d|nμ(d)=[n=1]

又:

λ(n)=(1)a1+a2+a3..+ak

其中:

n=i=1kPaii

所以:

(λI)(n)=d|nλ(d)=x1=0a1x2=0a2...xk=0ak(1)x1+x2+...xk=(x0+x1..xa1)...(x0+x1..xak)x=1

显然。当且仅当。a1,a2,..ak都为偶数时。上式不为0

所以:
h(n)=(λI)(n)=[n]

所以:
h=λIfλ=τ

杜教筛计算之。。。(卡常数)

另一种姿态:

f(n)=a|nb|n[ab]

g(n,k)=a|nb|n[gcd(a,b)=k]

G(n,k)=k|dg(n,d)=[k|n]a|nkb|nk1=τ2(nk)[k|n]

反演有:

f(n)=g(n,1)=d>=1μ(d)τ2(nd)[d|n]=d|nμ(d)τ2(nd)

所以
fI=τ2

这也就是说:
τ(n2)=μ2τ=μτ2τ(n2)λ=e

姿态很多。。。使用:
fI=τ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;
}

没有更多推荐了,返回首页