初等数论与组合数学初步

引言

数论部分是 ACM 竞赛中最难的部分之一,本文仅作入门学习参考之用,文中省略了某些复杂定理的证明,详细证明过程请参考其他资料,文中的扩展欧几里得算法和中国剩余定理并不常用,供学有余力的同学学习

1. 同余

1.1 带余除法

1.1.1 定义 带余除法 设 a, b 是整数,且 b>0 ,则存在非负整数 q, r ,使得

$$a = bq + r$$

且 $0≤ r < b$, 称 q 为商, r 为余数

显然带余除法中的商和余数都是唯一的,在下文中将商记为 $a/b$ ,将余数记为 $a%b$,“/” 与 “%”的运算优先级与乘除法相同,当然,在 C 语言中二者分别对应 $a/b$ 与 $a%b$

注意到被除数 a 是可以为负数的,但如无特殊说明,下文不考虑负数

1.2 同余方程

1.2.1 定义 同余方程 设 a,b 是整数,若它们除以正整数 m 的余数相同,则称 a,b 对于模 m 同余,记为

$$a ≡ b (mod m)$$

1.3 “+”、“-”、“×”的取余运算

1.3.1 性质 加法的取余运算

$$(a+b)%c=(a%c+b%c)%c$$

$$(a+b)%c=((a/c)c+a%c+(b/c)c+b%c)%c=(a%c+b%c)%c$$

用类似的方法可以证得

1.3.2 性质 减法的取余运算

$$(a−b)%c=(a%c−b%c+c)%c$$

1.3.3 性质 乘法的取余运算

$$(ab)%c=(a%c)(b%c)%c$$

注意到 1.3.2 性质中等式右边加了 c ,是因为 a 可能是负数,为了便于写代码,将其处理成非负数

1.4 逆元

虽然取余运算对于“+”、“-”、“×”不难,但通常情况下

$$\frac{a}{b}%c ≠ \frac{a%c}{b%c}%c$$

如何计算 $\frac{a}{b}%c$ ?我们通常要找一个逆元 $b^{−1}$ ,使得 $bb^{−1} ≡ 1 ( mod c)$ ,那么就有

$$\frac{a}{b}%c = ab^{−1}%c$$

如果 c 是素数,根据下面的定理

1.4.1 费马小定理 设 b 是一个整数,c 是一个素数,那么

$$b^{c−1} ≡ 1 (mod c)$$

定理的证明需要引入简化剩余系的概念,略显复杂,此处省略

将上式改写一下,得到

$$bb^{c−2} ≡ 1 (mod c)$$

因此取 $b^{−1} = b^{c−2}$ 即可,一般需要用快速幂计算,但要注意,与除数不能为 0 类似,要保证 $b % c ≠ 0$

至于 c 不是素数的情况,有时逆元甚至不存在,如果存在,可用扩展欧几里得算法计算,将在下面 3 中介绍,如果不存在,可以考虑不求逆元,直接使用高精度计算 $\frac{a}{b}$

2. 最大公因数与最小公倍数

2.1 最大公因数

顾名思义,最大公因数就是公因数中最大的那个,我们记 a, b 的最大公因数为 $gcd(a,b)$,有如下性质

  • 2.1.1 性质 gcd(a,b)=gcd(b,a)
  • 2.1.2 性质 gcd(a,b)=gcd(a−b,b)(a≥b)
  • 2.1.3 性质 gcd(a,b)=gcd(a%b,b)
  • 2.1.4 性质 gcd(a,b,c)=gcd(gcd(a,b),c)

2.1.1 性质是显然的,2.1.2 性质是辗转相减法的原理,2.1.3 性质可以视为2.1.2 性质的“一步到位”版本,2.1.4 性质指出多个数的最大公因数可以递推地进行求解

2.2 辗转相除法

根据2.1.3 性质,得到辗转相除法的参考代码模板

typedef long long ll;
ll gcd(ll a, ll b){
    return b?gcd(b,a%b):a;
}

注意当 $b ≠ 0$ 时,返回值为 $gcd(b,a%b)$ 而不是 $gcd(a%b,b)$ ,否则会不断递归导致栈溢出

2.3 最小公倍数

顾名思义,最小公倍数就是公倍数中最小的那个,我们记 a, b 的最小公倍数为 $lcm(a,b)$,有如下性质

2.3.1 性质

$$lcm(a,b) = \frac{ab}{gcd(a,b)}$$

下面是最小公倍数的参考代码模板

ll lcm(ll a, ll b){
    return a/gcd(a,b)*b;
}

注意是先除后乘,避免在中间过程中数据超出64位整数的范围

