[luogu]P3768 简单的数学题(莫比乌斯反演,杜教筛)

最新推荐文章于 2021-06-20 22:45:38 发布

Hardict

最新推荐文章于 2021-06-20 22:45:38 发布

阅读量342

点赞数

CC 4.0 BY-SA版权

分类专栏：数论-各种数论函数文章标签：莫比乌斯反演杜教筛

本文链接：https://blog.youkuaiyun.com/Hardict/article/details/82940301

数论-各种数论函数同时被 2 个专栏收录

9 篇文章

订阅专栏

math

8 篇文章

订阅专栏

题意

求 $\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;
}