Luogu P6825 「EZEC-4」求和

题目描述

$T$ 组数据，对于每组数据，给定 $n,p$，计算下面的式子对 $p$ 取余的结果。

$$\sum_{i=1}^n\sum_{i=1}^n\gcd(i, j)^{i+j}$$

其中 $1 \leq T \leq 2,1 \leq \sum n \leq 1.5 \times 10^6, 2 \leq p \leq 2^{31}-1 \text{且 p 为质数}$。

$\mathcal{Solution}$

推一下式子：

$$\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^{i+j} \\ \sum_{d=1}^n\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor } \sum_{j=1}^{\left \lfloor \frac{n}{d} \right \rfloor } [\gcd(i,j)=1]d^{id+jd}$$

利用 $\mu*1=\epsilon$：

$$\sum_{d=1}^n\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor } \sum_{j=1}^{\left \lfloor \frac{n}{d} \right \rfloor } \sum_{k|i,k|j} \left ( \mu(k)\times d^{id+jd} \right )$$

把枚举 $k$ 提到前面，原先的枚举 $i,j$ 就变成了枚举 $k$ 的倍数 $ik,jk$。

$$\sum_{d=1}^n\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor} \mu(k)\sum_{i=1}^{\left \lfloor \frac{n}{kd} \right \rfloor } \sum_{j=1}^{\left \lfloor \frac{n}{kd} \right \rfloor } d^{ikd+jkd}$$

设 $T=kd$：

$$\sum_{d=1}^n\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor } \mu(k) \sum_{i=1}^{\left \lfloor \frac{n}{T} \right \rfloor } \sum_{j=1}^{\left \lfloor \frac{n}{T} \right \rfloor } d^{iT}\times d^{jT} \\ \sum_{d=1}^n\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor } \mu(k) \left ( \sum_{i=1}^{\left \lfloor \frac{n}{T} \right \rfloor }d^{iT} \right ) \left ( \sum_{j=1}^{\left \lfloor \frac{n}{T} \right \rfloor } d^{jT} \right ) \\ \sum_{d=1}^n\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor } \mu(k) \left ( \sum_{i=1}^{\left \lfloor \frac{n}{T} \right \rfloor }{(d^{T})}^i \right )^2$$

平方里面是个等比数列求和，可以在 $\mathcal{O}(\log \left \lfloor \frac{n}{T} \right \rfloor)$ 的时间复杂度内解决，至于 $d^T$，可以在每枚举到一个 $d$ 的时候快速幂一下 $d^d$，然后在枚举 $k$ 一步一步乘。

枚举 $d,k$ 是调和级数 $\mathcal{O}(n\log n)$，后面复杂度就不会分析了，有大佬证明出是 $\mathcal{O}(n\log n)$ 的。

如何在 $\mathcal{O}(\log n)$ 的时间内算 $n$ 项等比数列的和？

设 $S_n$ 为等比数列的前 $n$ 项和，那么有：

$$S_n=\left\{\begin{matrix} S_{\frac{n}{2}}(1+q^{\frac{n}{2}}),n\bmod2=0 \\ S_{\frac{n-1}{2}}(1+q^{\frac{n-1}{2}})+q^n,n\bmod2=1 \end{matrix}\right.$$

实际上就是在分治快速幂算 $q^n$ 的时候顺带算了一下 $S_n$，注意到在算 $S_{n}$ 的时候把 $q^n$ 算出来后传给上一层，这样复杂度就为 $\mathcal{O}(\log n)$ 了，不必再快速幂计算 $q^n$。

$\mathcal{Code}$

#include<iostream>
#include<cstdio>
#include<algorithm>
namespace do_while_true {
	#define ld double
	#define ll long long
	#define re register
	#define pb push_back
	#define fir first
	#define sec second
	#define pp std::pair<ll, ll>
	#define mp std::make_pair
	template <typename T>
	inline T Max(T x, T y) { return x > y ? x : y; }
	template <typename T>
	inline T Min(T x, T y) { return x < y ? x : y; }
	template <typename T>
	inline T Abs(T x) {	return x < 0 ? -x : x; }
	template <typename T>
	inline T& read(T& r) {
		r = 0; bool w = 0; char ch = getchar();
		while(ch < '0' || ch > '9') w = ch == '-' ? 1 : 0, ch = getchar();
		while(ch >= '0' && ch <= '9') r = r * 10 + (ch ^ 48), ch = getchar();
		return r = w ? -r : r;
	}
	template <typename T>
	inline T qpow(T x, T y, T mod) {
		re T sumq = 1; x %= mod;
		while(y) {
			if(y&1) sumq = sumq * x % mod;
			x = x * x % mod;
			y >>= 1;
		}
		return sumq;
	}
	char outch[110];
	int outct;
	template <typename T>
	inline void print(T x) {
		do {
			outch[++outct] = x % 10 + '0';
			x /= 10;
		} while(x);
		while(outct >= 1) putchar(outch[outct--]);
	}
}
using namespace do_while_true;

const int N = 1500010 + 10;
template <typename T, typename T2>
inline T Mod(T x, T2 p) { return x % p; }

ll ans;
bool vis[N];
int prime[N], ct, mu[N];

pp db(ll x, int n, ll p) {
	if(n == 0) return mp(1ll, 1ll);
	if(n == 1) return mp(x, x);
	if((n & 1) == 0) {
		pp temp = db(x, n / 2, p);
		return mp(temp.fir * ((1 + temp.sec) % p) % p, temp.sec * temp.sec % p);
	}
	pp temp = db(x, (n-1) / 2, p);
	ll qn = temp.sec * temp.sec % p * x % p;
	return mp((temp.fir * ((1 + temp.sec) % p) % p + qn) % p, qn);
}

void mian(int n, ll p) {
	ans = 0;
	for(int d = 1; d <= n; ++d) {
		ll dd = qpow(1ll * d, 1ll * d, p);
		ll kd = dd;
		for(int k = 1; k <= n / d; ++k) {
			pp temp = db(kd, n / (k * d), p);
			ans = (ans + temp.fir * temp.fir % p * Mod(mu[k] + p, p) % p) % p;
			kd = kd * dd % p;
		}
	}
	printf("%lld\n", ans);
}

int n[3], mx;
ll p[3];

void pre() {
	vis[1] = 1; mu[1] = 1;
	for(int i = 2; i <= mx; ++i) {
		if(!vis[i]) {
			prime[++ct] = i;
			mu[i] = -1;
		}
		for(int j = 1; j <= ct && 1ll * i * prime[j] <= mx; ++j) {
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) {
				mu[i * prime[j]] = 0;
				break;
			}
			mu[i * prime[j]] = mu[i] * mu[prime[j]];
		}
	}
}

signed main() {
	int T = 1;
	read(T);
	for(int i = 1; i <= T; ++i) read(n[i]), read(p[i]), mx = Max(mx, n[i]);
	pre();
	for(int i = 1; i <= T; ++i) mian(n[i], p[i]);
	fclose(stdin);
	return 0;
}