bzoj 4659: Lcm （反演）

题目描述

传送门

题目大意：
$1<=a<=A$
$1<=b<=B$
$\forall n>1 n\dagger 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)$是积性函数
然后我们用式子表示上面的要求

$$\sum_{p=1}^n \sum\limits_{i=1}^n\sum\limits_{j=1}^m {i*j\over p}*f(p)[gcd(i,j)=p]$$
$$\sum_{p=1}^n f(p)*p \sum\limits_{i=1}^{\lfloor {n\over p}\rfloor}\sum\limits_{j=1}^{\lfloor{m\over p}\rfloor } i*j*[gcd(i,j)=1]$$
$$\sum_{p=1}^n f(p)*p \sum\limits_{i=1}^{\lfloor {n\over p}\rfloor}\sum\limits_{j=1}^{\lfloor{m\over p}\rfloor } i*j*\sum\limits_{d|(i,j)}\mu(d)$$
$$\sum_{p=1}^n f(p)*p \sum\limits_{d=1}^n \mu(d)\sum\limits_{i=1}^{\lfloor {n\over p}\rfloor}[d|i]*i\sum\limits_{j=1}^{\lfloor{m\over p}\rfloor } [d|j]*j$$
$$\sum_{p=1}^n f(p)*p \sum\limits_{d=1}^n \mu(d)\sum\limits_{i=1}^{\lfloor {n\over p*d}\rfloor}d*i\sum\limits_{j=1}^{\lfloor{m\over p*d}\rfloor } d*j$$
$$\sum_{p=1}^n f(p)*p \sum\limits_{d=1}^n \mu(d)*d^2 sum({\lfloor {n\over p*d}\rfloor},{\lfloor{m\over p*d}\rfloor } )$$
其中 $sum(x,y)={x(x+1)\over 2}*{y(y+1)\over 2}$ ,令 $T=p*d$
$$\sum_{T=1}^n sum({\lfloor {n\over T}\rfloor},{\lfloor{m\over T}\rfloor } )\sum\limits_{d|T} \mu({T\over d})*({T\over d})^2 *d*f(d)$$
设 $h(d)=\sum\limits_{d|T} \mu({T\over d})*({T\over d})^2 *d*f(d)$，分析可知 $h(d)$是积性函数，那么我们就可以线性筛了
当 $x$是质数的时候,$h(x)=\mu(1)*1*x*f(x)+\mu(x)*x^2*f(x)=x-x*x$
互质的时候可以直接相乘，关键是不互质的时候，我们考虑质因数分解
$x=p1^{q1}*....*pk^{qk}$,计算一下可以发现，如果 $qi>3$那么 $f(pi^{qi})=0$
那么假设当前加入的质因子就是 $p1$，如果 $q1>=2$，那么 $f(p1*x)=0$
如果 $q1=1$ ，那么我们把 $p1^2$,原式就变成了 $h(x*p1)=h(p1^2)*h(p2^{q2}*...*pk^{qk})$
这两项是互质的，根据积性函数的性质可以直接相乘。 $h(p1^2)=\mu(p1)*p1^2*p1*f(p1)=-p1^3$

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 4000003
using namespace std;
const int p=1<<30;
int n,m,T,prime[N+3],pd[N+3],mu[N+3],cnt[N+3],g[N+3];
int f[N+3];
void init()
{
    mu[1]=1; f[1]=1;
    for (int i=2;i<=N;i++) {
        if (!pd[i]) {
            prime[++prime[0]]=i;
            mu[i]=-1; cnt[i]=1; g[i]=i; f[i]=i-i*i;
        }
        for (int j=1;j<=prime[0];j++) {
            int t=i*prime[j];
            if (t>N) break;
            pd[t]=1;
            if (i%prime[j]==0) {
                mu[t]=0;
                cnt[t]=cnt[i]+1; g[t]=g[i]; int pi=prime[j];
                int temp=i/g[i];
                if (cnt[t]==2) f[t]=-pi*pi*pi*f[temp];
                else f[t]=0;
                break;
            }
            mu[t]=-mu[i]; cnt[t]=1; g[t]=prime[j]; 
            f[t]=f[i]*f[prime[j]];
        }
    }
    for (int i=1;i<=N;i++) f[i]=f[i]+f[i-1];
}
int calc(int x,int y)
{
    int t1=(x+1)*x>>1; int t2=(y+1)*y>>1;
    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=0;
      for (int i=1,j;i<=n;i=j+1) {
        j=min(n/(n/i),m/(m/i));
        ans=ans+(f[j]-f[i-1])*calc(n/i,m/i);
      }
      printf("%d\n",(ans%p+p)%p);
    }
}