【莫比乌斯反演】[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 20101009LLconst 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;}