BZOJ 2154 Crash 的数字表格

显然题目就是要我们求 n∑i=1m∑j=1lcm(i,j)。
根据小学奥数，就是 n∑i=1m∑j=1ijgcd(i,j)。

n∑i=1m∑j=1ijgcd(i,j)=min(n,m)∑i=1i∑xi≤ni∑yi≤mi[gcd(x,y)=1]xy

敏锐地察觉到 [n=1] 这个东西可以替换成 ϵ(n)=∑d|nμ(d)。
那么有

min(n,m)∑i=1i∑xi≤ni∑yi≤mi[gcd(x,y)=1]xy=min(n,m)∑i=1i∑xi≤ni∑yi≤mixy∑d|gcd(x,y)μ(d)

根据套路，接着把 d 提出来。

min(n,m)∑i=1i∑xi≤ni∑yi≤mixy∑d|gcd(x,y)μ(d)=min(n,m)∑i=1imin(⌊ni⌋,⌊mi⌋)∑d=1μ(d)d2∑xd≤ni∑yd≤mixy=min(n,m)∑i=1imin(⌊ni⌋,⌊mi⌋)∑d=1μ(d)d2sum(⌊nid⌋,⌊mid⌋)

其中 sum(x,y)=x∑i=1y∑j=1ij=xy(x+1)(y+1)4。
那么预处理 μ(i)i2 的前缀和，做两次数论分块即可。

代码：

#include <cstdio>
#include <algorithm>
using namespace std;
const int N = 1e7;
const long long mod = 20101009;
int n,m;
int cnt,prime[N + 5],vis[N + 5],sum[N + 5];
long long ans;
long long get_sum(int n,int m)
{
    long long ret = 0;
    if(n > m)
        swap(n,m);
    for(register int l = 1,r;l <= n;l = r + 1)
    {
        r = min(n / (n / l),m / (m / l));
        ret = (ret + (long long)(sum[r] - sum[l - 1] + mod) * ((long long)(n / l) * (n / l + 1) / 2 % mod) % mod * ((long long)(m / l) * (m / l + 1) / 2 % mod) % mod);
    }
    return ret;
}
int main()
{
    sum[1] = 1;
    for(register int i = 2;i <= N;++i)
    {
        if(!vis[i])
            prime[++cnt] = i,sum[i] = -1;
        for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
                break;
            sum[i * prime[j]] = -sum[i];
        }
    }
    for(register int i = 1;i <= N;++i)
        sum[i] = (sum[i - 1] + (long long)i * i % mod * (sum[i] + mod)) % mod;
    scanf("%d%d",&n,&m);
    if(n > m)
        swap(n,m);
    for(register int l = 1,r;l <= n;l = r + 1)
    {
        r = min(n / (n / l),m / (m / l));
        ans = (ans + (long long)(l + r) * (r - l + 1) / 2 % mod * get_sum(n / l,m / l) % mod) % mod;
    }
    printf("%lld\n",ans);
}