BZOJ-4659: Lcm

Description

给出A，B，考虑所有满足l<=a<=A，l<=b<=B，且不存在n>1使得n^2同时整除a和b的有序数对(a，b)，求其lcm(a，b)之和。答案模2^30。
A,B<=4*10^6 多组数据：T<=2000

就是求$\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)*|\mu(gcd(i,j))|$
=$\sum\limits_{i=1}^n\sum\limits_{j=1}^m \frac{i*j}{gcd(i,j)}*|\mu(gcd(i,j))|$
枚举gcd:
=$\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}i*j*d*|\mu(d)|[gcd(i,j)==1]$
=$\sum\limits_{d=1}^nd*|\mu(d)|*\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}i*j*\sum\limits_{k|i,k|j}\mu(k)$
令f(x)=x(x+1)/2
=$\sum\limits_{d=1}^nd*|\mu(d)|*\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)*k^2*f(\lfloor\frac{n}{kd}\rfloor)*f(\lfloor\frac{m}{kd}\rfloor)$
=$\sum\limits_{D=1}^n\sum\limits_{d|D}d*|\mu(d)|*\mu(\frac{D}{d})*(\frac{D}{d})^2*f(\lfloor\frac{n}{D}\rfloor)*f(\lfloor\frac{m}{D}\rfloor)$
=$\sum\limits_{D=1}^nf(\lfloor\frac{n}{D}\rfloor*f(\lfloor\frac{m}{D}\rfloor)*D\sum\limits_{d|D}d*\mu(d)*|\mu(\frac{D}{d})|$
发现后面$\sum\limits_{d|D}d*\mu(d)*|\mu(\frac{D}{d})|$是奇性函数$x*\mu(x)$和$|\mu(x)|$的狄利克雷卷积，可以快筛。之后乘以D，求前缀和。对于每次询问，将n/D按照D分块即可。时间复杂度O(n+T√n)

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define U unsigned 
U int f[4001000],pri[701000],cnt;
bool v[4001000];
void initialize()
{
    f[1]=1;
    for(U int i=2;i<=4000000;++i)
    {
        if(!v[i]) pri[++cnt]=i,f[i]=1-i;
        for(U int j=1;i*pri[j]<=4000000&&j<=cnt;++j)
        {
            v[i*pri[j]]=1;
            if(i%pri[j]) f[i*pri[j]]=f[i]*f[pri[j]];
            else
            {
                if(i%(pri[j]*pri[j])==0) f[i*pri[j]]=0;
                else f[i*pri[j]]=-f[i/pri[j]]*pri[j];
                break;
            }
        }
    }
    for(U int i=1;i<=4000000;++i) f[i]=f[i]*i+f[i-1];
}
inline U int Min(U int x,U int y){return x<y ? x:y;}
inline U int calc(U int x){return x*(x+1)/2;}
int main()
{
    initialize();
    int t; scanf("%d",&t);
    U int n,m,ans,b,e;
    while(t--)
    {
        scanf("%u%u",&n,&m); if(n>m) swap(n,m);
        b=1; ans=0;
        while(b<=n)
        {
            e=Min(n/(n/b),m/(m/b));
            ans+=calc(n/b)*calc(m/b)*(f[e]-f[b-1]);
            b=e+1;
        }
        printf("%u\n",ans&0x3fffffff);
    }
    return 0;
}