Bzoj2154: Crash的数字表格

题意

\(求ans=\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i, j)\)
n,m<=10^7

Sol

\(原式=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i*j}{gcd(i, j)}\)
假设n < m，
\(则ans=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} d*i*j*[gcd(i,j)==1]\)
枚举d，看里面的
\(令x=\lfloor\frac{n}{d}\rfloor, y=\lfloor\frac{m}{d}\rfloor\)
\(设f(k)=\sum_{i=1}^{x}\sum_{i=1}^{y}i*j*[gcd(i, j)==1]\)
\(设g(i)=\sum_{i|d} f(d)\)即表示gcd是i及i的倍数的数对的乘积和
\(就是i^2*\frac{(\lfloor\frac{x}{i}\rfloor + 1)*\lfloor\frac{x}{i}\rfloor}{2}*\frac{(\lfloor\frac{y}{i}\rfloor + 1)*\lfloor\frac{y}{i}\rfloor}{2}\)
然后就可以莫比乌斯反演求出f数组，从而得到答案
但是暴力求显然跑不过10^7的点
所以可以用两个数论分块+前缀和优化，记得取模，会爆

代码

不用数论分块，暴力

# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int _(1e7), MOD(20101009);
 
IL ll Read(){
    char c = '%'; ll x = 0, z = 1;
    for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1;
    for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0';
    return x * z;
}
 
int prime[_], num, mu[_];
bool isprime[_];
 
IL void Prepare(){
	mu[1] = 1;
	for(RG int i = 2; i <= _; ++i){
		if(!isprime[i]){  prime[++num] = i; mu[i] = -1;  }
		for(RG int j = 1; j <= num && i * prime[j] <= _; ++j){
			isprime[i * prime[j]] = 1;
			if(i % prime[j]) mu[i * prime[j]] = -mu[i];
			else{  mu[i * prime[j]] = 0; break;  }
		}
	}
}
 
IL ll Calc(RG ll n, RG ll m){
	RG ll f = 0, g;
	for(RG ll i = 1; i <= n; ++i){
		RG ll x = n / i, y = m / i;
		g = i * i % MOD * (x * (x + 1) >> 1) % MOD * (y * (y + 1) >> 1) % MOD;
		(f += 1LL * mu[i] * g % MOD) %= MOD;
	}
	return (f + MOD) % MOD;
}
 
int main(RG int argc, RG char *argv[]){
	Prepare();
	RG ll n = Read(), m = Read(), ans = 0;
	if(n > m) swap(n, m);
	for(RG ll d = 1; d <= n; ++d) (ans += d * Calc(n / d, m / d) % MOD) %= MOD;
	printf("%lld\n", ans);
	return 0;
}

用数论分块

# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int _(1e7 + 10), MOD(20101009);

IL ll Read(){
    char c = '%'; ll x = 0, z = 1;
    for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1;
    for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0';
    return x * z;
}

int prime[_], num, mu[_], squ[_], s[_];
bool isprime[_];

IL void Prepare(){
	mu[1] = 1; RG int maxn = 1e7;
	for(RG int i = 2; i <= 1e7; ++i){
		if(!isprime[i]){  prime[++num] = i; mu[i] = -1;  }
		for(RG int j = 1; j <= num && i * prime[j] <= maxn; ++j){
			isprime[i * prime[j]] = 1;
			if(i % prime[j]) mu[i * prime[j]] = -mu[i];
			else{  mu[i * prime[j]] = 0; break;  }
		}
	}
}

IL ll Calc(RG ll n, RG ll m){
	RG ll f = 0, g, j;
	for(RG ll i = 1; i <= n; i = j + 1){
		RG ll x = n / i, y = m / i; j = min(n / (n / i), m / (m / i));
		g = (x * (x + 1) >> 1) % MOD * ((y * (y + 1) >> 1) % MOD) % MOD;
		(f += 1LL * (squ[j] - squ[i - 1]) % MOD * g % MOD) %= MOD;
	}
	return (f + MOD) % MOD;
}

int main(RG int argc, RG char *argv[]){
	Prepare();
	RG ll n = Read(), m = Read(), ans = 0, j;
	if(n > m) swap(n, m);
	for(RG int i = 1; i <= n; ++i)
		squ[i] = ((squ[i - 1] + 1LL * mu[i] * i * i % MOD) % MOD + MOD) % MOD, s[i] = (s[i - 1] + i) % MOD;
	for(RG ll d = 1; d <= n; d = j + 1){
		j = min(n / (n / d), m / (m / d));
		(ans += 1LL * ((s[j] - s[d - 1]) % MOD + MOD) % MOD * Calc(n / d, m / d) % MOD) %= MOD;
	}
	printf("%lld\n", ans);
	return 0;
}