BZOJ 2154 Crash的数字表格 (莫比乌斯反演)

Crash的数字表格

今天的数学课上，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。

题解

n < m
Ans=∑i=1n∑j=1mlcm(i,j)=∑i=1n∑j=1mijgcd(i,j)=∑d=1n∑i=1n∑j=1mijd[gcd(i,j)==d]=∑d=1n∑i=1n/d∑j=1m/dijd[gcd(i,j)==1]=∑d=1nd∑i=1n/d∑j=1m/dij[gcd(i,j)==1]令:Sum(n,m)=∑i=1n∑j=1mijf(x)=∑i=1n∑j=1mij[gcd(i,j)==x]g(x)=∑i=1n∑j=1mij[x∣gcd(i,j)]=x2∑i=1n/x∑j=1m/xij[1∣gcd(i,j)]=∑x∣dnf(d)则：f(x)=∑x∣dnμ(dx)g(d)f(1)=∑d=1nμ(d)d2Sum(nd,md)Ans=∑d=1nd∑i=1n/d∑j=1m/dij[gcd(i,j)==1]=∑d=1nd∑i=1n/dμ(i)i2Sum(nid,mid)\begin{aligned} Ans &= {\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}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{d}[gcd(i, j)==d]}\\ &={\sum_{d=1}^{n}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{ijd}[gcd(i, j)==1]}\\ &={\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{ij}[gcd(i, j)==1]}\\ 令:Sum(n,m) &={\sum_{i=1}^{n}\sum_{j=1}^{m}ij}\\ f(x) &= {\sum_{i=1}^{n}\sum_{j=1}^{m}ij[gcd(i,j)==x]}\\ g(x)&={\sum_{i=1}^{n}\sum_{j=1}^{m}ij[x|gcd(i,j)]}\\ &={x^2\sum_{i=1}^{n/x}\sum_{j=1}^{m/x}ij[1|gcd(i,j)]}\\ &={\sum_{x|d}^{n}f(d)}\\ 则：f(x)&={\sum_{x|d}^{n}\mu(\frac{d}{x})g(d)}\\ f(1)&={\sum_{d=1}^{n}\mu(d)d^2{Sum(\frac{n}{d},\frac{m}{d})}}\\ Ans &={\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{ij}[gcd(i, j)==1]}\\ &=\sum_{d=1}^{n}d{\sum_{i=1}^{n/d}\mu(i)i^2{Sum(\frac{n}{id},\frac{m}{id})}}\\ \end{aligned} Ans令:Sum(n,m)f(x)g(x)则：f(x)f(1)Ans=i=1∑nj=1∑mlcm(i,j)=i=1∑nj=1∑mgcd(i,j)ij=d=1∑ni=1∑nj=1∑mdij[gcd(i,j)==d]=d=1∑ni=1∑n/dj=1∑m/dijd[gcd(i,j)==1]=d=1∑ndi=1∑n/dj=1∑m/dij[gcd(i,j)==1]=i=1∑nj=1∑mij=i=1∑nj=1∑mij[gcd(i,j)==x]=i=1∑nj=1∑mij[x∣gcd(i,j)]=x2i=1∑n/xj=1∑m/xij[1∣gcd(i,j)]=x∣d∑nf(d)=x∣d∑nμ(xd)g(d)=d=1∑nμ(d)d2Sum(dn,dm)=d=1∑ndi=1∑n/dj=1∑m/dij[gcd(i,j)==1]=d=1∑ndi=1∑n/dμ(i)i2Sum(idn,idm)

μ(i)i2\mu(i)i^2μ(i)i2用前缀和预处理
求f(x}的时候用分块可以在n{\sqrt{n}}n完成
求Sum(x,y)的时候用分块也可以在n{\sqrt{n}}n完成
所以一次询问可以在O(n)的时间解决
如果是多次询问（10000次询问）
需要继续化简
令T=id，考虑对T进行分块Ans=∑d=1nd∑i=1n/dμ(i)i2Sum(nT,mT)=∑T=1nSum(nT,mT)∑i∣Tnμ(i)i2Ti\begin{aligned} 令T&=id，考虑对T进行分块\\ Ans &=\sum_{d=1}^{n}d{\sum_{i=1}^{n/d}\mu(i)i^2{Sum(\frac{n}{T},\frac{m}{T})}}\\ &={\sum_{T=1}^{n}Sum(\frac{n}{T}, \frac{m}{T})\sum_{i|T}^{n}\mu(i)i^2\frac{T}{i}} \end{aligned} 令TAns=id，考虑对T进行分块=d=1∑ndi=1∑n/dμ(i)i2Sum(Tn,Tm)=T=1∑nSum(Tn,Tm)i∣T∑nμ(i)i2iT
∑i∣Tnμ(i)i2Ti\sum_{i|T}^{n}\mu(i)i^2\frac{T}{i}∑i∣Tnμ(i)i2iT是积性函数可以用线性筛预处理，令t = i × prime[j]
- 线性筛的时候如果i % prime[j] = 0，则t对应的莫比乌斯函数为0，prime[j]的贡献为g[i] * prime[j]
- i % prime[j] != 0，符合积形函数直接相乘
- i 为素数，g[i]=i−i2g[i] = i - i^2g[i]=i−i2
  Sum(x,y)可以在n\sqrt{n}n完成

