洛谷-1829 [国家集训队]Crash的数字表格 / JZPTAB

题目背景
提示：原 P1829 半数集问题 已经迁移至 P1028 数的计算
题目描述
今天的数学课上，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只想知道表格里所有数的和mod20101009的值。
输入格式
输入的第一行包含两个正整数，分别表示N和M。
输出格式
输出一个正整数，表示表格中所有数的和mod20101009的值。

输入输出样例
输入 #1
4 5

输出 #1
122

说明/提示
100%的数据满足N, M≤ 10^7。

解释:推公式呗~转大佬的来学习一下
Solution
易知原式等价于
$\sum _{i=1}^n\sum _{j=1}^m\frac{i\cdot j}{\gcd (i,j)}$
枚举最大公因数 d，显然两个数除以 d 得到的数互质
$\sum _{i=1}^n\sum _{j=1}^m\sum _{d|i,d|j,\gcd (\frac{i}{d},\frac{j}{d})=1}\frac{i\cdot j}{d}$
非常经典的 $\gcd$ 式子的化法
$\sum _{d=1}^{n}d\cdot \sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[\gcd (i,j)=1]\cdot i\cdot j$
后半段式子中，出现了互质数对之积的和，为了让式子更简洁就把它拿出来单独计算。于是我们记
$\mathrm{sum}(n,m)=\sum _{i=1}^n\sum _{j=1}^m[\gcd (i,j)=1]\cdot i\cdot j$
接下来对 $\mathrm{sum}(n,m)$ 进行化简。首先枚举约数，并将 $[\gcd (i,j)=1]$ 表示为 $\epsilon (\gcd (i,j))\epsilon (gcd(i,j))\sum _{d=1}^n\sum _{d|i}^n\sum _{d|j}^m\mu (d)\cdot i\cdot j$
指上式中的 i,ji,j ），显然式子可以变为
$\sum _{d=1}^n\mu (d)\cdot d^2\cdot \sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }i\cdot j$
观察上式，前半段可以预处理前缀和；后半段又是一个范围内数对之和，记
$g(n,m)=\sum _{i=1}^n\sum _{j=1}^{m}i\cdot j=\frac{n\cdot (n+1)}{2}\times \frac{m\cdot (m+1)}{2}$
至此
$\mathrm{sum}(n,m)=\sum _{d=1}^n\mu (d)\cdot d^2\cdot g(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )$
我们可以 $\lfloor \frac{n}{\lfloor \frac{n}{d}\rfloor }\rfloor$
数论分块求解 $\mathrm{sum}(n,m)$函数。
在求出 $\mathrm{sum}(n,m)$ 后，回到定义 $\mathrm{sum}$ 的地方，可得原式为
$\sum _{d=1}^{n}d\cdot \mathrm{sum}(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )$
可见这又是一个可以数论分块求解的式子！
最后分块就好了

#include <cstdio>
#include <algorithm>
using namespace std; 
const int N=1e7;
const int mod=20101009;
int n,m,mu[N+5],p[N/10+5],sum[N+5];
bool flg[N+5];
void init() {
    mu[1]=1;
    int tot=0,k=min(n,m);
    for(int i=2;i<=k;++i) {
        if(!flg[i]) p[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot&&i*p[j]<=k;++j) {
            flg[i*p[j]]=1;
            if(i%p[j]==0) {mu[i*p[j]]=0;break;}
            mu[i*p[j]]=-mu[i];
        }
    }
    for(int i=1;i<=k;++i) sum[i]=(sum[i-1]+1LL*i*i%mod*(mu[i]+mod))%mod;
}
int Sum(int x,int y) {
    return (1LL*x*(x+1)/2%mod)*(1LL*y*(y+1)/2%mod)%mod;
}
int func(int x,int y) {
    int res=0;
    for(int i=1,j;i<=min(x,y);i=j+1) {
        j=min(x/(x/i),y/(y/i));
        res=(res+1LL*(sum[j]-sum[i-1]+mod)*Sum(x/i,y/i)%mod)%mod;
    }
    return res;
}
int solve(int x,int y) {
    int res=0;
    for(int i=1,j;i<=min(x,y);i=j+1) {
        j=min(x/(x/i),y/(y/i));
        res=(res+1LL*(j-i+1)*(i+j)/2%mod*func(x/i,y/i)%mod)%mod;
    }
    return res;
}
int main() {
    scanf("%d%d",&n,&m);
    init();
    printf("%d\n",solve(n,m));
}