luogu1829 Crash的数字表格

题目链接

problem

给出(n,m(n,mle10^7)),求(sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j))

(lcm(i,j))表示i和j的最小公倍数

solution

设(nle m)

[sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j)\ =sumlimits_{i=1}^nsumlimits_{j=1}^mfrac{ij}{gcd(i,j)}\ =sumlimits_{d=1}^ndsumlimits_{i=1}^{lfloor frac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}ij[gcd(i,j)=1]\ =sumlimits_{d=1}^ndsumlimits_{i=1}^{lfloor frac{n}{d} floor}sumlimits_{j=1}^{lfloorfrac{m}{d} floor}ijsumlimits_{k|i,k|j}mu(k)\ =sumlimits_{d=1}^ndsumlimits_{k=1}^{lfloorfrac{n}{d} floor}k^2mu(k)sumlimits_{i=1}^{lfloorfrac{n}{dk} floor}isumlimits_{j=1}^{lfloorfrac{m}{dx} floor}j ]

令(t=dx)
原式=(sumlimits_{t=1}^nsumlimits_{k|t}k^2mu(k)frac{t}{k}sumlimits_{i=1}^{lfloorfrac{n}{t} floor}isumlimits_{j=1}^{lfloorfrac{m}{t} floor}j)

发现后面的两个(sum)都可以(O(1))计算。然后就是如何处理前面(sumlimits_{k|t}k^2mu(k)frac{t}{k})的问题了。

显然(k^2mu(k))是积性函数，设(f(n)=n^2mu(n))。那么前面这一块其实就是(f*Id (k))。因为积性函数卷积性函数还是积性函数。所以前面这一块就是一个积性函数。线性筛即可。

那么这个函数到底该怎么筛呢。

按照套路，设(g=f*Id)先观察(g(q^p))的值，发现(g(q^p)=q^p-q^{p+1})。

所以筛的方式与筛(varphi)类似。

然后就可以(O(n))做了。

其实发现上式可以数论分块，那么瓶颈其实在预处理。所以此题可以出成多次询问的版本。

code

#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
#include<ctime>
#include<cmath>
using namespace std;
typedef long long ll;
const int N = 1e7 + 5,mod = 20101009;
ll read() {
	ll x = 0,f = 1;char c = getchar();
	while(c < '0' || c > '9') {
		if(c == '-') f = -1; c = getchar();
	}
	while(c >= '0' && c <= '9') {
		x = x * 10 + c - '0'; c = getchar();
	}
	return x * f;
}
int vis[N],tot,pri[N];
ll f[N];
inline ll calc(ll x) {
	return (x * (x + 1) / 2) % mod;
}
int main() {
	ll n = read(),m = read();
	if(n > m) swap(n,m);
	f[1] = 1;
	for(int i = 2;i <= n;++i) {
		if(!vis[i]) { pri[++tot] = i;f[i] = (i - 1ll * i * i) % mod; }
		
		for(int j = 1;j <= tot && pri[j] * i <= n;++j) {
			vis[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				f[i * pri[j]] = 1ll * f[i] * pri[j] % mod;
				break;
			}
			f[i * pri[j]] = f[i] * f[pri[j]];
		}
	}
	ll ans = 0;
	for(int i = 1;i <= n;++i) {
		ans += f[i] * calc(n / i) % mod * calc(m / i) % mod;
		ans %= mod;
	}
	cout<<(ans + mod) % mod;
	return 0;
}