1. 程式人生 > >[luogu]P3768 簡單的數學題(莫比烏斯反演,杜教篩)

[luogu]P3768 簡單的數學題(莫比烏斯反演,杜教篩)

題意

\sum_{i=1}^N \sum_{j=1}^N ijgcd(i,j),答案膜素數p

資料範圍:n\leq 10^{10},p\leq 1.1\times 10^{9}

題解

f[d]=\sum \sum ij[gcd(i,j)==d],F[d]=\sum \sum ij[d|gcd(i,j)]

\\ ans=\sum_{i=1}^N df[d]=\sum_{i=1}^N d\sum_{d|x} \mu (\frac{x}{d}) F[x]\\=\sum_{d=1}^N d\sum_{i=1}^{\frac{N}{d}}\mu (i)F[di]=\sum_{j=1}^N F[j] \sum_{x|j} x\mu (\frac{x}{j})

如果所以尤拉函式和狄利克雷卷積,可以知道(\mu * id)(N)=\varphi (N)(我一開始也沒反應過來,可以設為(\mu *f)進一步推導)

\\ F[d]=\sum\sum ij[d|gcd(i,j)]=\sum_{i=1}^{\frac{N}{d}} (id) \sum_{j=1}^{\frac{N}{d}} (jd)=d^2 S_1(\frac{N}{d})\quad (S_1(N)=1+2+...+N)

ans=\sum_{j=1}^N S_1(\frac{N}{j}) j^2\varphi(j),S_1可以整除分塊求,問題就變為化解j^2\varphi(j)(特指杜教篩操作)

g(x)=x^2,g為積性函式,(j^2\varphi(j) * g)(N) = N^2 \sum_{d|N} \varphi(d)=N^3

由杜教篩

\sum_{i=1}^N i^3 =\sum_{i=1} \sum_{d|i} (d^2\varphi(d) * \frac{i^2}{d^2})=\sum_{i=1} g(x)S_2(\frac{N}{i})\quad (S_2(N)=\sum_{i=1}^N i^2\varphi(i))

那麼S_2(N)+\sum_{i=2}^N i^2S_2(\frac{N}{i})=\sum_{i=1}^N i^3

就可以快速求得字首積了(預處理範圍取(10^{10})^{\frac{2}{3}})

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <unordered_map>
using namespace std;

using LL = long long;
const int MAXN = 8e6 + 5;
LL N, P;
bool noprime[MAXN];
int prime[MAXN], cnt_p, phi[MAXN];
void Euler_Sieve(int top);
LL pre[MAXN];
unordered_map<LL, LL> sum;
LL inv_2, inv_6;
LL getpre(LL);
LL qpow(LL, LL);
LL getsum(LL, int);

int main(){
  cin >> P >> N;
  inv_2 = (P + 1) / 2, inv_6 = qpow(6, P - 2);
  Euler_Sieve(MAXN - 5);
  LL ans = 0;
  LL l, r;
  for(l = 1; l <= N; l = r + 1){
    r = N / (N / l);
    LL tmp = getpre(r) - getpre(l - 1);
    tmp = (tmp % P + P) % P;

    LL x = getsum(N / l, 3);
    tmp = tmp * x % P;

    ans = (ans + tmp) % P;
  }
  cout << ans << endl;
  return 0;
}

LL getpre(LL n){
  if(n <= MAXN - 5) return pre[n] % P;
  if(sum[n]) return sum[n];
  LL l, r, res;
  res = getsum(n, 3);
  for(l = 2; l <= n; l = r + 1){
    r = n / (n / l);
    LL tmp = (getsum(r, 2) - getsum(l - 1, 2) + P) % P;
    res -= getpre(n / l) * tmp % P; res %= P;
  }
  res = (res + P) % P;
  return sum[n] = res;
}

LL qpow(LL x, LL n){
  LL res = 1;
  while(n){
    if(n & 1) res = res * x % P;
    x = x * x % P;
    n >>= 1;
  }
  return res;
}

void Euler_Sieve(int top){
  int i, j;
  phi[1] = 1;
  for(i = 2; i <= top; i++){
    if(!noprime[i])
      prime[++cnt_p] = i, phi[i] = i - 1;
    for(j = 1; j <= cnt_p && prime[j] * i <= top; j++){
      noprime[prime[j] * i] = true;
      if(i % prime[j] == 0){
        phi[prime[j] * i] = phi[i] * prime[j];
        break;
      }
      phi[prime[j] * i] = phi[i] * (prime[j] - 1);
    }
  }

  for(j = 1; j <= top; j++)
    pre[j] = (pre[j - 1] + 1LL * phi[j] * j % P * j % P) % P;
}

LL getsum(LL n, int op){
  n %= P;
  if(op == 1)
    return n * (n + 1) % P * inv_2 % P;
  else if(op == 2)
    return n * (n + 1) % P * (2 * n + 1) % P * inv_6 % P;
  else
    return getsum(n, 1) * getsum(n, 1) % P;
}