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;
}
