LOJ 6229 LCM / GCD (杜教筛+Moebius)

链接：

https://loj.ac/problem/6229

题意：

\[F(n)=\sum_{i=1}^n\sum_{j=1}^i\frac{\mathrm{lcm}(i,j)}{\mathrm{gcd}(i,j)}\]

让你求 \(F(n) \bmod1000000007\)。

题解：

设\(\begin{align} f(n)=\sum_{i=1}^n\frac{lcm(i,n)}{gcd(i,n)}&=\sum_{i=1}^n\frac{n*i}{(i,n)^2}\\ &=\sum_{i=1}^n\sum_{d|n}[(i,n)=d]\frac{n*i}{d^2}\\ &=\sum_{d|n}\sum_{i=1}^{[\frac nd]}[(i,\frac nd)=1]\frac{n*i}d\\ &=\sum_{d|n}d\sum_{i=1}^d[(i,d)=1]*i\\ &=\frac 12(1+\sum_{d|n}d^2\varphi(d)) \end{align}\)。

即求 \(\sum_{i=1}^n\sum_{d|i}d^2\varphi(d)=\sum_{i=1}^n\sum_{d=1}^{[\frac ni]}d^2\varphi(d)\)。

令 \(\phi'(n)=\sum_{i=1}^ni^2\varphi(i)\)。

因为 \(\sum_{d|n}d^2\varphi(d)*(\frac nd)^2=n^2\sum_{d|n}\varphi(d)=n^3\)。

所以，

\(\begin{align} \sum_{i=1}^ni^3=[\frac{n(n+1)}{2}]^2&=\sum_{i=1}^n\sum_{d|i}d^2\varphi(d)*(\frac id)^2\\ &=\sum_{i=1}^ni^2\sum_{d=1}^{[ \frac ni]}d^2\varphi(d)\\ &=\sum_{i=1}^ni^2\phi'([\frac ni]) \end{align}\)。

所以得到：\(\phi'(n)=[\frac{n(n+1)}{2}]^2-\sum_{i=2}^ni^2\phi'([\frac ni])\)。

可以杜教筛先预处理前 \(n^{2/3}\)，原问题可以在复杂度\(O(n^{2/3}log(n))\)内解决。

整合一下，就是：

推公式可以得到( 结合公式4 )：\(ans=\sum_{d=1}^n\sum_{i=1}^{\lfloor{n\over d}\rfloor}\sum_{j=1}^i ij[\gcd(i,j)=1]\)。

因为存在恒等式：\(\sum_{i=1}^ni[\gcd(i,n)=1]={[n=1]+n\varphi(n)\over 2}\)。

所以有：\(ans={n\over 2}+{1\over 2}\sum_{d=1}^n\sum_{i=1}^{\lfloor{n\over d}\rfloor}i^2\varphi(i)\)。

考虑 \(\sum_{i=1}^{n}i^2\varphi(i)\)出现的次数，可以得到： \(ans={n\over 2}+{1\over 2}\sum_{i=1}^ni^2\varphi(i)\lfloor{n\over i}\rfloor\)。

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int maxn = 1e6+100;
const int mod = 1e9+7;
int n;
int p[maxn],phi[maxn],pre[maxn];

int inv2,inv6;
ll qpower(ll a,ll b,ll mod)
{
  ll res = 1;
  while(b>0) {
    if(b&1) res = res * a % mod;
    b >>= 1;
    a = a * a % mod;
  }
  return res;
}
void init(int n)
{
  phi[1]=1;
  for(int i=2;i<=n;i++)
  {
    if(p[i]==0) p[++*p]=i,phi[i]=i-1;
    for(int j=1;j<=*p && 1LL*p[j]*i<=n;j++)
    {
      p[p[j]*i]=1;
      if(i%p[j]) phi[i*p[j]]=phi[i]*phi[p[j]];
      else
      {
        phi[i*p[j]]=phi[i]*p[j];
        break;
      }
    }
  }
  for(int i=1;i<=n;i++) {
    pre[i]=(pre[i-1]+1LL*phi[i]*i%mod*i)%mod;
  }
}
map<ll,int> mp;
int calcinv2(ll l,ll r)
{
  l %= mod;
  r %= mod;
  return (r - l + 1) * (l + r) % mod * inv2 % mod;
}
int calcinv6(ll n)
{
  n %= mod;
  return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}
int calc2(ll l,ll r)
{
 return (calcinv6(r) - calcinv6(l-1) ) % mod;
}
int calc3(ll n)
{
  return 1LL * calcinv2(1 , n) * calcinv2(1 , n) % mod;
}
int S(ll n)
{
  if(n <= 1e6) return pre[n];
  if(mp.count(n)) return mp[n];
  int res = calc3(n);
  for(ll i = 2, j; i <= n ;i = j + 1) {
    j = n / (n / i);
    res = (res - 1LL * calc2(i,j) * S(n / i)) % mod;
  }
  return mp[n] = res;
}
int main(int argc, char const *argv[]) {

  ll n;
  std::cin >> n;
  init(1000000);// 2/3
  inv2 = qpower(2,mod-2,mod);
  inv6 = qpower(6,mod-2,mod);
  int ans = 0;
  int last = 0;
  for(ll i = 1, j; i <= n; i = j + 1) {
    j = n /( n / i );
    int cur = S(j);
    ans = (ans + 1LL * (cur - last) * ( n / i)) % mod;
    last  = cur;
  }
  ans = (ans + n) % mod * inv2 % mod;
  std::cout << (ans + mod) % mod << '\n';
  cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.\n";
  return 0;
}

原文地址：https://www.cnblogs.com/LzyRapx/p/8459266.html