【莫比乌斯反演】[BZOJ2154]Crash的数字表格

我膜拜PoPoQQQ

http://www.lydsy.com/JudgeOnline/problem.php?id=2154

∑ni=1∑mj=1lcm(i,j)=∑ni=1∑mj=1i×jgcd(i,j)\sum_{i=1}^n\sum_{j=1}^mlcm(i, j)=\sum_{i=1}^n\sum_{j=1}^m\frac{i\times j}{gcd(i, j)}

令gcd(i,j)=dgcd(i, j)=d那么∑ni=1∑mj=1i×jgcd(i,j)=∑min(n,m)d=1d2×g(n,m,d)d\sum_{i=1}^n\sum_{j=1}^m\frac{i\times j}{gcd(i, j)}=\sum_{d=1}^{min(n, m)}\frac{d^2\times g(n,m,d)}{d}

这里令g(n,m,d)g(n, m, d)表示为在

g(n,m,d)=∑⌊nd⌋i=1∑⌊md⌋j=1[gcd(i,j)==1]i×j=∑⌊nd⌋i=1∑⌊md⌋j=1∑t|i,t|jμ(t)i×jg(n,m,d)=\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}[gcd(i,j)==1]i\times j=\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}\sum_{t|i,t|j} \mu(t) i\times j

那么就有∑min(n,m)d=1d2×g(n,m,d)d=∑min(n,m)d=1d×g(n,m,d)\sum_{d=1}^{min(n, m)}\frac{d^2\times g(n,m,d)}{d}=\sum_{d=1}^{min(n, m)}d\times g(n,m,d)

那么令Sum(a,b)=∑ai=1∑bj=1i×j=(1+2....+a)×(1+2.....+b)=a(a+1)2×b(b+1)2Sum(a,b)=\sum_{i=1}^{a}\sum_{j=1}^{b}i\times j=(1+2....+a)\times(1+2.....+b)=\frac{a(a+1)}{2}\times\frac{b(b+1)}{2}

因为g(n,m,d)=∑⌊nd⌋i=1∑⌊md⌋j=1∑t|i,t|jμ(t)i×j=∑min(n,m)dt=1μ(t)×t2×Sum(⌊ntd⌋,⌊mtd⌋)g(n, m, d)=\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}\sum_{t|i,t|j} \mu(t) i\times j=\sum_{t=1}^{\frac{min(n, m)}{d}}\mu(t)\times t^2\times Sum(\lfloor \frac{n}{td}\rfloor, \lfloor \frac{m}{td}\rfloor)

[code]#include <cstdio>
#include <iostream>
#include <cstring>
#include <climits>
#include <algorithm>
using namespace std;
#define mod 20101009LL
const int MAXN = 10000000;
bool notprime[MAXN+10];
int prime[MAXN+10], mu[MAXN+10], sum[MAXN+10];
long long Max;
void Init(){
    mu[1] = 1;
    int tmp;
    for(int i=2;i<=Max;i++){
        if(!notprime[i]){
            mu[i] = -1;
            prime[++prime[0]] = i;
        }
        for(int j=1;j<=prime[0]&&(tmp=prime[j]*i)<=Max;j++){
            notprime[tmp] = true;
            if(i%prime[j] == 0){
                mu[tmp] = 0;
                break;
            }
            mu[tmp] = -mu[i];
        }
    }
    for(long long i=1;i<=Max;i++)
        sum[i]=(sum[i-1]+(i*i*mu[i])%mod)%mod;
}
long long Sum(long long n, long long m){
    return (n*(n+1)/2%mod)*(m*(m+1)/2%mod)%mod;
}
long long F(long long n, long long m){
    long long ret = 0, last;
    for(long long i=1;i<=n;i=last+1){
        last = min(n/(n/i), m/(m/i));
        ret += (sum[last] - sum[i-1]) * Sum(n/i, m/i)%mod;
        ret %= mod;
    }
    return ret;
}
long long solve(int n, int m){
    long long ret = 0, last;
    for(long long i=1;i<=n;i=last+1){
        last = min(n/(n/i),  m/(m/i));
        ret += (i+last)*(last-i+1)/2 * F(int(n/i), int(m/i))%mod;
        ret %= mod;
    }
    return ret;
}
int main(){
    long long a, b;
    cin>>a>>b;
    if(a > b) swap(a, b);
    Max = a;
    Init();
    cout<<(solve(a, b)%mod+mod)%mod<<endl;

    return 0;
}