[51nod1227]平均最小公倍数（莫比乌斯反演+杜教筛）

题意

求 $sum_{i=a}^b sum_{j=1}^i frac{lcm(i,j)}{i}$.

分析

只需要求出前缀和，

$$egin{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} ight floor} sum_{j=1}^i j cdot [gcd(i,j)=1]
end{aligned}$$

其后面部分提出来，即求 $sum_{i=1}^n icdot [gcd(i,n)=1]$，对于这种一个值固定的gcd求和有一个套路，即倒序两两配对：

若 $n=1$，和为1；

若 $n>1$，因为 $gcd(i, n) = gcd(n-i, n)$ 且 $displaystyle  sum_{i=1}^n icdot [gcd(i,n)=1] = sum_{i=1}^{n-1} icdot [gcd(i,n)=1]$，

$displaystyle sum_{i=1}^{n-1} i cdot [gcd(i,n)=1] + sum_{i=n-1}^1icdot [gcd(i,n)=1] = nvarphi (n)$，所以和为 $ nvarphi (n) /2$.

综合得 $displaystyle sum_{i=1}^n icdot [gcd(i,n)=1] = frac{nvarphi (n)+[n=1]}{2}$.

具体实现上，$[i=1]$ 只成立 $n$，除2可以提出来，即 原式 = $displaystyle frac{1}{2}(sum_{d=1}^n sum_{i=1}^{left lfloor frac{n}{d} ight floor} ivarphi (i) + n)$.

现在唯一得问题是如何求 $displaystyle  S(n) = sum_{i=1}^n ivarphi (i)$.

根据杜教筛，

设 $displaystyle S(n) = sum_{i=1}^n f(i)$，$f(n) = nvarphi (n), g(n) = n$.

$displaystyle f*g = sum_{d|n} d varphi (d) cdot frac{n}{d} = nsum_{d|n}varphi (d) = n^2$.

因此 $displaystyle S(n) = sum_{i=1}^n i^2 - sum_{d=2}^n dcdot S(left lfloor frac{n}{d} ight floor)$.

再最外层套个数论分块即可。

代码

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <map>
#include<unordered_map>
using namespace std;
const int maxn = 2000010;
typedef long long ll;
const ll mod = 1000000007;
const ll inv2 = (mod+1)>>1;
const ll inv6 = 166666668;  //166666668
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+1) % mod * (2*x%mod+1) % mod * inv6 % mod;  //%敲成*，浪费一个小时
  for (ll i = 2, j; i <= x; i = j + 1) {
    j = x / (x / i);
    ret =(ret - S_phi_d(x / i) * (i+j) % mod * (j-i+1) % mod * inv2 % mod + mod) % mod;
  }
  return mp_phi_d[x] = ret;
}
void initPhi_d()
{
  phi[1] = 1;
  for (int i = 2; i < maxn; i++) {
    if (!vis[i]) pri[++tot] = i, phi[i] = i-1;
    for (int j = 1; j <= tot && i * pri[j] < maxn; j++) {
      vis[i * pri[j]] = true;
      if (i % pri[j] == 0)
      {
          phi[i * pri[j]] = phi[i] * pri[j] % mod;
          break;
      }
      else
      {
          phi[i * pri[j]] = phi[i] * phi[pri[j]] % mod;
      }
    }
  }
  for (int i = 1; i < maxn; i++) sum_phi_d[i] = (sum_phi_d[i - 1] + phi[i]*i) % mod;
}
ll solve(ll n)
{
    ll res = 0;
    for(ll i = 1, j;i <= n;i = j+1)
    {
        j = n / (n / i);
        res = (res + S_phi_d(n/i) * (j-i+1) % mod) % mod;
    }
    return (res+n)*inv2%mod;
}


int main() {
  initPhi_d();
  scanf("%lld%lld", &a, &b);
  printf("%lld
", (solve(b)-solve(a-1)+mod) % mod);
  return 0;
}

参考链接：

1. https://blog.csdn.net/FromATP/article/details/74999989

2. https://www.cnblogs.com/owenyu/p/7397687.html