[BZOJ2154]Crash的数字表格（莫比乌斯反演）

题目描述

传送门

题解

是时候又写一大坨推导公式了。
首先设

$$s(x,y)=\sum\limits_{i=1}^x\sum\limits_{j=1}^y ij[gcd(i,j)=1]$$
并且 $x<y$
利用莫比乌斯反演变形 $$s(x,y)=\sum\limits_{i=1}^x\sum\limits_{j=1}^y ij\sum\limits_{d=1}^x[d|i][d|j]\mu(d)$$
令i=id j=jd $$s(x,y)=\sum\limits_{d=1}^x\mu(d)d^2\sum\limits_{i=1}^{x\over d}\sum\limits_{j=1}^{y\over d} ij$$

再设

$$sum(x,y)=\sum\limits_{i=1}^x\sum\limits_{j=1}^y ij$$
利用等差数列求和公式
$$sum(x,y)={{x(x+1)}\over 2}*{{y(y+1)}\over 2}$$
将sum函数代入s函数得 $$s(x,y)=\sum\limits_{d=1}^x\mu(d)d^2 sum(\lfloor{x\over d}\rfloor,\lfloor{y\over d}\rfloor)$$

好我们现在回到题目

$$\sum\limits_{i=1}^n\sum\limits_{j=1}^mlcm[i,j]$$
设 $n<m$
那么
$$\sum\limits_{i=1}^n\sum\limits_{j=1}^m{{ij}\over gcd(i,j)}$$
$$\sum\limits_{i=1}^n\sum\limits_{j=1}^mij\sum\limits_{d=1}^n[d|i][d|j][gcd(i,j)=d]{1\over d}$$
令i=id j=jd
$$\sum\limits_{d=1}^n\sum\limits_{i=1}^{n\over d}id\sum\limits_{j=1}^{m\over d}jd[gcd(i,j)=1]{1\over d}$$
$$\sum\limits_{d=1}^nd\sum\limits_{i=1}^{n\over d}\sum\limits_{j=1}^{m\over d}ij[gcd(i,j)=1]$$

将s函数代入这个式子得到

$$\sum\limits_{d=1}^nd s(\lfloor{n\over d}\rfloor,\lfloor{m\over d}\rfloor)$$

sum函数可以在$O(1)$的时间得到
s函数最多有$\sqrt n+\sqrt m$种取值，可以在$O(\sqrt n+\sqrt m)$的时间得到
那么答案也只有$\sqrt n+\sqrt m$种取值，同样可以在$O(\sqrt n+\sqrt m)$的时间得到，所以时间复杂度是接近$O(n)$的

代码

#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
#define LL long long
const int max_n=1e7+5;
const LL Mod=20101009;
LL n,m,x,y,xx,yy,k,jj,sum,ans,s[max_n];
int mu[max_n],prime[max_n],p[max_n];
inline void get_mu(){
    mu[1]=1;
    for (int i=2;i<=n;++i){
        if (!p[i]){
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=prime[0]&&i*prime[j]<=n;++j){
            p[i*prime[j]]=1;
            if (i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for (LL i=1;i<=n;i++) s[i]=(LL)(s[i-1]+(LL)(i*i*mu[i])%Mod)%Mod;  
}
inline LL Sum(LL x,LL y){
    LL w1=((x+1)*x/2)%Mod;
    LL w2=((y+1)*y/2)%Mod;
    return w1*w2%Mod;
}
int main(){
    scanf("%lld%lld",&n,&m); if (n>m) swap(n,m);
    get_mu(); jj=0;
    for (LL i=1;i<=n;i=jj+1){
        jj=min(n/(n/i),m/(m/i));
        x=n/i,y=m/i; if (x>y) swap(x,y);
        k=0; sum=0;
        for (LL j=1;j<=x;j=k+1){
            k=min(x/(x/j),y/(y/j));
            xx=x/j; yy=y/j;
            sum=(sum+((s[k]-s[j-1])%Mod*Sum(xx,yy)%Mod)%Mod)%Mod;
        }
        ans=(ans+(((i+jj)*(jj-i+1)/2)%Mod*sum%Mod)%Mod)%Mod;
    }
    printf("%lld\n",(ans+Mod)%Mod);
}