Description
今天的數學課上，Crash小朋友學習了最小公倍數(Least Common Multiple)。對於兩個正整數a和b，LCM(a, b)表示能同時被a和b整除的最小正整數。例如，LCM(6, 8) = 24。回到家後，Crash還在想著課上學的東西，為了研究最小公倍數，他畫了一張NM的表格。每個格子裡寫了一個數字，其中第i行第j列的那個格子裡寫著數為LCM(i, j)。一個45的表格如下： 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。
Output
輸出一個正整數，表示表格中所有數的和mod 20101009的值。
Sample Input
4 5
Sample Output
122
【資料規模和約定】
100%的資料滿足N, M ≤ 10^7。

是bzoj2693的弱化版
不用最後一步推導就能A了
詳細見那道題題解：https://blog.csdn.net/a1035719430/article/details/84961774

#include<bits/stdc++.h>
using namespace std;
#define rep(i,j,k) for(int i = j;i <= k;++i)
#define repp(i,j,k) for(int i = j;i >= k;--i)
#define rept(i,x) for(int i = linkk[x];i;i = e[i].n)
#define
 
 P pair<int,int>
#define Pil pair<int,ll>
#define Pli pair<ll,int>
#define Pll pair<ll,ll>
#define pb push_back 
#define pc putchar
#define mp make_pair 
#define file(k) memset(k,0,sizeof(k))
#define ll long long
#define N 10000000
int rd()
{
	int num = 0;char c = getchar();bool flag =
 
 true;
	while(c < '0'||c > '9') {if(c == '-') flag = false;c = getchar();}
	while(c >= '0' && c <= '9') num = num*10+c-48,c = getchar();
	if(flag) return num;else return -num;
}
const int p = 20101009;
int n,m,inv;
int f[N+10000],prime[N+10000],miu[N+10000];
inline int mo(int a){return a>=p?a-p:a;}
inline int calc(int a,int b){return (a+=b)>=p?a-=p:a;}
inline int mul(int a,int b){return 1ll*a*b%p;}
inline int del(int a,int b){return (a-=b)<0?a+=p:a;}
void pre()
{
	miu[1] = 1;f[1] = miu[1]*1;
	rep(i,2,N)
	{
		if(!f[i]) miu[i] = p-1,f[i] = calc(f[1],mul(miu[i],i)),prime[++prime[0]] = i;
		rep(j,1,prime[0])
		{
			if(i*prime[j]>N) break;
			miu[i*prime[j]] = i%prime[j]?mo(p-miu[i]):0;
			f[i*prime[j]] = i%prime[j]?mul(f[i],f[prime[j]]):f[i];
			if(i%prime[j] == 0) break;
		}
	}
	rep(i,1,N) f[i] = mul(f[i],i);
	rep(i,1,N) f[i] = calc(f[i],f[i-1]);
}
void work()
{
	n = rd();m = rd();
	if(n < m) swap(n,m);
	int i = 1;
	int ans = 0;
	while(i<=m)
	{
		int j = min(n/(n/i),m/(m/i)),tmpn = n/i,tmpm = m/i;
		int v = mul(inv,mul(1+tmpn,tmpn));v = mul(v,mul(inv,mul(1+tmpm,tmpm)));
		ans = calc(ans,mul(v,del(f[j],f[i-1])));
		i = j+1;
	}
	printf("%d\n",ans);
}
int main()
{
	inv = (p+1)/2;
	pre();
	work();
	return 0;
}