bzoj 2154 莫比乌斯反演求lcm的和

 题目大意：

表格中每一个位置(i,j)填的值是lcm(i,j) , 求n*m的表格值有多大

论文贾志鹏线性筛中过程讲的很好

最后的逆元我利用的是欧拉定理求解的

我这个最后线性扫了一遍，勉强过了，效率不是很高。。。

/*bzoj 2154*/
#include <bits/stdc++.h>

using namespace std;
#define ll long long
#define N 10000000
const int MOD = 20101009;
int mu[N+5] , prime[N+5] , f[N+5] , tot;
bool check[N+5];

void init()
{
    mu[1] = 1 , f[1] = 1;
    for(int i=2 ; i<=N ; i++){
        if(!check[i]){
            prime[tot++] = i;
            mu[i] = -1;
            f[i] = 1-i;
        }
        for(int j=0 ; j<tot ; j++){
            if((ll)prime[j]*i>N) break;
            check[prime[j]*i] = true;
            if(i%prime[j]){
                mu[prime[j]*i] = -mu[i];
                f[prime[j]*i] = ((ll)f[prime[j]]*f[i])%MOD;
            }else{
                f[prime[j]*i] = f[i];
                break;
            }
        }
    }
}

int q_pow(int b)
{
    ll ans = 1 , a=4;
    while(b)
    {
        if(b&1) ans = (ans*a)%MOD;
        a = (a*a)%MOD;
        b>>=1;
    }
    return (int)ans;
}

int solve(int n , int m)
{
    if(n>m) swap(n , m);
    ll ret = 0;
    for(int k=1 ; k<=n ; k++){
        ll a = n/k , b = m/k;
        ret = (ret+(k*(1+a)%MOD)*a%MOD*b%MOD*(1+b)%MOD*f[k]%MOD)%MOD;
    }
    ret = ((ll)ret*q_pow(MOD-2)%MOD+MOD)%MOD;
    return (int)ret;
}

int main()
{
  //  freopen("in.txt" , "r" , stdin);
    init();
    int n , m;
   /* for(int i=2 ; i<=100000 ; i++)
        if(!check[i] && MOD%i==0) cout<<i<<" yes "<<endl;*/
    while(~scanf("%d%d" , &n , &m)){
        printf("%d
" , solve(n , m));
    }
    return 0;
}