Crash的数字表格（HYSBZ-2154）

Problem Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。

一个4*5的表格如下：

1	2	3	4	5
2	2	6	4	10
3	6	3	12	15
4	4	12	4	20

看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

100%的数据满足N, M ≤ 10^7。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Examples

Input

4 5

Output

122

思路：

题目本质是给出 n、m 两个数，求 

由于 ，因此 

将 gcd(i,j) 提出来，有：

把 d 提出来，有：

设 ，则有：

设 ，那么，

再设 

可以看出，f(d)、g(d) 符合莫比乌斯反演的倍数关系，即：

对 g(d) 进行化简，有：

那么，

因此，需要预处理 u(i)*i^2 的前缀和，然后利用整除分块求出后面的 ij 的和算出 f(1) 后，再计算 res

Source Program

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<string>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<utility>
#include<stack>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<bitset>
#define PI acos(-1.0)
#define INF 0x3f3f3f3f
#define LL long long
#define Pair pair<int,int>
LL quickPow(LL a,LL b){ LL res=1; while(b){if(b&1)res*=a; a*=a; b>>=1;} return res; }
LL quickModPow(LL a,LL b,LL mod){ LL res=1; a=a%mod; while(b){if(b&1)res=(a*res)%mod; a=(a*a)%mod; b>>=1;} return res; }
LL getInv(LL a,LL mod){ return quickModPow(a,mod-2,mod); }
const double EPS = 1E-10;
const int MOD = 20101009;
const int N = 12000000+5;
const int dx[] = {-1,1,0,0,-1,-1,1,1};
const int dy[] = {0,0,-1,1,-1,1,-1,1};
using namespace std;

int mu[N];
int prime[N];
bool bprime[N];
int cnt;
LL sum[N];
void getMu(int n){//线性筛求莫比乌斯函数
    cnt=0;
    mu[1]=1;//根据定义，μ(1)=1
    memset(bprime,false,sizeof(bprime));
 
    for(int i=2;i<=n;i++){//求2~n的莫比乌斯函数
        if(!bprime[i]){
            prime[++cnt]=i;//存储质数
            mu[i]=-1;//i为质数时，μ(1)=-1
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){//枚举i之前的素数个数
            bprime[i*prime[j]]=true;//不是质数
            if(i%prime[j])//i不是prime[j]的整数倍时，i*prime[j]就不会包含相同质因子
                mu[i*prime[j]]=-mu[i];//mu[k]=mu[i]*mu[prime[j]]，因为prime[j]是质数，mu值为-1
            else{
                mu[i*prime[j]]=0;
                break;//留到后面再筛
            }
        }
    }
}
LL cal(LL n,LL m){
    LL res=0;
    LL minn=min(n,m);
    for(LL left=1,right;left<=minn;left=right+1){
        right=min(n/(n/left),m/(m/left));
        LL G=( ((1+n/left)*(n/left)/2)%MOD ) * ( ((1+m/left)*(m/left)/2)%MOD ) %MOD;
        res=( res+((sum[right]-sum[left-1])%MOD*G%MOD)%MOD )%MOD;
    }
    return (res+MOD)%MOD;
}
int main(){
    getMu(12000000);
    LL n,m;
    scanf("%lld%lld",&n,&m);
    
    LL minn=min(n,m);
    for(LL i=1;i<=minn;i++){
        LL temp=(((i%MOD*i%MOD)%MOD)%MOD*mu[i]%MOD)%MOD;
        sum[i]=(sum[i-1]+temp+MOD)%MOD;
    }
    
    LL res=0;
    for(LL left=1,right;left<=minn;left=right+1){
        right=min(n/(n/left),m/(m/left));
        LL temp=(left+right)*(right-left+1)/2%MOD;
        res=(res+(cal(n/left,m/left)%MOD*temp%MOD)%MOD)%MOD;
    }
    printf("%lld\n",res);
    // system("pause");
    return 0; 
}