bzoj4659: Lcm

链接

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

题解

  

ans=i=1nj=1mlcm(i,j)[gcd(i,j)]    =i=1nj=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)xd|xmu(xd)xdμ(d)2

  然后令 f(x)=xd|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 long
using 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 long
using 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;
}
  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值