整除分块

最后更新于：2023年2月24日 下午

ACMer 考虑算法时总会优先考虑时间复杂度，这里介绍几个优美的根号复杂度的算法（原来这个叫整除分块）

可能需要 这里的基础知识

\(s(n) = \sum_{i=1}^{n} \lfloor \frac{n}{i} \rfloor\)

由于 \(\lfloor \frac{n}{i} \rfloor\) 的取值个数不会超过 \(2\sqrt{n}\)，所以可能存在 \(O(\sqrt{n})\) 的算法：

1
2
3
4
5
6
7
8
LL getsum(LL n){ // The code is simple and easy to understand
    LL sum = 0;
    for(LL i=1,j;i<=n;i=j+1){
        j = n/(n/i);
        sum += (j-i+1)*(n/i);
    }
    return sum;
}

事实上， \(s(n)\) 表示图像 \(xy=1\) 下方正整点的个数

\(\sigma_k(n) = \sum_{d|n} d^k\)

  1. \(\sigma_0(n)\) 表示正因子个数
  2. \(\sigma_1(n)\) 表示正因子之和
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
LL mypow(LL x,LL n){
    LL r = 1;
    while(n){
        if(n&1) r=r*x;
        n>>=1;  x=x*x;
    }
    return r;
}
LL getr(LL n,LL k){
    LL r = 0,d;
    for(d=1;d*d<n;++d){
        if(n%d==0)  r += mypow(d,k) + mypow(n/d,k);
    }
    if(d*d == n) r+=mypow(d,k);
    return r;
}

\(f_k(n) = \sum_{i=1}^n \sigma_k(i)\)

\[ f_k(n) = \sum_{i=1}^n \sigma_k(n) = \sum_{i=1}^n \sum_{d|i} d^k = \sum_{d=1} d^k \sum_i ^n[d|i] = \sum_{d=1} ^n d^k \lfloor \frac{n}{d} \rfloor \]

如果我们已经得到 \(ts[n] = \sum_{i=1}^n i^k\) 类似问题一，我们有以下 C++ 代码：

1
2
3
4
5
6
7
8
LL getf(LL n){
    LL sum = 0;
    for(LL i=1,j;i<=n;i=j+1){
        j = n/(n/i);
        sum += (ts[j]-ts[i-1])*(n/i);
    }
    return sum;
}

事实上我们有： \[ ts[n] = \left \lbrace \begin{array}{ll} \frac{n(n+1)}{2} & k=1 \\ \frac{n(n+1)(2n+1)}{6} & k=2 \\ \frac{n^2(n+1)^2}{4} & k=3 \\ \end{array} \right. \]

更一般的我们有 \[ 1^p+2^p+ \dots + n^p = \sum _{k=1} ^p \; (\; \sum_{j=1} ^ {k} (-1)^{k-j} C_k^j j^p \;) \; C_{n+1} ^{k+1} \]

\(g(n) = \sum_{i=1}^n\psi(i)\)

这里 \(\psi(n)\) 是欧拉函数表示小于 \(n\) 且和 \(n\) 互素的数的个数。 Euler's 乘积公式: \[ \psi(n) = n \prod _{p|n}( 1-\frac{1}{p} ) \] 我们现在开始计算 \(g(n)\) \[ g(n) = \sum_{i=1}^n \psi(i) = \sum_{1 \leq x \leq y \leq n , \gcd(x,y)=1} 1 \] 我们定义 \[ g_k(n) = \sum_{1 \leq x \leq y \leq n , \gcd(x,y)=k} 1 = g(\lfloor \frac{n}{k} \rfloor) \] 并且 \(\sum_{i=1}^n g_k(i) = \sum_{1 \leq x \leq y \leq n} 1 = \frac{n(n+1)}{2}\) ，于是我们知道 \[ g(n) = \frac{n(n+1)}{2} - \sum_{k=2}^n g(\lfloor \frac{n}{k} \rfloor) \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
const int N=1000006;
LL ans[N];
void init(){
    memset(ans,-1,sizeof(ans));
    ans[1]=1;ans[2]=2;
}
LL getans(int n){
    if(n<N&&ans[n]!=-1) return ans[n];
    LL r = n*(n+1)/2;
    for(int i=2,j;i<=n;i=j+1){
        j = n/(n/i);
        r -= (j-i+1)*getans(n/i);
    }
    if(n<N) ans[n]=r;
    return r;
}