// 单次询问 O(n)
#include <bits/stdc++.h>
#define LL  long long
#define P pair<int, int>
#define lowbit(x) (x & -x)
#define mem(a, b) memset(a, b, sizeof(a))
#define rep(i, a, n) for (int i = a; i <= n; ++i)
#define maxn 10000006
#define mid ((l + r) >> 1)
#define lc rt<<1
#define rc rt<<1|1
using namespace std;int mod = 20101009;
int mo[maxn], prime[maxn], g[maxn];
bool vis[maxn];
int Sum(int x) {return (1ll * (1 + x) * x / 2) % mod;
}
void init() {mo[1] = 1;vis[1] = 1;int len = 0;for (int i = 2; i < maxn; ++i) {if (!vis[i]) {prime[len++] = i;mo[i] = -1;}for (int j = 0; j < len && 1ll * i * prime[j] < maxn; ++j) {int t = i * prime[j];vis[t] = 1;if (i % prime[j] == 0) {mo[t] = 0;break;}mo[t] = -mo[i];}}
}int solve(int n, int m) {int i = 1, j;int sum = 0;while (i <= n) {j = min(n / (n/i), m / (m/i));sum += (1ll * (g[j] - g[i-1]) * Sum(n/i) % mod) * Sum(m/i) % mod;sum %= mod;i = j + 1;}return sum;
}int main() {#ifndef ONLINE_JUDGE// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);
#endifinit();for (int i = 1; i < maxn; ++i) {g[i] = 1ll * i * i * mo[i] % mod;g[i] = (g[i] + g[i-1]) % mod;}int n, m;scanf("%d %d", &n, &m);if (n > m)    swap(n, m);int i = 1, j;int ans = 0;while (i <= n) {j = min(n/(n/i), m/(m/i));int t = 1ll * (i + j) * (j - i + 1) / 2 % mod;ans += 1ll * solve(n/i, m/i) * t % mod;ans %= mod;i = j + 1;}printf("%d\n", (ans + mod) % mod);return 0;}

// 单次询问O(sqrt(n))
#include <bits/stdc++.h>
#define LL  long long
#define P pair<int, int>
#define lowbit(x) (x & -x)
#define mem(a, b) memset(a, b, sizeof(a))
#define rep(i, a, n) for (int i = a; i <= n; ++i)
#define maxn 10000006
#define mid ((l + r) >> 1)
#define lc rt<<1
#define rc rt<<1|1
using namespace std;int mod = 20101009;
int mo[maxn], prime[maxn], g[maxn];
bool vis[maxn];
int Sum(int x) {return (1ll * (1 + x) * x / 2) % mod;
}
void init() {// mo[1] = 1;g[1] = 1;vis[1] = 1;int len = 0;for (int i = 2; i < maxn; ++i) {if (!vis[i]) {prime[len++] = i;// mo[i] = -1;g[i] = (i - 1ll * i * i % mod) % mod;}for (int j = 0; j < len && 1ll * i * prime[j] < maxn; ++j) {int t = i * prime[j];vis[i * prime[j]] = 1;if (i % prime[j] == 0) {// mo[t] = 0;g[t] = 1ll * g[i] * prime[j] % mod;break;}// mo[t] = -mo[i];g[t] = 1ll * g[i] * g[prime[j]] % mod;}}for (int i = 2; i < maxn; ++i) {g[i] = (g[i] + g[i-1]) % mod;}
}int main() {#ifndef ONLINE_JUDGE// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);
#endifinit();int n, m;scanf("%d %d", &n, &m);if (n > m)    swap(n, m);int i = 1, j;int ans = 0;while (i <= n) {j = min(n/(n/i), m/(m/i));ans += (1ll * Sum(n/i) * Sum(m/i) % mod) * (g[j] - g[i-1]) % mod;ans %= mod;i = j + 1;}printf("%d\n", (ans + mod) % mod);return 0;
}