P3768 简单的数学题 题解

P3768 简单的数学题

题意：

输入一个整数 \(n\) 和一个整数 \(p\)，你需要求出：

\[\left(\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\right) \bmod p \]

思路：

此处可跳过

---

\(\displaystyle\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)=\sum_{i=1}^n\sum_{j=1}^n ij \sum_{d=1}^{n}[\gcd(i,j)=d]d\)
\(\displaystyle =\sum_{d=1}^{n}d^3\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor} ij [\gcd(i,j)=1]\)
令 \(\displaystyle G(n)=\sum_{i=1}^{n}\sum_{j=1}^{n} ij [\gcd(i,j)=1]\)，
\(\displaystyle G(n)=\sum_{d=1}^{n}d^2\cdot \mu(d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}i\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor} j\)
然后可以用杜教筛实现，不过这个有两次分块，最终时间复杂度为 \(O(n^{\frac{3}{4}}\log n)\)，会超时。

---

换种写法：
\(\displaystyle\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)=\sum_{i=1}^n\sum_{j=1}^n ij \sum_{d=1}^{n}[d|i][d|j]\varphi(d)\)
\(\displaystyle =\sum_{d=1}^{n}\varphi(d)\cdot d^2(\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}i)^2\)
用杜教筛实现，最终时间复杂度为 \(O(n^{\frac{2}{3}}\log n)\)

代码

点击查看代码

#include <bits/stdc++.h>
#define int long long
#define ll long long
#define ull unsigned long long
#define double long double
#define PII pair<int,int>
using namespace std;
mt19937_64 rng(time(0));

const int N = 4e6 + 10, M = 2e3 + 10;
int mod = 1e9 + 7;
const double EPS = 1e-8;


ll qpow(ll x, int n) {
	ll ans = 1;
	while (n) {
		if (n & 1) {
			ans *= x;
			ans %= mod;
		}
		x *= x;
		x %= mod;
		n >>= 1;
	}
	return ans;
}
//int fac[N];
//int inv[N];
//void init(){
//	fac[0]=inv[0]=1;
//	for(int i=1;i<N;i++){
//		fac[i]=fac[i-1]*i%mod;
//	}
//	inv[N-1]=qpow(fac[N-1],mod-2);
//	for(int i=N-2;i>=0;i--){
//		inv[i]=inv[i+1]*(i+1)%mod;
//	}
//}
//
//
//int C(int n,int m){
//	return fac[n]*inv[m]%mod*inv[n-m]%mod;
//}
constexpr int MAXN = 5e6 + 10;
long long T, n, pri[MAXN], cur, mu[MAXN], sum_mu[MAXN];
bool vis[MAXN];
map<long long, long long> mp_mu;

int inv6, inv2;
int s2(int x) {
	return (x % mod) * ((x + 1) % mod) % mod * ((2 * x + 1) % mod) % mod * inv6 % mod;
}
int s3(int x) {
	return (x % mod * ((x + 1) % mod) % mod * inv2 % mod) * (x % mod * ((x + 1) % mod) % mod * inv2 % mod) % mod;
}
long long S_mu(long long x) {  // 求mu的前缀和
	if (x < MAXN) return sum_mu[x];
	if (mp_mu[x]) return mp_mu[x];  // 如果map中已有该大小的mu值，则可直接返回
	long long ret = (long long)s3(x);
	for (long long i = 2, j; i <= x; i = j + 1) {
		j = x / (x / i);
		ret += mod - S_mu(x / i) * (s2(j) - s2(i - 1) + mod) % mod;
		ret %= mod;
	}
	return mp_mu[x] = ret;  // 路径压缩，方便下次计算
}


void init() {
	inv6 = qpow(6, mod - 2);
	inv2 = qpow(2, mod - 2);
	mu[1] = 1;
	for (int i = 2; i < MAXN; i++) {
		if (!vis[i]) {
			pri[++cur] = i;
			mu[i] = (i - 1) * i % mod * i % mod;

		}
		for (int j = 1; j <= cur && i * pri[j] < MAXN; j++) {
			vis[i * pri[j]] = true;
			if (i % pri[j])
				mu[i * pri[j]] = (mu[i] * mu[pri[j]]) % mod ;
			else {
				mu[i * pri[j]] = (mu[i] * pri[j]) % mod * ( pri[j] % mod) % mod * ( pri[j] % mod) % mod;
				break;
			}
		}
	}
	for (int i = 1; i < MAXN; i++)
		sum_mu[i] = (sum_mu[i - 1] + mu[i]) % mod; // 求mu数组前缀和
}

int H(int n) {
	int res = 0;
	int l = 1, r;
	while (l <= n) {
		r = n / (n / l);
		int t = (n / l % mod) * ((n / l + 1) % mod) % mod * inv2 % mod;
		t *= t;
		t %= mod;
		res += (S_mu(r) - S_mu(l - 1) + mod) % mod * t % mod;
		res %= mod;
		l = r + 1;
	}
	return res;
}



void solve() {

	cin >> mod >> n;
	init();

	cout << H(n) << endl;

}


signed main() {
	ios::sync_with_stdio(false);
	cin.tie(0);
	int t = 1;
//	cin >> t;

	while (t--) {
		solve();
	}

	return 0;
}