题目链接
题目描述
输入一个整数n和一个整数p，你需要求出
[sum_{i=1}^nsum_{j=1}^n (icdot jcdot gcd(i,j)) mod p
]
其中(gcd(a,b))表示(a)与(b)的最大公约数
输入
一行两个整数(p,n)
输出
一行一个整数，为题目中所求值
样例
样例输入
998244353 2000
样例输出
883968974
数据范围
(nleq 10^{10})
(5 imes 10^8 leq p leq 1.1 imes 10^9)
(p)为质数（但貌似也可以不是？又不用求逆元）
题解
自己想出来的题！但是连(WA)两发就是因为杜教筛写挂了……
先不考虑取余，我们化一下题目中的式子，枚举(gcd)（警告！多公式）。
[sum_{i=1}^nsum_{j=1}^n icdot jcdot gcd(i,j)
]
[sum_{d=1}^{n}dsum_{i=1}^{leftlfloor frac{n}{d}
ight
floor}sum_{j=1}^{leftlfloor frac{n}{d}
ight
floor}[iperp j]icdot j cdot d^2
]
[sum_{d=1}^{n}d^3sum_{i=1}^{left lfloor frac{n}{d}
ight
floor}sum_{j=1}^{leftlfloor frac{n}{d}
ight
floor}[iperp j]icdot j
]
[sum_{d=1}^{n}d^3sum_{p=1}^{left lfloor frac{n}{d}
ight
floor}mu(p)p^2cdot Big(frac{(1+leftlfloor frac{n}{dp}
ight
floor)leftlfloor frac{n}{dp}
ight
floor}{2}Big)^2
]
额，现在可以使用分块优化做到(O(n))了，但是这完全不能胜任数据范围，我们换个角度，设(dp=T)，枚举(T)会有什么结果？
[sum_{T=1}^{n}Big(frac{(1+leftlfloor frac{n}{T}
ight
floor)leftlfloor frac{n}{T}
ight
floor}{2}Big)^2sum_{d|T}d^3cdot mu(frac{T}{d})(frac{T}{d})^2
]
[sum_{T=1}^{n}Big(frac{(1+leftlfloor frac{n}{T}
ight
floor)leftlfloor frac{n}{T}
ight
floor}{2}Big)^2 T^2sum_{d|T}dcdot mu(frac{T}{d})
]
现在好像反而变成(O(nlog n))或(O(nsqrt{n}))了，别急，我们看看第二层的求和的意义——狄利克雷卷积，这是(Id)函数与(mu)函数的狄利克雷卷积，其值就等于(varphi)。
[sum_{T=1}^{n}Big(frac{(1+leftlfloor frac{n}{T}
ight
floor)leftlfloor frac{n}{T}
ight
floor}{2}Big)^2 T^2varphi(T)
]
现在，我们只需要快速求出一个东西即可——(T^2varphi(T))，前面的部分可以分块优化，我们急需解决的就是这个函数(f(T)=T^2varphi(T))的前缀和(F(T))。显然，这是一个积性函数。
杜教筛的公式：
[sum_{i=1}^{n}(f*g)(i)=sum_{i=1}^{n}sum_{d|i}f(d)cdot g(frac{i}{d})=sum_{i=1}^{n}g(i)sum_{j=1}^{leftlfloor frac{n}{i}
ight
floor}f(j)
]
于是我们需要一个函数与(f)卷起来，我们根据套路或枚举发现(T^2)项很恼人，于是尝试把这一项消掉，于是想到了(g(x)=x^2)。
[sum_{i=1}^{n}sum_{d|i}d^2varphi(d)cdot (frac{i}{d})^2=sum_{i=1}^{n}i^2sum_{j=1}^{leftlfloor frac{n}{i}
ight
floor}f(j)
]
[sum_{i=1}^{n}i^2sum_{d|i}varphi(d)=sum_{i=1}^{n}i^2F(leftlfloor frac{n}{i}
ight
floor)
]
根据公式(sum_{d|i}varphi(d)=i)，继续变形
[sum_{i=1}^{n}i^3=F(n)+sum_{i=2}^{n}i^2F(leftlfloor frac{n}{i}
ight
floor)
]
[F(n)=sum_{i=1}^{n}i^3-sum_{i=2}^{n}i^2F(leftlfloor frac{n}{i}
ight
floor)
]
由于(p(i)=i^3)和(q(i)=i^2)的前缀和都有公式，我们可以对右边进行分块优化，就可以杜教筛了！这道题圆满解决，时间复杂度(O(n^{frac{2}{3}}))。
不过有些小细节要注意，比如模数乘(2)可能会爆(int)，(n^2)可能会爆(long long)，需要先取模再平方
(Code:)
#include <map>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define N 5000005
#define ll long long
map<ll, ll>Phi;
ll n, mod, g[N];
int p[N], h[N], phi[N], cnt;
ll sqr(ll x)
{
ll a = 2 * x + 1, b = x + 1, c = x;
if (b % 2 == 0)b /= 2;
else c /= 2;
if (a % 3 == 0)a /= 3;
else
if (b % 3 == 0)b /= 3;
else c /= 3;
a %= mod, b %= mod, c %= mod;
return a * b % mod * c % mod;
}
ll seq(ll x)
{
ll a = x + 1, b = x;
if (a % 2 == 0)a /= 2;
else b /= 2;
a %= mod, b %= mod;
return a * b % mod;
}
ll vas(ll x)
{
ll a = seq(x);
return a * a % mod;
}
ll G(ll x)
{
if (x <= N - 5)
return g[int(x)];
if (Phi.find(x) != Phi.end())
return Phi[x];
ll ans = vas(x);
ll lst = 1;
for (ll i = 2; i <= x; i++)
{
i = x / (x / i);
ll w = (sqr(i) - sqr(lst)) % mod;
ans = (ans - w * G(x / i) % mod) % mod;
lst = i;
}
if (ans < 0)
ans += mod;
Phi.insert(make_pair(x, ans));
return ans;
}
ll Ans(ll x)
{
ll ans = 0, lst = 0;
for (ll i = 1; i <= x; i++)
{
i = x / (x / i);
ll z = seq(x / i);
z = z * z % mod;
ans = (ans + z * (G(i) - G(lst)) % mod) % mod;
lst = i;
}
if (ans < 0)
ans += mod;
return ans;
}
int main()
{
phi[1] = 1;
for (int i = 2; i <= N - 5; i++)
{
if (!h[i])
{
phi[i] = i - 1;
p[++cnt] = i;
}
for (int j = 1; j <= cnt; j++)
{
if (i * p[j] > N - 5)
break;
h[i * p[j]] = 1;
if (i % p[j] == 0)
phi[i * p[j]] = phi[i] * p[j];
else
phi[i * p[j]] = phi[i] * (p[j] - 1);
}
}
cin >> mod >> n;
for (int i = 1; i <= N - 5; i++)
g[i] = (g[i - 1] + 1ll * phi[i] * i % mod * i % mod) % mod;
cout << Ans(n) << '
';
}