bzoj2154 Crash的数字表格（莫比乌斯反演）

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;
}