【莫比乌斯反演】[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)
令gcd(i,j)=d那么∑ni=1∑mj=1i×jgcd(i,j)=∑min(n,m)d=1d2×g(n,m,d)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×j
那么就有∑min(n,m)d=1d2×g(n,m,d)d=∑min(n,m)d=1d×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)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⌋)

#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;
}