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

 1 /*bzoj 2154*/
 2 #include <bits/stdc++.h>
 3 
 4 using namespace std;
 5 #define ll long long
 6 #define N 10000000
 7 const int MOD = 20101009;
 8 int mu[N+5] , prime[N+5] , f[N+5] , tot;
 9 bool check[N+5];
10 
11 void init()
12 {
13     mu[1] = 1 , f[1] = 1;
14     for(int i=2 ; i<=N ; i++){
15         if(!check[i]){
16             prime[tot++] = i;
17             mu[i] = -1;
18             f[i] = 1-i;
19         }
20         for(int j=0 ; j<tot ; j++){
21             if((ll)prime[j]*i>N) break;
22             check[prime[j]*i] = true;
23             if(i%prime[j]){
24                 mu[prime[j]*i] = -mu[i];
25                 f[prime[j]*i] = ((ll)f[prime[j]]*f[i])%MOD;
26             }else{
27                 f[prime[j]*i] = f[i];
28                 break;
29             }
30         }
31     }
32 }
33 
34 int q_pow(int b)
35 {
36     ll ans = 1 , a=4;
37     while(b)
38     {
39         if(b&1) ans = (ans*a)%MOD;
40         a = (a*a)%MOD;
41         b>>=1;
42     }
43     return (int)ans;
44 }
45 
46 int solve(int n , int m)
47 {
48     if(n>m) swap(n , m);
49     ll ret = 0;
50     for(int k=1 ; k<=n ; k++){
51         ll a = n/k , b = m/k;
52         ret = (ret+(k*(1+a)%MOD)*a%MOD*b%MOD*(1+b)%MOD*f[k]%MOD)%MOD;
53     }
54     ret = ((ll)ret*q_pow(MOD-2)%MOD+MOD)%MOD;
55     return (int)ret;
56 }
57 
58 int main()
59 {
60   //  freopen("in.txt" , "r" , stdin);
61     init();
62     int n , m;
63    /* for(int i=2 ; i<=100000 ; i++)
64         if(!check[i] && MOD%i==0) cout<<i<<" yes "<<endl;*/
65     while(~scanf("%d%d" , &n , &m)){
66         printf("%d\n" , solve(n , m));
67     }
68     return 0;
69 }