bzoj2154(莫比乌斯反演)

//
//  main.cpp
//  bzoj2154
//
//  Created by New_Life on 16/7/7.
//  Copyright © 2016年 chenhuan001. All rights reserved.
//

#include <iostream>
#include <string.h>
#include <stdio.h>
using namespace std;

#define N 10001000
#define MOD 20101009

//－－莫比乌斯反演函数－－//
//说明：利用线性素数筛选顺便求了个mu
//注释部分为求从区间[1,b]和区间[1,d]中取两个数，互质对数O（n^0.5）
//复杂度:O(n)
int mu[N];
long long sum[N];
int prime[N];
bool mark[N];

void mobus()
{
    int pcnt=0;
    memset(mark,0,sizeof(mark));
    mu[1] = 1;
    for(int i=2;i<N;i++)
    {
        if(mark[i] == 0)
        {
            prime[pcnt++] = i;
            mu[i] = -1;
        }
        for(int j=0;j<pcnt && i*prime[j]<N;j++)
        {
            int tmp = i*prime[j];
            mark[tmp] = 1;
            if( i%prime[j] == 0 )
            {
                mu[tmp] = 0;
                break;
            }
            
            mu[tmp] = mu[i]*-1;
        }
    }
    for(int i=1;i<N;i++)
    {
        sum[i] += sum[i-1]+(long long)mu[i]*i*i;
        sum[i] %= MOD;
    }
}



long long gaobili(long long b,long long d)
{
    if(b<=0||d<=0) return 0;
    long long m = min(b,d);
    long long ans = 0;
    while(m>=1)
    {
        long long tb = b/( b/m +1 )+1;
        long long td = d/( d/m +1 )+1;
        //前进的最大位置
        long long tm = max(tb,td);
        ans += (sum[m]-sum[tm-1])*(((b/m+1)*(b/m)/2)%MOD)%MOD*(((d/m+1)*(d/m)/2)%MOD)%MOD ;
        ans %= MOD;
        m = tm-1;
    }
    return ans;
}
//等差数列求和模板，[a1,a1+d,...,an]
long long allsum(long long a1,long long an,long long n)
{
    if(n%2==0)
        return (a1+an)*(n/2);
    else return ((a1+an)/2)*n;
}

int main(int argc, const char * argv[]) {
    mobus();
    int b,d;
    while(scanf("%d%d",&b,&d)!=EOF)
    {
        int m = min(b,d);
        long long ans = 0;
        while(m>=1)
        {
            int tb = b/( b/m +1 )+1;
            int td = d/( d/m +1 )+1;
            //前进的最大位置
            int tm = max(tb,td);
            ans += allsum(tm,m,m-tm+1)%MOD*gaobili(b/m, d/m)%MOD;
            ans %= MOD;
            m = tm-1;
        }
        cout<<(ans+MOD)%MOD<<endl;
    }
    return 0;
}
/*
 4 5
 
 */