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

链接：

https://loj.ac/problem/6229

题意：

$$F(n)=sum_{i=1}^nsum_{j=1}^ifrac{mathrm{lcm}(i,j)}{mathrm{gcd}(i,j)}$$

让你求 (F(n) mod1000000007)。

题解：

设(egin{align} f(n)=sum_{i=1}^nfrac{lcm(i,n)}{gcd(i,n)}&=sum_{i=1}^nfrac{n*i}{(i,n)^2}\ &=sum_{i=1}^nsum_{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}dsum_{i=1}^d[(i,d)=1]*i\ &=frac 12(1+sum_{d|n}d^2varphi(d)) end{align})。

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

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

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

所以，

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

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

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

整合一下，就是：

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

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

所以有：(ans={nover 2}+{1over 2}sum_{d=1}^nsum_{i=1}^{lfloor{nover d} floor}i^2varphi(i))。

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

其中，(sum_{i=1}^ni^2 = frac{ncdot(n+1)cdot(2n+1)}{6})，(varphi(i))为欧拉函数。

代码：

#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 << '
';
  cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.
";
  return 0;
}