BZOJ2693 jzptab 莫比乌斯反演

Posted 2020-09-01 WDZRMPCBIT

题意：给定N,M，求$\\sum\\limits_{i = 1}^N {\\sum\\limits_{j = 1}^M {lcm(i,j)} }$，多组询问

题解：

设$F(x,y) = \\sum\\limits_{i = 1}^x {\\sum\\limits_{j = 1,\\gcd (i,j) = 1}^y {ij} }$  $S(x,y) = \\frac{{x(x + 1)}}{2}\\frac{{y(y + 1)}}{2}$

我们枚举最大公因数d，则有\\[\\sum\\limits_{i = 1}^N {\\sum\\limits_{j = 1}^M {lcm(i,j)} }  = \\sum\\limits_{i = 1}^N {\\sum\\limits_{j = 1}^M {\\frac{{ij}}{{\\gcd (i,j)}}} }  = \\frac{{{d^2}F(\\left\\lfloor {\\frac{N}{d}} \\right\\rfloor ,\\left\\lfloor {\\frac{M}{d}} \\right\\rfloor )}}{d} = dF(\\left\\lfloor {\\frac{N}{d}} \\right\\rfloor ,\\left\\lfloor {\\frac{M}{d}} \\right\\rfloor )\\]

反演得到\\[F(x,y) = \\sum\\limits_{i = 1}^{\\min (x,y)} {{i^2}\\mu (i)S(\\left\\lfloor {\\frac{x}{i}} \\right\\rfloor ,\\left\\lfloor {\\frac{y}{i}} \\right\\rfloor )} \\]

继续优化\\[Ans = \\sum\\limits_{d = 1}^{\\min (N,M)} {[d\\sum\\limits_{i = 1}^{\\min (N,M)} {{i^2}\\mu (i)S(\\left\\lfloor {\\frac{x}{{di}}} \\right\\rfloor ,\\left\\lfloor {\\frac{y}{{di}}} \\right\\rfloor )} ]}  = \\sum\\limits_{D = 1}^{\\min (N,M)} S (\\left\\lfloor {\\frac{N}{D}} \\right\\rfloor ,\\left\\lfloor {\\frac{M}{D}} \\right\\rfloor )\\sum\\limits_{i|D} {\\frac{D}{i}{i^2}\\mu (i)} \\]

后面那个和式可以通过线性筛预处理出来，后面分块就好

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
using namespace std;
#define ll long long

const int P=20101009;
const int MAXN=10000000+2;
int T;
ll N,M,mu[MAXN],prime[MAXN],f[MAXN];
bool flag[MAXN]; 

void Pre(ll N){
    f[1]=1;
    for(int i=2,j=0;i<=N;i++){
        if(!flag[i]) prime[++j]=i,f[i]=(i-(ll)i*i)%P;
        for(ll k=1;k<=j && i*prime[k]<=N;k++){
            flag[i*prime[k]]=1;
            if(i%prime[k]==0){
                f[i*prime[k]]=(f[i]*prime[k])%P;
                break;
            }
            f[i*prime[k]]=(f[prime[k]]*f[i])%P;
        }
    }
    for(int i=1;i<=N;i++) f[i]=(f[i]+f[i-1])%P;
}

ll Calc(ll x,ll y){
    ll ret1,ret2;
    ret1=((x*(x+1))>>1)%P;
    ret2=((y*(y+1))>>1)%P;
    return (ret1*ret2)%P;
}

int main(){
    Pre(MAXN);

    cin >> T;
    while(T--){
        ll ans=0;
        scanf("%lld %lld",&N,&M);
        if(N<M) swap(N,M);

        for(ll i=1,j;i<=M;i=j+1){
            j=min(N/(N/i),M/(M/i));
            ans=(ans+((f[j]-f[i-1])*Calc(N/j,M/j))%P)%P;
        }

        printf("%lld\\n",(ans+P)%P);
    }

    return 0;
}

View Code