「Luogu3768」简单的数学题-杜教筛

求
$\sum_{i = 1}^{n} \sum_{j = 1}^{n} i j gcd (i, j)) m o d p$
其中 $n < 10^{10}$ 。

题解🔗

$a n s = \sum_{i = 1}^{n} \sum_{j = 1}^{n} i j gcd (i, j)) m o d p = \sum_{d = 1}^{n} d \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} \sum_{j = 1}^{⌊ \frac{n}{d} ⌋} i j d^{2} [gcd (i, j) == 1] = \sum_{d = 1}^{n} d^{3} \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} \sum_{j = 1}^{⌊ \frac{n}{d} ⌋} i j [gcd (i, j) == 1]$

如果我们令

$f (d) = \sum_{i = 1}^{n} \sum_{j = 1}^{n} i j [gcd (i, j) == d] g (d) = \sum_{d | k} f (k) = d^{2} \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} \sum_{j = 1}^{⌊ \frac{n}{d} ⌋} i j g (d) = d^{2} {[\frac{⌊ \frac{n}{d} ⌋ (⌊ \frac{n}{d} ⌋ + 1)}{2}]}^{2}$

令 $s u m (x) = \frac{x (x + 1)}{2}$ ，原式化为：

$g (d) = d^{2} \cdot s u m (⌊ \frac{n}{d} ⌋)^{2}$

就有：

$f (d) = \sum_{d | k} μ (k) g (\frac{k}{d}) f (1) = \sum_{i = 1}^{n} μ (i) g (i) f (1) = \sum_{i = 1}^{n} μ (i) i^{2} s u m (⌊ \frac{n}{i} ⌋)^{2}$

那么：

$a n s = \sum_{d = 1}^{n} d^{3} \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} μ (i) i^{2} s u m (⌊ \frac{n}{d i} ⌋)^{2}$

枚举 $i d = T$ ，则有

$a n s = \sum_{T = 1}^{n} s u m (⌊ \frac{n}{T} ⌋)^{2} \sum_{d | T} d^{3} μ (\frac{T}{d}) \times {(\frac{T}{d})}^{2} = \sum_{T = 1}^{n} T^{2} s u m (⌊ \frac{n}{T} ⌋)^{2} \sum_{d | T} d μ (\frac{T}{d})$

有 $i d * μ = φ$ ，所以

$a n s = \sum_{T = 1}^{n} s u m (⌊ \frac{n}{T} ⌋)^{2} T^{2} φ (T)$

令 $f (x) = x^{2} φ (x)$ ，我们就有
$a n s = \sum_{T = 1}^{n} s u m (⌊ \frac{n}{T} ⌋)^{2} f (T)$

注意到 $⌊ \frac{n}{T} ⌋$ 只有根号个取值，所以我们想要处理出 $f (T)$ 的前缀和。

杜教筛：

$S (n) = \sum_{i = 1}^{n} h (i) - \sum_{d = 2}^{n} g (d) S (⌊ \frac{n}{d} ⌋)$

如果我们令 $g (n) = n^{2}$ ，那么

$h (i) = (g * f) (i) = \sum_{d | i} f (d) g (\frac{i}{d}) = \sum_{d | i} d^{2} φ (d) {(\frac{i}{d})}^{2} = \sum_{d | i} φ (d) i^{2} = i^{3}$

又因为
$\sum_{i = 1}^{n} i^{3} = {[\frac{n (n + 1)}{2}]}^{2}$

所以我们就有

$S (n) = \sum_{i = 1}^{n} h (i) - \sum_{d = 2}^{n} g (d) S (⌊ \frac{n}{d} ⌋) = {[\frac{n (n + 1)}{2}]}^{2} - \sum_{d = 2}^{n} d^{2} S (⌊ \frac{n}{d} ⌋)$

综上：

$a n s = \sum_{T = 1}^{n} s u m (⌊ \frac{n}{T} ⌋)^{2} f (T) S (n) = {[\frac{n (n + 1)}{2}]}^{2} - \sum_{d = 2}^{n} d^{2} S (⌊ \frac{n}{d} ⌋)$

代码实现一定要多取模…

代码🔗

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#include <bits/stdc++.h>
#define ll long long
using namespace std;

const int MAXN = 8500000;

ll pow(ll x,ll k,ll p){
  ll ans = 1;
  for(ll i = k;i;i>>=1,x = (x*x)%p) if(i & 1) ans = (ans*x)%p;
  return ans;
}

ll inv(ll x,ll p){return pow(x,p-2,p);}

ll p,n,rev_6,rev_2;
map<ll,ll> s;
ll pre_s[MAXN];

void sieve(){
  static ll prime[MAXN],phi[MAXN],cnt = 0;
  static int vis[MAXN];
  phi[1] = 1;
  for(int i = 2;i<=MAXN-1;i++){
    if(vis[i] == 0){
      prime[++cnt] = i;
      phi[i] = i-1;
    }
    for(int j = 1;i * prime[j] <= MAXN-1 && j<=cnt;j++){
      vis[i*prime[j]] = 1;
      if(i % prime[j] == 0){
        phi[i*prime[j]] = phi[i] * (prime[j]);
        break;
      }
      else{
        phi[i*prime[j]] = phi[i] * (prime[j]-1);
      }
    }
  }
  for(int i = 1;i<=MAXN-1;i++){
    pre_s[i] = phi[i]%p * i % p * i % p;
    pre_s[i] += pre_s[i-1];
    pre_s[i] %= p;
  }
}

ll sum(ll x){x%=p;return x%p *(x+1)%p *rev_2%p;}
ll sums(ll x){x%=p;return x%p *(x+1)%p *(x+x+1)%p *rev_6%p;}

ll S(ll n){
  if(n < MAXN) return pre_s[n];
  if(s.count(n))return s[n];
  ll ans = sum(n)%p;
  ans = ans*ans%p;
  for(ll l = 2,r;l <= n;l = r+1){
    r = (n/(n/l));
    ans -= S(n/l)%p * (sums(r) - sums(l-1)+p)%p;
    ans = ans % p;
  }
  return s[n] = ans % p;
}

void init(){
  scanf("%lld %lld",&p,&n);
  rev_6 = inv(6,p),rev_2 = inv(2,p);
}

ll calc(ll n){
  ll ans = 0;
  for(ll l = 1,r;l <= n;l = r+1){
    r = (n/(n/l));
    ll tmp = (S(r) - S(l-1) + p) % p;
    ll summ = sum(n/l) %p;
    ans += summ%p * summ%p * tmp%p;
    ans %= p;
  }
  return ans;
}

int main(){
  init();
  sieve();
  printf("%lld\n",calc(n));
  return 0;
}

「Luogu3768」简单的数学题-杜教筛

题解🔗

代码🔗

See Also

Comments