P1829 [國家集訓隊]Crash的數字表格 / JZPTAB

>題目連結<

為了書寫方便，本文中涉及的除法除特殊說明以外均向下取整

題目描述

今天的數學課上，Crash 小朋友學習了最小公倍數（Least Common Multiple）。對於兩個正整數 \(a\) 和 \(b\)，\(lcm(a,b)\) 表示能同時整除 \(a\) 和 \(b\) 的最小正整數。例如，\(lcm(6,8)=24\)。

回到家後，Crash 還在想著課上學的東西，為了研究最小公倍數，他畫了一張 \(n \times m\) 的表格。每個格子裡寫了一個數字，其中第 \(i\) 行第 \(j\) 列的那個格子裡寫著數為
\(lcm(i,j)\)

。

看著這個表格，Crash 想到了很多可以思考的問題。不過他最想解決的問題卻是一個十分簡單的問題：這個表格中所有數的和是多少。當 \(n\) 和 \(m\) 很大時，Crash 就束手無策了，因此他找到了聰明的你用程式幫他解決這個問題。由於最終結果可能會很大，Crash 只想知道表格裡所有數的和 \(\bmod 20101009\) 的值。

輸入格式

輸入包含一行兩個整數，分別表示 \(n\) 和 \(m\)。

輸出格式

輸出一個正整數，表示表格中所有數的和 \(\bmod 20101009\)的值。

輸入輸出樣例

輸入

4 5

輸出

122

---

解析

題目求

\[\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j) \]

不妨令 \(n\leq m\)，正確性顯然。

由小學奧數，我們可知：\(lcm(i,j)=\frac{ij}{gcd(i,j)}\)

所以題目求的是

\[\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{gcd(i,j)} \]

於是開始愉快(並不)推式子：令\(gcd(i,j)=d\)

\[\begin{aligned} & \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}[gcd(i,j)=d]\frac{ij}{d}\\ & =\sum_{d=1}^{n}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)=1]\frac{id\cdot jd}{d}\\ & =\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)=1]ij \end{aligned}\\ \]

後面這一團和GCD結論的形式很相似，我們可以依樣畫葫蘆推一下。

對於上界 \(n/d,m/d\) 不要想太多，可以當成 \(p,q\) 或其他變數推就是了。

\[\text{設：}f(x)=\sum_{i=1}^{p}\sum_{j=1}^{q}[gcd(i,j)=x]ij\\ F(x)=\sum_{x|d}^{p}f(d)\\ 由莫比烏斯反演得到:\\ f(x)=\sum_{x|d}\mu(\frac{d}{x})F(d)\\ f(1)=\sum_{i=1}^{p}\mu(i)F(i) \]

現在看 \(F(x)\) 如何求：

\[\begin{aligned} F(x) & =\sum_{x|d}^{p}f(d)\\ & =\sum_{x|d}^{p}\sum_{i=1}^{p}\sum_{j=1}^{q}[gcd(i,j)=d]ij\\ & =\sum_{i=1}^{p}\sum_{j=1}^{q}[x|gcd(i,j)]ij\\ & =x^2\sum_{i=1}^{p/x}\sum_{j=1}^{q/x}ij\\ & =x^2\sum_{i=1}^{p/x}i\cdot \frac{(1+\frac{q}{x})\frac{q}{x}}{2}\\ & =x^2\cdot \frac{(1+\frac{p}{x})\frac{p}{x}}{2}\cdot\frac{(1+\frac{q}{x})\frac{q}{x}}{2} \end{aligned} \]

代入莫比烏斯反演式子後：

\[Ans=\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\mu(i)\cdot i^2\cdot\frac{(1+\frac{n}{id})\frac{n}{id}}{2}\cdot\frac{(1+\frac{m}{id})\frac{m}{id}}{2}\\ \text{令}g(n,m)=\sum_{i=1}^{n}\mu(i)\cdot i^2\cdot\frac{(1+\frac{n}{i})\frac{n}{i}}{2}\cdot\frac{(1+\frac{m}{i})\frac{m}{i}}{2}\\ Ans=\sum_{d=1}^{n}d\cdot g(\frac{n}{d},\frac{m}{d}) \]

對於 \(g(\frac{n}{d},\frac{m}{d})\) ，我們可以數論分塊求。

而對於整個式子，也是可以數論分塊的。

總複雜度 \(O(n+m)\) 。

code:(三年OI一場空，不開long long、unsigned long long見祖宗)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const ll N=1e7+10,mod=20101009;

int primes[N],tot=0;
ll mu[N],sum[N];
bool mp[N];

void init(int n)
{
	mu[1]=mp[1]=1;
	for(int i=2; i<=n; i++)
	{
		if(!mp[i]) primes[++tot]=i,mu[i]=-1;
		for(int j=1; i*primes[j]<=n; j++)
		{
			int x=primes[j]*i;
			mp[x]=1;
			if(i%primes[j]==0)
			{
				mu[x]=0;
				break;
			}
			mu[x]=-mu[i];
		}
	}

	for(int i=1; i<=n; i++)
		sum[i]=(sum[i-1]+1LL*i*i%mod*(mu[i]+mod))%mod;
}

inline ll func(ll n,ll m)
{
	return (1LL*(n+1)*n/2%mod)*1LL*((m+1)*m/2%mod)%mod;
}

ll solve(ll n,ll m)
{
	if(n>m) swap(n,m);
	ll res=0;
	for(int i=1,j=0; i<=n; i=j+1)
	{
		j=min(n/(n/i),m/(m/i));
		res=(res+1LL*(sum[j]-sum[i-1]+mod)*func(n/i,m/i)%mod)%mod;
	}
	return res;
}

int main()
{
	init(1e7);
	ll n,m;
	scanf("%lld%lld",&n,&m);
	if(n>m) swap(n,m);
	ll ans=0;
	for(int i=1,j=0; i<=n; i=j+1)
	{
		j=min(n/(n/i),m/(m/i));
		ans=(ll)(ans+1LL*(j-i+1)*(i+j)/2%mod*solve(n/i,m/i)%mod)%mod;
	}
	printf("%lld",ans);
	return 0;
}