bzoj 4659: Lcm （反演）

题目描述

传送门

题目大意：

1<=a<=A

1<=b<=B

∀n>1n†gcd(a,b)

求满足条件的所有a,b的lcm(a,b)的和。

题解

设 f(x) 表示x是否含有平方因子，如果不含平方因子， f(x)=1 ，否则 f(x)=0

那么当a,b互质的时候, f(ab)=f(a)∗f(b) 那么 f(x) 是积性函数

然后我们用式子表示上面的要求

∑p=1n∑i=1n∑j=1mi∗jp∗f(p)[gcd(i,j)=p]

∑p=1nf(p)∗p∑i=1⌊np⌋∑j=1⌊mp⌋i∗j∗[gcd(i,j)=1]

∑p=1nf(p)∗p∑i=1⌊np⌋∑j=1⌊mp⌋i∗j∗∑d|(i,j)μ(d)

∑p=1nf(p)∗p∑d=1nμ(d)∑i=1⌊np⌋[d|i]∗i∑j=1⌊mp⌋[d|j]∗j

∑p=1nf(p)∗p∑d=1nμ(d)∑i=1⌊np∗d⌋d∗i∑j=1⌊mp∗d⌋d∗j

∑p=1nf(p)∗p∑d=1nμ(d)∗d2sum(⌊np∗d⌋,⌊mp∗d⌋)

其中 sum(x,y)=x(x+1)2∗y(y+1)2 ,令 T=p∗d

∑T=1nsum(⌊nT⌋,⌊mT⌋)∑d|Tμ(Td)∗(Td)2∗d∗f(d)

设 h(d)=∑d|Tμ(Td)∗(Td)2∗d∗f(d) ，分析可知 h(d) 是积性函数，那么我们就可以线性筛了

当 x 是质数的时候,h(x)=μ(1)∗1∗x∗f(x)+μ(x)∗x2∗f(x)=x−x∗x

互质的时候可以直接相乘，关键是不互质的时候，我们考虑质因数分解

x=p1q1∗....∗pkqk ,计算一下可以发现，如果 qi>3 那么 f(piqi)=0

那么假设当前加入的质因子就是 p1 ，如果 q1>=2 ，那么 f(p1∗x)=0

如果 q1=1 ，那么我们把 p12 ,原式就变成了 h(x∗p1)=h(p12)∗h(p2q2∗...∗pkqk)

这两项是互质的，根据积性函数的性质可以直接相乘。 h(p12)=μ(p1)∗p12∗p1∗f(p1)=−p13

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 4000003
using namespace std;
const int p=<<;
int n,m,T,prime[N+],pd[N+],mu[N+],cnt[N+],g[N+];
int f[N+];
void init()
{
    mu[]=; f[]=;
    for (int i=;i<=N;i++) {
        if (!pd[i]) {
            prime[++prime[]]=i;
            mu[i]=-; cnt[i]=; g[i]=i; f[i]=i-i*i;
        }
        for (int j=;j<=prime[];j++) {
            int t=i*prime[j];
            if (t>N) break;
            pd[t]=;
            if (i%prime[j]==) {
                mu[t]=;
                cnt[t]=cnt[i]+; g[t]=g[i]; int pi=prime[j];
                int temp=i/g[i];
                if (cnt[t]==) f[t]=-pi*pi*pi*f[temp];
                else f[t]=;
                break;
            }
            mu[t]=-mu[i]; cnt[t]=; g[t]=prime[j]; 
            f[t]=f[i]*f[prime[j]];
        }
    }
    for (int i=;i<=N;i++) f[i]=f[i]+f[i-];
}
int calc(int x,int y)
{
    int t1=(x+)*x>>; int t2=(y+)*y>>;
    return t1*t2;
}
int main()
{
    freopen("a.in","r",stdin);
    //freopen("my.out","w",stdout);
    scanf("%d",&T); init();
    while (T--)
    {
      scanf("%d%d",&n,&m);
      if (n>m) swap(n,m);
      int ans=;
      for (int i=,j;i<=n;i=j+) {
        j=min(n/(n/i),m/(m/i));
        ans=ans+(f[j]-f[i-])*calc(n/i,m/i);
      }
      printf("%d\n",(ans%p+p)%p);
    }
}