LG P3768 簡單的數學題
阿新 • • 發佈:2021-07-18
\(\text{Problem}\)
求
\[\left(\sum_{i=1}^n \sum_{j=1}^n i j \gcd(i,j)\right) \bmod p \]\(n \le 10^{10},5 \times 10^8 \le p \le 1.1 \times 10^9\) 且 \(p \in \mathbb{P}\)
\(\text{Solution}\)
顯然走尤拉反演
\[\begin{aligned} \sum_{i=1}^n \sum_{j=1}^n i j \gcd(i,j) &= \sum_{i=1}^n \sum_{j=1}^n i j \sum_{d|\gcd(i,j)} \varphi(d) \\ &= \sum_{d=1}^n \varphi(d) \sum_{d|i} i \sum_{d|j} j \\ &= \sum_{d=1}^n \varphi(d) \sum_{i=1}^{\lfloor \frac n d \rfloor} i \sum_{j=1}^{\lfloor \frac n d \rfloor} j \\ &= \sum_{d=1}^n d^2 \varphi(d) S^3 (\lfloor \frac n d \rfloor) \end{aligned} \]\(\varphi\)
這個式子顯然可以數論分塊(非常顯然)
那麼重點就是求 \(S(n)=\sum_{d=1}^n d^2 \varphi(d)\)
杜教篩即可
即考慮卷積 \(f * g\),記 \(f(n)=n^2 \varphi(n)\),令 \(g = ID^2\)
\[(f * g)(n) = \sum_{d|n} f(d) g(\frac{n}{d}) = \sum_{d|n} d^2 \varphi(d) \left(\frac{n}{d}\right)^2 = n^2 \sum_{d|n} \varphi(d) = n^3 \]那麼
仍然可以數論分塊,利用平方和與立方和公式快速計算
\(\text{Code}\)
#include<cstdio> #include<tr1/unordered_map> #define LL long long #define maxn 5000000 #define N 5000005 using namespace std; int vis[N], phi[N], prime[N], totp; LL P, n, inv6, sf[N]; tr1::unordered_map<LL, LL> SF; inline LL fpow(LL x, LL y) { LL res = 1; for(; y; y >>= 1) { if (y & 1) res = res * x % P; x = x * x % P; } return res; } inline LL S2(LL n) { n %= P; return n * (n + 1) % P * (n * 2 + 1) % P * inv6 % P; } inline LL S3(LL n) { n %= P; return n * (n + 1) / 2 % P * (n * (n + 1) / 2 % P) % P; } inline void sieve() { vis[1] = phi[1] = 1; for(register int i = 2; i <= maxn; i++) { if (!vis[i]) prime[++totp] = i, phi[i] = i - 1; for(register int j = 1; j <= totp && prime[j] * i <= maxn; j++) { vis[i * prime[j]] = 1; if (i % prime[j]) phi[i * prime[j]] = phi[i] * phi[prime[j]]; else{phi[i * prime[j]] = phi[i] * prime[j]; break;} } } for(register int i = 1; i <= maxn; i++) sf[i] = (sf[i - 1] + (LL)i * i % P * phi[i] % P) % P; } LL SumF(LL n) { if (n <= maxn) return sf[n]; if (SF[n]) return SF[n]; LL res = S3(n), r; for(register LL l = 2; l <= n; l = r + 1) { r = n / (n / l); res = (res - (S2(r) - S2(l - 1) + P) % P * SumF(n / l) % P + P) % P; } return SF[n] = res; } int main() { scanf("%lld%lld", &P, &n); sieve(); LL ans = 0, r; inv6 = fpow(6, P - 2); for(register LL l = 1; l <= n; l = r + 1) { r = n / (n / l); ans = (ans + S3(n / l) * (SumF(r) - SumF(l - 1) + P) % P) % P; } printf("%lld\n", ans); }