3. 扩展欧几里得算法

3.1 扩展欧几里得算法

考虑方程

$$ax + by = c$$

其中 $a, b, c$ 是已知的正整数,如何求出方程的解呢?

3.1.1 定理 上述方程有解的充要条件是 $gcd(a,b)|c$ (c 是 gcd(a,b) 的倍数)

证明过程暂时省略

我们先将问题简化,求方程

$$ax + by = d$$

的解,其中 a,b 是正整数, $d = gcd(a, b)$

前文中所述的辗转相除法又名欧几里得算法,扩展欧几里得算法就是在此基础上实现的,下面举例来描述算法的过程

求方程 $36x + 28y = 4$ 的解

  • 36 = 28 × 1 + 8
  • 28 = 8 × 3 + 4
  • 8 = 4 × 2

将前两步改写一下:

  • 8 = 36 − 28 × 1
  • 4 = 28 − 8 × 3

再代入一下,得

4 = 28 − (36 − 28 × 1) × 3 = 28 × 4 + 36 × (−3)

因此 $x=−3,y=4$ 是方程的一组解,方程的任一解可以表示成

image_1dkg7t576jve1627d7d1ddqhn91t.png-5.5kB

更一般地,方程 $ax+by=d,d=gcd(a,b)$ 的所有解为

image_1dkg7tk7bqu2v4e1f8m1b16142p2a.png-7.9kB

其中 $x_{0}, y_{0}$ 是一组特解

回到最开始的问题,方程 $ax+by=c,gcd(a,b)|c$ 的所有解为

image_1dkg7u3391qn91rlf1s1v1v9911h32n.png-9kB

其中 $x_{0}, y_{0}$ 是方程 $ax + by = d, d = gcd(a, b)$ 的一组特解

下面是参考代码模板:

ll ext_gcd(ll a,ll b,ll& x,ll& y){
    ll d = a;
    if (!b){
        x = 1;y = 0;
    }else{
        d = ext_gcd(b,a%b,y,x);
        y -= a/b*x;
    }
    return d;
}

3.2 扩展欧几里得算法与逆元

回到1.4中挖的坑,要找到 $b^{−1}$ 使得 $bb^{−1} ≡ 1( mod c)$ ,实质上是求解方程 $bx + cy = 1$ 中的 $x$ ,当然只有 $gcd(b,c) = 1$ 时才有解,否则逆元不存在

4. 同余方程组与中国剩余定理

有神就有光,有方程就有方程组

image_1dkg7loht12dqiqd13nq1o1222vm.png-13.8kB

其中 $a_{1}, a_{2}, …, a_{n}$ 是整数, $m_{1}, m_{2}, …, m_{n}$ 是正整数

4.1.1 中国剩余定理(孙子定理) 设上述方程组中 $m_{1}, m_{2}, …, m_{n}$ 两两互质,则方程组的通解为

$$x = k \prod_{i=1}^{n}m_{i} + \sum_{i=1}^{n}a_{i}M_{i}M_{i}^{-1}$$

其中 $M_{i} = \prod_{j≠i} m_{j}$

定理的证明过程暂时省略

下面是参考代码模板,需要调用前面的扩展欧几里得算法模板

ll Sunzi(ll *m,ll *a,int len){
    ll lcm = 1;
    ll ans = 0;
    for (int i=0;i<len;i++){
        ll k0,ki;
        ll d = ext_gcd(lcm,m[i],k0,ki);
        if ((a[i]-ans)%d!=0) return -1;
        else {
            ll t = m[i]/d;
            k0 = ( k0*(a[i]-ans)/d%t + t)%t;
            ans = k0*lcm + ans;
            lcm = lcm/d*m[i];
        }
    }
    return ans;
}

5. 素数

素数是只有 1 和它本身两个因数的数,1 不是素数

5.1 素数的判断

bool isPrime(ll n){
    if(n==1)return false;
    for(ll i=2;i*i<=n;i++)
        if(n%i==0)return false;
    return true;
}

这是最简单的素数的判断的参考代码模板,复杂度为 $O(\sqrt{n})$

原理其实很简单,对于一个大于 1 的整数,如果 x 是它的一个大于 $\sqrt{n}$ 的因子,那么 $\frac{n}{x}$ 是它的小于 $\sqrt{n}$ 的因子

在大多数情况下,这种判断方式的复杂度已经足够小了,如果要吹毛求疵追求更低的复杂度,可以使用下面的方法

bool isPrime(ll n){
    if(n==2||n==3||n==5)return 1;
    if(n%2==0||n%3==0||n%5==0||n==1)return 0;
    ll c=7,a[8]={4,2,4,2,4,6,2,6};
    while(c*c<=n)for(auto i:a){if(n%c==0)return 0;c+=i;}
    return 1;
}

