题意
求 $\sum_{i=a}^b \sum_{j=1}^i \frac{lcm(i,j)}{i}$.
分析
只需要求出前缀和，
$$\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^i \frac{lcm(i,j)}{i} &= \sum_{i=1}^n \sum_{j=1}^i \frac{j}{gcd(i,j)} \\
&= \sum_{d=1}^n \sum _{i=1}^n \sum_{j=1}^i \frac{j}{d} \cdot [gcd(i,j)=1] \\
&= \sum_{d=1}^n \sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor} \sum_{j=1}^i j \cdot [gcd(i,j)=1]
\end{aligned}$$
其后面部分提出来，即求 $\sum_{i=1}^n i\cdot [gcd(i,n)=1]$，对于这种一个值固定的gcd求和有一个套路，即倒序两两配对：
若 $n=1$，和为1；
若 $n>1$，因为 $gcd(i, n) = gcd(n-i, n)$ 且 $\displaystyle \sum_{i=1}^n i\cdot [gcd(i,n)=1] = \sum_{i=1}^{n-1} i\cdot [gcd(i,n)=1]$，
$\displaystyle \sum_{i=1}^{n-1} i \cdot [gcd(i,n)=1] + \sum_{i=n-1}^1i\cdot [gcd(i,n)=1] = n\varphi (n)$，所以和为 $ n\varphi (n) /2$.
综合得 $\displaystyle \sum_{i=1}^n i\cdot [gcd(i,n)=1] = \frac{n\varphi (n)+[n=1]}{2}$.
具体实现上，$[i=1]$ 只成立 $n$，除2可以提出来，即 原式 = $\displaystyle \frac{1}{2}(\sum_{d=1}^n \sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor} i\varphi (i) + n)$.
现在唯一得问题是如何求 $\displaystyle S(n) = \sum_{i=1}^n i\varphi (i)$.
根据杜教筛，
设 $\displaystyle S(n) = \sum_{i=1}^n f(i)$，$f(n) = n\varphi (n), \ g(n) = n$.
$\displaystyle f*g = \sum_{d|n} d \varphi (d) \cdot \frac{n}{d} = n\sum_{d|n}\varphi (d) = n^2$.
因此 $\displaystyle S(n) = \sum_{i=1}^n i^2 - \sum_{d=2}^n d\cdot S(\left \lfloor \frac{n}{d} \right \rfloor)$.
再最外层套个数论分块即可。
代码
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <map>
#include<unordered_map>
using namespace std;
const int maxn = ;
typedef long long ll;
const ll mod = ;
const ll inv2 = (mod+)>>;
const ll inv6 = ; //
ll T, a, b, pri[maxn], tot, phi[maxn], sum_phi_d[maxn];
bool vis[maxn];
unordered_map<ll, ll> mp_phi_d; //可换成unordered_map,约快3倍
ll S_phi_d(ll x) {
if (x < maxn) return sum_phi_d[x];
if (mp_phi_d[x]) return mp_phi_d[x];
ll ret = x * (x+) % mod * (*x%mod+) % mod * inv6 % mod; //%敲成*，浪费一个小时
for (ll i = , j; i <= x; i = j + ) {
j = x / (x / i);
ret =(ret - S_phi_d(x / i) * (i+j) % mod * (j-i+) % mod * inv2 % mod + mod) % mod;
}
return mp_phi_d[x] = ret;
}
void initPhi_d()
{
phi[] = ;
for (int i = ; i < maxn; i++) {
if (!vis[i]) pri[++tot] = i, phi[i] = i-;
for (int j = ; j <= tot && i * pri[j] < maxn; j++) {
vis[i * pri[j]] = true;
if (i % pri[j] == )
{
phi[i * pri[j]] = phi[i] * pri[j] % mod;
break;
}
else
{
phi[i * pri[j]] = phi[i] * phi[pri[j]] % mod;
}
}
}
for (int i = ; i < maxn; i++) sum_phi_d[i] = (sum_phi_d[i - ] + phi[i]*i) % mod;
}
ll solve(ll n)
{
ll res = ;
for(ll i = , j;i <= n;i = j+)
{
j = n / (n / i);
res = (res + S_phi_d(n/i) * (j-i+) % mod) % mod;
}
return (res+n)*inv2%mod;
} int main() {
initPhi_d();
scanf("%lld%lld", &a, &b);
printf("%lld\n", (solve(b)-solve(a-)+mod) % mod);
return ;
}
参考链接：