1. 程式人生 > 實用技巧 >bzoj4916 - 神犇和蒟蒻(杜教篩)

bzoj4916 - 神犇和蒟蒻(杜教篩)

題目


\(A=\sum\limits^N_{i=1}{\mu(i^2)}\)
\(B=\sum\limits^N_{i=1}{\varphi(i^2)}\)

題解

易知,\(A=1\)

由於\(\varphi(i^2)=i\varphi(i)\)\(B=\sum\limits^N_{i=1}{i\varphi(i)}\)

\(f(i)=i\varphi(i)\),發現\((f \ast id)(n)=\sum\limits_{d|n}{d\varphi(d)\frac{n}{d}}=n\sum\limits_{d|n}\varphi(d)=n^2\)

由杜教篩,最終得到

\(B=S(n)=\frac{n(n+1)(2n+1)}{6}-\sum\limits_{i=2}^{n}{iS(\lfloor\frac{n}{i}\rfloor)}\)

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <unordered_map>
using namespace std;
const int N = 2000010;
typedef long long ll;
const int M = 1e9 + 7;
const int r2 = 500000004;
const int r6 = 166666668;
ll T, n;
unordered_map<ll, ll> mp_phi2;
ll sum_phi2[N];
int phi[N], prime[N], cnt;
bool vis[N];


ll S_phi2(ll x) {
  if (x == 0) return 0;
  if (x < N) return sum_phi2[x];
  if (mp_phi2[x]) return mp_phi2[x];
  ll ret = x * (x + 1) % M * (2 * x + 1) % M * r6 % M;
  for (ll i = 2, j; i <= x; i = j + 1) {
    j = x / (x / i);
    ret -= S_phi2(x / i) * (i + j) % M * (j - i + 1) % M * r2 % M;
    ret %= M;
  }
  return mp_phi2[x] = (ret + M) % M;
}



int main() {
  int n;
  cin >> n;
  cout << 1 << endl;
  phi[1] = 1;
  for(int i = 2; i < N; i++) {
    if(!vis[i]) {
      prime[cnt++] = i;
      phi[i] = i - 1;
    }
    for(int j = 0; j < cnt; j++) {
      int p = prime[j];
      if(1ll * p * i > N) break;
      vis[p * i] = 1;
      if(i % p == 0) {
        phi[i * p] = p * phi[i];
        break;
      }
      phi[i * p] = phi[p] * phi[i];
    }
  }
  for(int i = 1; i < N; i++) {
    sum_phi2[i] = sum_phi2[i - 1] + 1ll * i * phi[i] % M;
    sum_phi2[i] %= M;
  } 
  cout <<  S_phi2(n) << endl;
}