这种方法至少能保证int范围内的数是正确的,当 n 过大时可能会出错,笔者暂时还不知道最小的会出错的数是多少

5.2 素数筛

如果要求出不超过 n 的所有素数,素数筛是最好的选择,下面是一种朴素的筛法

void getPrime(bool p[],int n){
    for(int i=1;i<=n;i++)p[i]=true;
    p[1]=false;
    for(int i=2;i<=n;i++){
        if(p[i]){
            for(int j=i+i;j<=n;j+=i)p[j]=false;
        }
    }
}

这种方法的原理是从小到大将素数的倍数筛掉,复杂度为 $O(nlogn)$ ,注意到每个合数如果有多个素因子,那么就会被重复筛掉,造成复杂度的浪费,因此,用下面的方法可以保证每个合数只被它最小的素因子筛掉一遍,以 $O(n)$ 的复杂度解决上述问题

ll getPrime(ll n,bool vis[],ll prime[]){
    ll tot=0;
    for(ll i=1;i<=n;i++)vis[i]=0;
    for(ll i=2;i<=n;i++){
        if(!vis[i])prime[tot++]=i;
        for(ll j=0;j<tot;j++){
            if(prime[j]*i>n)break;
            vis[prime[j]*i]=1;
            if(i%prime[j]==0)break;
        }
    }
    return tot;
}

上述代码可能略显难懂,笔者也不知如何讲清原理,请读者仔细体会

6. 组合数学初步

6.1 组合数

4.1.1 定义 组合数 在 $n(n≥0)$ 个不同元素中选取 $m(0≤m≤n)$ 个元素,不同的取法记为

$$C_{n}^{m} = \frac{n!}{m!(n-m)!}$$

6.2 杨辉三角

组合数与杨辉三角中的数字是一一对应的

image_1dkg7ml6h10241hn4dte1vf99et13.png-11.5kB
杨辉三角的自然数形式

image_1dkg7narn18n41id71gnvbqd0i1g.png-23.5kB

杨辉三角的组合数形式

按照上面的写法,杨辉三角的第 $n$ 行第 $m$ 列即为 $C_{n+m−2}^{m−1}$

注意到上图中每个数等于其左边的数与上边的数(如果有的话)之和,这就是杨辉恒等式

$$C_{n}^{m} = C_{n−1}^{m-1} + C_{n−1}^{m}$$

6.3 组合数的取余计算

在 ACM 竞赛中,我们常常需要计算 $C_{n}^{m}%p$ ,可以参考下面两种方法

如果 $n, m$ 不大,可以开 $O(n^{2})$ 的空间,可以利用杨辉恒等式来预处理组合数表

const ll mo=1e9+7;
ll C[1005][1005];
void getC(int n){
    for(int i=0;i<=n;i++){
        for(int j=0;j<=i;j++){
            if(j==0 || j==i)
                C[i][j]=1;
            else
                C[i][j]=(C[i-1][j-1]+C[i-1][j])%mo;
        }
    }
}

如果 $n, m$ 比较大,可以开 $O(n)$ 的空间,可以利用前文所述的逆元来求解,当然,要保证 $p$ 是素数

const ll mo=1e9+7;
ll C(ll n,ll m){
    static ll M=0,inv[N],mul[N],invMul[N];
    while(M<=n){
        if(M){
            inv[M]=M==1?1:(mo-mo/M)*inv[mo%M]%mo;
            mul[M]=mul[M-1]*M%mo;
            invMul[M]=invMul[M-1]*inv[M]%mo;
        }
        else mul[M]=1,invMul[M]=1;
        M++;
    }
    return mul[n]*invMul[m]%mo*invMul[n-m]%mo;
}

上面的代码中用 $O(n)$ 的复杂度处理了 $[1,n]$ 的逆元(原理在此不深究,当然读者可以使用上文中提到的用快速幂求逆元的方法,在大多数情况下足够了),处理 $Q$ 次 $n, m ≤ N$ 的询问的总复杂度为 $O(N+Q)$

附录

题目列表

  • codeforces 1108B Divisors of Two Integers 想法题
  • codeforces 299A Ksusha and Array 最大公因数
  • codeforces 1033B Square Difference 素数
  • codeforces 822A I’m bored with life 最大公因数
  • codeforces 131C The World is a Theatre 组合数
  • codeforces 898B Proper Nutrition 扩展欧几里得

这些题目应该都不难

……

Last modification:September 11, 2019
如果觉得我的文章对你有用,请随意赞赏