题解P3768 简单的数学题
非常典型的一道杜教筛求前缀和题目
首先我们先顺着题目给出的式子往下面推
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^nijgcd(i,j)=\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^nijd[gcd(i,j)=d]\\=\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}ij[gcd(i,j)=1]
\]
用莫比乌斯反演,令 \(g(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\)
\[\Rightarrow ans=\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\frac{n}{d}}\mu(i)i^2g(\frac{n}{id})
\]
令 \(T=id\),尝试交换循环顺序
\[ans=\sum\limits_{T=1}^ng(\frac{n}{T})T^2\sum\limits_{d|T}d\mu(\frac{T}{d})
\]
将 \(\sum\limits_{d|T}d\mu(\frac{T}{d})\) 写成迪利克雷卷积的形式 \(\mu*id=\phi\)
\[\Rightarrow ans=\sum\limits_{T=1}^ng(\frac{n}{T})T^2\phi(T)
\]
令 \(f(n)=n^2\phi(n)\)
\[ans=\sum\limits_{T=1}^ng(\frac{n}{T})f(T)
\]
对于 \(g\) 可以使用数论分块,考虑如何快速求 \(S(n)=\sum\limits_{i=1}^nf(i)\)
首先搬出杜教筛式子 \(g(1)S(n)=\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)\)
\[(f*g)(n)=\sum\limits_{d|n}g(d)f(\frac{n}{d})=\sum\limits_{d|n}\frac{n^2}{d^2}\phi(\frac{n}{d})g(d)
\]
不妨令 \(g(i)=i^2\)
\[(f*g)(n)=\sum\limits_{d|n}n^2\phi(\frac{n}{d})=n^2\sum\limits_{d|n}\phi(\frac{n}{d})=n^3
\]
\[\Rightarrow S(n)=\sum\limits_{i=1}^ni^3-\sum\limits_{d=2}^ng(d)S(\lfloor\frac{n}{d}\rfloor)
\]
\[ans=\sum\limits_{T=1}^ng(\frac{n}{T})T^2\phi(T)
\]
前面那坨可以数论分块,后面那一坨可以杜教筛 \(O(n^{\frac{2}{3}})\) 处理前缀和,注意取模
#include<iostream>
#include<cstdio>
#include<unordered_map>
using namespace std;
typedef long long ll;
const int N = 8000001;
bool st[N];
int prime[N], tot;
ll sum[N], p, n, phi[N], inv2, inv6;
unordered_map<ll, ll> S;
void init()
{
phi[1] = 1;
for (int i = 2; i < N; i++)
{
if (!st[i])
prime[++tot] = i, phi[i] = i - 1;
for (int j = 1; i * prime[j] < N; j++)
{
st[i * prime[j]] = true;
if (i % prime[j] == 0)
{
phi[i * prime[j]] = (phi[i] * prime[j]) % p;
break;
}
phi[i * prime[j]] = (phi[i] * (prime[j] - 1)) % p;
}
}
for (int i = 1; i < N; i++)
sum[i] = (sum[i - 1] + (ll)phi[i] * i % p * i % p) % p;
}
ll Sum1(ll x)
{
ll res = (((x + 1) % p * x % p) % p) * inv2 % p;
return res * res % p;
}
ll Sum2(ll x)
{
ll r1 = (x % p * (x + 1) % p) % p;
return ((r1 * ((x + x + 1) % p)) % p * inv6) % p;
}
ll Djf(ll x)
{
if (x < N)
return sum[x];
if (S[x])
return S[x];
ll res = Sum1(x);
for (ll l = 2, r; l <= x; l = r + 1)
{
r = min(x, (x / (x / l)));
res -= (Sum2(r) - Sum2(l - 1))* (Djf(x / l)) % p;
res %= p;
}
return S[x] = (res + p) % p;
}
ll qmi(ll a, ll b)
{
ll res = 1;
while (b)
{
if (b & 1)
res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}
int main()
{
scanf("%lld%lld", &p, &n);
init();
inv2 = qmi(2, p - 2), inv6 = qmi(6, p - 2);
ll ans = 0;
for (ll l = 1, r; l <= n; l = r + 1)
{
r = min(n, n / (n / l));
ans = (ans + (Sum1(n / l) % p * ((Djf(r) - Djf(l - 1)) % p)) % p + p) % p;
}
printf("%lld", ans);
return 0;
}
