bzoj2154 Crash的数字表格

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。 Output

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

$$\begin{align} & \sum_{i=1}^n\sum_{j=1}^m lcm(i,j) \\ = & \sum_{i=1}^n\sum_{j=1}^m \frac{ij} {\gcd(i,j)} \\ = & \sum_{d=1}^{\min(n,m)} \sum_{i=1}^{\lfloor\frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor} \frac{id*jd}d e(\gcd(i,j)) \\ = & \sum_{d=1}^{\min(n,m)}d \sum_{i=1}^{\lfloor\frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor} ij\sum_{x|i,x|j}\mu(x) \\ = & \sum_{d=1}^{\min(n,m)}d \sum_{x=1}^{\min(\lfloor\frac nd\rfloor,\lfloor\frac md \rfloor)} \mu(x)\sum_{i=1}^{\lfloor\frac n{dx} \rfloor}\sum_{j=1}^{\lfloor \frac m{dx} \rfloor} ix*jx \\ = & \sum_{d=1}^{\min(n,m)}d \sum_{x=1}^{\min(\lfloor\frac nd\rfloor,\lfloor\frac md \rfloor)} \mu(x)*x^2\frac{\lfloor \frac n{dx}\rfloor(1-\lfloor \frac n{dx}\rfloor)\lfloor \frac m{dx}\rfloor(1-\lfloor \frac m{dx}\rfloor)}4 \\ \end{align}$$

把

$$\sum_{x=1}^{\min(\lfloor\frac nd\rfloor,\lfloor\frac md \rfloor)} \mu(x)*x^2\frac{\lfloor \frac n{dx}\rfloor(1-\lfloor \frac n{dx}\rfloor)\lfloor \frac m{dx}\rfloor(1-\lfloor \frac m{dx}\rfloor)}4$$
看成关于 $\lfloor\frac nd\rfloor$和 $\lfloor\frac md\rfloor$的函数，对 $d$进行第一次下底函数分块，需要预处理$id$前缀和。求上面的函数的时候对 $x$进行第二次下底函数分块，需要预处理$\mu*id^2$的前缀和。
复杂度 $O(\sqrt n*\sqrt n)=O(n)$。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LL long long
const int p=20101009,inv=15075757;
int prm[10000010],have[10000010],mu[10000010],sf[10000010],sid[10000010],tot;
int calc(int n,int m)
{
    int ret=0,x;
    for (int i=1;i<=n;i=x+1)
    {
        x=min(n/(n/i),m/(m/i));
        ret=(ret+(LL)(sf[x]-sf[i-1])*(n/i)%p*(1+n/i)%p*(m/i)%p*(1+m/i)%p*inv%p+p)%p;
    }
    return ret;
}
int main()
{
    int m,n,x,ans=0;
    scanf("%d%d",&n,&m);
    if (n>m) swap(n,m);
    sf[1]=sid[1]=1;
    for (int i=2;i<=m;i++)
    {
        if (!have[i])
        {
            prm[++tot]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=tot&&(LL)i*prm[j]<=m;j++)
        {
            have[i*prm[j]]=1;
            if (i%prm[j]) mu[i*prm[j]]=-mu[i];
            else
            {
                mu[i*prm[j]]=0;
                break;
            }
        }
        sf[i]=(sf[i-1]+(LL)i*i*mu[i]%p+p)%p;
        sid[i]=(sid[i-1]+i)%p;
    }
    for (int i=1;i<=n;i=x+1)
    {
        x=min(n/(n/i),m/(m/i));
        ans=(ans+(LL)(sid[x]-sid[i-1])*calc(n/i,m/i)%p+p)%p;
    }
    printf("%d\n",ans);
}