bzoj4659: Lcm

链接

　　http://www.lydsy.com/JudgeOnline/problem.php?id=4659

题解

　　

ans=∑i=1n∑j=1mlcm(i,j)[gcd(i,j)无平方因子]　　　　=∑i=1n∑j=1mlcm(i,j)μ(gcd(i,j))2

　　这里用μ(gcd(i,j))2意思就是，当gcd(i,j)有平方因子的时候就乘1，否则就乘0。
　　令g(x)=(x)(x+1)2
　　你就各种推啊推，最后
　　

ans=∑x=1ng(⌊nx⌋)g(⌊mx⌋)x∑d|xmu(xd)xdμ(d)2

　　然后令f(x)=x∑d|xmu(xd)xdμ(d)2
　　O(nlog2n)预处理，O(TN−−√)处理询问。
　　复杂度太高了，对f()打个表，发现那些只有一个质因子的数的函数值在次数>3的时候就等于0了，对这些数特殊处理之后就可以线性筛
　　复杂度降到O(n)

代码

//线性筛+莫比乌斯反演（nlogn） #include <cstdio>#include <algorithm>#define maxn 5000001#define lim 1152921504606846976ll#define mod 1073741824ll#define ll long longusing namespace std;ll N, f[maxn];int prime[maxn], mu[maxn], h[maxn], g[maxn], q[maxn][2];bool mark[maxn];void init(){    int i, j;    mu[1]=1;    for(i=2;i<=N;i++)    {        if(!mark[i])prime[++prime[0]]=i, mu[i]=-1;        for(j=1;j<=prime[0] and i*prime[j]<=N;j++)        {            mark[i*prime[j]]=1;            if(i%prime[j]==0){mu[i*prime[j]]=0;break;}            mu[i*prime[j]]=-mu[i];        }    }    for(i=1;i<=N;i++)h[i]=mu[i]*i, mu[i]=mu[i]>0?mu[i]:-mu[i];    for(i=1;i<=N;i++)for(j=i;j<=N;j+=i)    {        f[j]+=h[j/i]*mu[i];        if(f[j]>lim or -f[j]>lim)f[j]%=mod;    }    for(i=1;i<=N;i++)f[i]=(f[i-1]+i*f[i])%mod;    for(i=1;i<=N;i++)g[i]=((ll)i*(i+1)>>1)%mod;}void work(ll n, ll m){    ll i,  last, ans=0;    if(n>m)swap(n,m);    for(i=1;i<=n;i=last+1)    {        last=min(n/(n/i),m/(m/i));        ans+=g[n/i]*g[m/i]%mod*(f[last]-f[i-1]);        if(ans>lim)ans%=mod;    }    printf("%lld\n",(ans%mod+mod)%mod);}inline ll read(ll x=0){    char c=getchar();    while(c<48 and c>57)c=getchar();    while(c>=48 and c<=57)x=(x<<1)+(x<<3)+c-48, c=getchar();    return x;}int main(){    ll T, i;    for(T=read(),i=1;i<=T;i++)    {        q[i][0]=read(), q[i][1]=read();        N=max((int)N,max(q[i][0],q[i][1]));    }    init();    for(i=1;i<=T;i++)work(q[i][0],q[i][1]);    return 0;}

//线性筛+莫比乌斯反演 O(N) #include <cstdio>#include <algorithm>#include <cmath>#define maxn 5000001#define lim 1152921504606846976ll#define mod 1073741824ll#define ll long longusing namespace std;ll N, f[maxn];int prime[maxn], mu[maxn], g[maxn], q[maxn][2], pr[maxn];bool mark[maxn];void init(){    int i, j, t;    mu[1]=1; f[1]=1;    for(i=2;i<=N;i++)    {        if(!mark[i])prime[++prime[0]]=i, f[i]=-i+1, pr[i]=i;        else        {            if(i==pr[i])            {                t=int(sqrt(i));                if(!mark[t] and t*t==i)f[i]=-(int)sqrt(i);            }            else f[i]=f[i/pr[i]]*f[pr[i]];        }        for(j=1;j<=prime[0] and i*prime[j]<=N;j++)        {            mark[i*prime[j]]=1, pr[i*prime[j]]=prime[j];            if(i%prime[j]==0){pr[i*prime[j]]=pr[i]*prime[j];break;}        }    }    for(i=1;i<=N;i++)f[i]=(f[i-1]+i*f[i])%mod;    for(i=1;i<=N;i++)g[i]=((ll)i*(i+1)>>1)%mod;}void work(ll n, ll m){    ll i,  last, ans=0;    if(n>m)swap(n,m);    for(i=1;i<=n;i=last+1)    {        last=min(n/(n/i),m/(m/i));        ans+=g[n/i]*g[m/i]%mod*(f[last]-f[i-1]);        if(ans>lim)ans%=mod;    }    printf("%lld\n",(ans%mod+mod)%mod);}inline ll read(ll x=0){    char c=getchar();    while(c<48 and c>57)c=getchar();    while(c>=48 and c<=57)x=(x<<1)+(x<<3)+c-48, c=getchar();    return x;}int main(){    ll T, i;    for(T=read(),i=1;i<=T;i++)    {        q[i][0]=read(), q[i][1]=read();        N=max((int)N,max(q[i][0],q[i][1]));    }    init();    for(i=1;i<=T;i++)work(q[i][0],q[i][1]);    return 0;}