\(f(n) = \sum_{i=1}^n \sum_{j=1}^n [\gcd(i,j)==1]\)

\[ f(n) = \sum_{i=1}^n \sum_{j=1}^n [\gcd(i,j)==1] = \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{d|i,d|j} \mu(d) = \sum_{d=1}^n \mu(d) (\lfloor \frac{n}{d} \rfloor)^2 \]

实际上画图可知，\(f(n)\) 还有一个表达式：

\[ f(n)= 2 \sum_{i=1}^{n} \phi(i) - 1 \]

我们有不妨设 \(n \leq m\)

\(r(n,m)=\sum_{i=1}^n \sum_{j=1}^m f(\gcd(i,j))\)

\[ \begin{aligned} r(n,m) &=\sum_{i=1}^n \sum_{j=1}^m f(gcd(i,j)) \\ &= \sum_{d=1}^n f(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{d} \rfloor} \sum_{t|i,t|j} \mu(t) \\ &= \sum_{d=1}^n f(d) \sum_{t=1}^{\lfloor \frac{n}{d} \rfloor} \mu(t) \lfloor \frac{n}{dt} \rfloor \lfloor \frac{m}{dt} \rfloor \\ &= \sum_{G=1}^n \lfloor \frac{n}{G} \rfloor \lfloor \frac{m}{G} \rfloor (f \star \mu)(G) \end{aligned} \]

特别地，\(g(n, m) = \sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)==1] = \sum_{d=1}^n \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \mu(d)\)

Problem hdu5608

已知 \(n^2 -3n+2 = \sum_{d|n} f(d)\) 计算 \(h(n) = \sum_{i=1}^n f(i) \mod 10^9+7\)

由于 \[ \sum_{i=1}^n \sum_{d|i} f(d) = \sum_{i=1}^n f(i) \lfloor \frac{n}{i} \rfloor = \sum_{i=1}^n h(\lfloor \frac{n}{i} \rfloor) \] 我们知道 \[ \sum_{i=1}^n \sum_{d|i} f(d) = \sum_{i=1}^n i^2-3i+2 = \sum_{i=1}^n (i-1)(i-2) = \frac{n(n-1)(n-2)}{3} \] 所以 \[ h(n) = \frac{n(n-1)(n-2)}{3} - \sum_{i=2}^n h(\lfloor \frac{n}{i} \rfloor) \]

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
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1000006;
const int M = 1e9+7;
const int inv3 = (M+1)/3;
int ans[N];
map<int,int> mp;
void init(){
    for(int i=0;i<N;++i){
        ans[i] = LL(i-1)*(i-2)%M;
    }
    for(int i=1;i<N;++i){  // Pretreatment acceleration
        for(int j=i<<1;j<N;j+=i){
            ans[j] -= ans[i];
            if(ans[j] < 0) ans[j] += M;
        }
    }
    for(int i=2;i<N;++i){
        ans[i] += ans[i-1];
        if(ans[i] > M)  ans[i] -= M;
    }
}
int getans(int n){
    if(n<N) return ans[n];
    map<int,int>::iterator it = mp.find(n); //Memory search
    if (it != mp.end())  return it->second;
    int r = LL(n)*(n-1)%M*(n-2)%M*inv3%M;
    for(int i=2,j;i<=n;i=j+1){
        j = n/(n/i);
        r -= LL(j-i+1)*getans(n/i)%M;
        if(r<0) r+=M;
    }
    mp.insert(pair<int,int>(n,r));
    return r;
}
int main(){
    init();
    int T,n;
    cin>>T;
    while(T--){
        scanf("%d",&n);
        cout<<getans(n)<<endl;
    }
    return 0;
}