[国家集训队]Crash的数字表格 / JZPTAB

可以知道,题目中所求的是

∑i=1n∑j=1mlcm(i,j)=∑i=1n∑j=1mijgcd(i,j)\sum\limits_{i=1}^n\sum\limits_{j=1}^mlcm(i,j)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{ij}{gcd(i,j)}

枚举所有约数

∑i=1n∑j=1mijgcd(i,j)=∑d=1min(n,m)∑i=1n∑j=1mijd[gcd(i,j)=d]=∑d=1min(n,m)d∑i=1⌊nd⌋∑j=1⌊md⌋ij[gcd(i,j)=1]\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{ij}{gcd(i,j)}=\sum\limits_{d=1}^{min(n,m)}\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{ij}{d}[gcd(i,j)=d]\\=\sum\limits_{d=1}^{min(n,m)}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij[gcd(i,j)=1]

令 sum(n,m)=∑i=1n∑j=1mij[gcd(i,j)=1]f(k)=∑i=1n∑j=1mij[gcd(i,j)=k],F(k)=∑i=1n∑j=1mij[k∣gcd(i,j)]sum(n,m)=\sum\limits_{i=1}^n\sum\limits_{j=1}^mij[gcd(i,j)=1]\\f(k)=\sum\limits_{i=1}^n\sum\limits_{j=1}^mij[gcd(i,j)=k],F(k)=\sum\limits_{i=1}^n\sum\limits_{j=1}^mij[k|gcd(i,j)]

⇒F(k)=∑k∣df(d)⇒f(k)=∑k∣dμ(⌊dk⌋)F(d)=∑k∣dμ(⌊dk⌋)∑i=1n∑j=1mij[d∣gcd(i,j)]=∑k∣dμ(⌊dk⌋)d2∑i=1⌊nd⌋∑j=1⌊md⌋ij\Rightarrow F(k)=\sum\limits_{k|d}f(d)\Rightarrow f(k)=\sum\limits_{k|d}\mu(\lfloor\frac{d}{k}\rfloor)F(d)\\=\sum\limits_{k|d}\mu(\lfloor\frac{d}{k}\rfloor)\sum\limits_{i=1}^n\sum\limits_{j=1}^mij[d|gcd(i,j)]=\sum\limits_{k|d}\mu(\lfloor\frac{d}{k}\rfloor)d^2\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij

令 g(i,j)=i×(i+1)2×j×(j+1)2g(i,j)=\frac{i\times(i+1)}{2}\times\frac{j\times(j+1)}{2}

⇒f(k)=∑k∣dμ(⌊dk⌋)d2g(⌊nd⌋,⌊md⌋)\Rightarrow f(k)=\sum\limits_{k|d}\mu(\lfloor\frac{d}{k}\rfloor)d^2g(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)

⇒sum(n,m)=f(1)=∑d=1min(n,m)μ(d)d2g(⌊nd⌋,⌊md⌋)\Rightarrow sum(n,m)=f(1)=\sum\limits_{d=1}^{min(n,m)}\mu(d)d^2g(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)

ans=∑d=1min(n,m)dsum(⌊nd⌋,⌊md⌋)ans=\sum\limits_{d=1}^{min(n,m)}dsum(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)

推到这就很明显可以用整除分块了,时间复杂度 O(n)O(n),注意取模和long long

代码

#include<iostream>
#include<cstdio>
#define MOD (20101009)
#define g(x, y) ((((ll)(x + 1) * (x) / 2 % MOD) * ((ll)(y + 1) * (y) / 2 % MOD)) % MOD)
using namespace std;
const int N = 1e7 + 5;
typedef long long ll;

bool st[N];
int prime[N], mu[N];
ll s[N];

void init()
{
int tot = 0;
mu[1] = 1;
for (int i = 2; i < N; i++)
{
if (!st[i])
prime[++tot] = i, mu[i] = -1;
for (int j = 1; i * prime[j] < N; j++)
{
st[i * prime[j]] = true;
if (i % prime[j] == 0)
break;
mu[i * prime[j]] = -mu[i];
}
}
for (int i = 1; i < N; i++)
s[i] = ((ll)s[i - 1] + (ll)mu[i] * i * i + MOD) % MOD;
}

ll sum(int n, int m)
{
int k = min(n, m);
ll res = 0;
for (int l = 1, r; l <= k; l = r + 1)
{
r = min(k, min(n / (n / l), m / (m / l)));
res = (res + (ll)(s[r] - s[l - 1] + MOD) * g(n / l, m / l) % MOD) % MOD;
}
return res;
}

int main()
{
init();
int n, m;
scanf("%d%d", &n, &m);
int k = min(n, m);
ll res = 0;
for (int l = 1, r; l <= k; l = r + 1)
{
r = min(k, min(n / (n / l), m / (m / l)));
res += (ll)(l + r) * (r - l + 1) / 2 * sum(n / l, m / l);
res %= MOD;
}
printf("%lld", res);
return 0;
}