bzoj 2694: Lcm （反演）

题目描述

传送门

题目大意：
定义整数a,b,求所有满足条件的lcm(a,b)的和：
$1<=a<=A$
$1<=b<=B$
$\forall n>1, n^2 \dagger gcd(a,b)$
输出答案对2^30取模。

题解

要求gcd(a,b)不能含平方因子，所以gcd(a,b)一定是mu不等于0的数。
那么我们设所有满足条件的数为p，n<=m
就可以将上述条件转化成式子

$$\sum\limits_{p}\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(a,b)[gcd(a,b)=p]$$
$$\sum\limits_{p}\sum\limits_{i=1}^n\sum\limits_{j=1}^m{i*j\over p}[gcd(a,b)=p]$$
$$\sum\limits_{p} p \sum\limits_{d=1}^n \mu(d) *d^2 sum({\lfloor {n\over d*p}\rfloor},{\lfloor {m\over d*p}\rfloor})$$
其中 $sum(x,y)={x*(x+1)\over 2}*{y*(y+1)\over 2}$
设T=d*p
$$\sum\limits_{T=1}^n sum({\lfloor {n\over T}\rfloor},{\lfloor {m\over T}\rfloor})\sum\limits_{p|T} \mu({T\over p})*({T\over p})^2*p$$
这个式子的话，要求 $h(T)=\sum\limits_{p|T} \mu({T\over p})*({T\over p})^2*p$
我们其实可以直接枚举mu不等于0的数，然后用它更新它的倍数，这样的时间复杂度应该均摊是不超过 $O(n\log n)$
还有这道题的模数非常特殊，是2的整数次幂，所以我们直接用int自然溢出，最后取模即可。

代码

#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];
int f[N+3];
void init()
{
    mu[1]=1;
    for (int i=2;i<=N;i++) {
        if (!pd[i]) {
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        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;
                break;
            }
            mu[t]=-mu[i];
        }
    }
    for (int i=1;i<=N;i++) 
      if (mu[i]){
        int t=i;
        for (int j=1;j*t<=N;j++) f[j*t]=f[j*t]+mu[j]*j*j*i;
    }
    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);
    }
}