基本算数

一个整数可以用物体的个数来表示。比如用8个点我们就能表示数字8。但这样的话对于大的数字我们就必须用非常多的点,造成了不方便。因此为了更方便地表示一个数,我们通常要选择一个进制。我们在日常生活中常用的是10进制,在计算机中常用的是2进制,等等。

数\(N\)在\(a\)进制下需要\(\log_a N\)位,在\(b\)进制下需要\(\log_b N\),二者的比值\(\log_b^a\)是常数。因此我们可以认为,对进制的选择是不影响算法的时间复杂度的。

加法

从我们在小学时做的列竖式计算中就可以看出,逐位计算并进位的复杂度是正比于数的位数的。设数的位数为\(n\),则复杂度是\(O(n)\)。同时,仅仅读入数据就要消耗\(O(n)\),因此这就是加法的最优复杂度了。

乘法

按照列竖式的方法,乘法的复杂度是\(O(n^2)\)。

我们也可以通过这样的递归方法来完成乘法:

lua数学运算 数学的运算方法_乘法逆元

递归次数正比于\(y\)的位数\(n\),每次完成加法的复杂度也是\(n\),总的复杂度依然是\(O(n^2)\)。

乘法本质上是多项式乘法,通过快速傅里叶变换是可以优化到\(O(n \log n)\)的。我们在讨论“分治”算法时还会来讨论这个问题。

除法

我们想算出\(x\)除以\(y\)的商\(q\)和余数\(r\)。考虑递归,假设已经算出了\(\lfloor x/2\rfloor\)除以\(y\)的余数\(q'\)和\(r'\)。则有\(\lfloor x/2 \rfloor=q' \cdot y + r'\)。两边同时乘以2,得到\(2\lfloor x/2 \rfloor=2q' \cdot y + 2r'\)。如果\(x\)是偶数,那么有\(x=2q' \cdot y + 2r'\),否则有\(x = 2q' \cdot y + 2r'+1\)。注意如果\(2r'\)或\(2r'+1\)超出了\(y\),要减一个\(y\)。这样就推出了\(q,r\)。

递归的次数为\(n\),减\(y\)的复杂度为\(n\),总复杂度是\(O(n^2)\)。

模运算

对于\(x\)除以\(N\)的商为\(q\),余数为\(r\),可以写出\(x = qN+r\),其中\(r\)被限制在\(\{0,\cdots,N-1\}\)内。如果对于\(y\)也有\(y=q'N+r\),那么就称\(x\)和\(y\)是同余的。记作

\(x \equiv y \pmod N\)

由此可见\(x-y\)一定是\(N\)的倍数,即\(x-y \equiv 0 \pmod N\)。

如果\(x \equiv x' \pmod N\),同时\(y \equiv y' \pmod N\)。那么此时有\(x-x' = sN\), \(y-y' = tN\)。那么有\(x+y=(s+t)N+x'+y'\),所以\(x+y \equiv x'+y' \pmod N\)。而\(xy = (x'+sN)(y'+tN) = x'y'+(tx'+sy'+stN)N\),所以\(xy \equiv x'y' \pmod N\)。这告诉我们,在模运算中加法和乘法是可以被同余的项代换的。

模运算中由于每个数都被限制在\(\{0,\cdots,N-1\}\)内,因此加法不会超过\(2N-2\),因此加法以后最多只需要减一次\(N\)就能保证模运算的范围。两数相乘不会超过\((N-1)^2\),因此需要做一次除法(区域)才能保证范围,但我们注意到数字的位数最多只会扩展一倍。

注意,在模运算下两个都不为0的数相乘可能得到0。比如5*5在模25意义下为0。

快速幂

在模意义下计算\(x^y\)要用到快速幂算法。核心是:

lua数学运算 数学的运算方法_lua数学运算_02

递归的次数是\(n\),每次乘法用去\(n^2\),复杂度为\(O(n^3)\).

欧几里得算法(GCD)

对于\(x \geq y\),有\(gcd(x,y)=gcd(x\mod y,y)\)。一直这样递归下去,直到\(x\)成为\(y\)的倍数为止,我们就求出了最大公约数。这称为求解最大公约数的辗转相除法。要证明这个关系,只需证明辗转相减法\(gcd(x,y)=gcd(x-y,y)\),因为模运算只是减法运算的一种“捷径”。这是自然的,因为\(x,y\)都可以看作是由它们的公约数这个基本积木拼成的,它们的差依然是由同样的基本积木组成的。严格地,设\(x,y\)有一个公约数为\(k\),那么有\(x = s \cdot k\),\(y = t \cdot k\)。于是\(x-y = (s-t) \cdot k\)。因此\(x,y\)的任意一个公约数一定也是\(x-y,y\)的一个公约数。所以一定有\(gcd(x,y) \leq gcd(x-y,y)\)。反之同理,\(x-y,y\)的任意一个公约数一定也是\(x,y\)的一个公约数,所以又有\(gcd(x-y,y) \leq gcd(x,y)\)。联立即得\(gcd(x,y)=gcd(x-y,y)\)。这样我们就证明了辗转相除法。

为了分析辗转相除法的复杂度,我们来估计一下\(x \mod y\)的大小。如果\(y \leq x/2\),那么\(x \mod y < y \leq x/2\);如果\(y>x/2\),那么\(x \mod y = x-y < x/2\)。因此我们每一次至少把一个数缩小了一半,递归的次数是\(O(n)\)的。而模运算的计算复杂度就是除法的复杂度\(O(n^2)\),因此总的复杂度是\(O(n^3)\)。

扩展欧几里得算法(EXGCD)

对于方程\(ax+by=d\),左边可以提取出\(a,b\)的最大公约数,因此\(d\)一定是\(gcd(a,b)\)的倍数。换句话说,\(ax+by\)构成的集合一定是“量子化”的,\(gcd(a,b)\)是最小的单位。如果\(d\)不是\(gcd(a,b)\)的倍数,那么方程一定是无解的。

而我们要说明,\(ax+by=gcd(a,b)\)一定是有解的。由欧几里得算法得,\(gcd(a,b)=gcd(b,a \mod b)\)。我们列一个新的方程\(bx'+(a\mod b)y'=gcd(b,a \mod b)\)。联立可得\(ax+by=bx'+(a\mod b)y'\)。注意到\(a \mod b = a-b\lfloor a/b \rfloor\)。因此有\(ax+by=ay'+b(x'-\lfloor a/b \rfloor y')\)。假设我们已经递归地解得了\(x',y'\),就可知\(x=y',y=x'-\lfloor a/b\rfloor y'\)是一组可行的解。而由于\(gcd(a,b)\)递归的尽头\(b=0\),此时\(x\)一定有解,因此返回去最初的\(x,y\)一定有解。这样求解方程的方法就成为扩展欧几里得算法。

我们通过Exgcd已经求出了\((x_0,y_0)\)是\(ax+by=gcd(a,b)\)的一组解,如何用它去找到方程的所有可行解呢?如果\(a(x_0+p)+b(y_0+q)=gcd(a,b)\)也是一组解,那么必须满足\(ap+bq=0\)(很像线性代数中的零空间的解)。反过来,所有满足\(ap+bq=0\)的\(p,q\)一定能够保证\((x+p,y+q)\)还是方程的解。因此我们只需解出\(ap+bq=0\)的所有可行解\((p,q)\)即可。设\(a,b\)除掉彼此的最大公约数得到\(a',b'\),\(ap+bq=0\)当且仅当\(a'p+b'q=0\)。要使得\(a'p\)是\(b'\)的倍数,\(a'p\)必须是\(lcm(a',b')\)的倍数,而\(gcd(a',b')=1\),所以\(lcm(a',b')=a'b'\),\(a'p\)必须是\(a'b'\)的倍数,即\(p\)必须是\(b'\)的倍数,因此必须满足\(p=\dfrac{kb}{gcd(a,b)}\)。而容易验证所有这样的\(p\)都是满足的。综上,\(ax+by=gcd(a,b)\)的全部解就是\((x_0+\dfrac{kb}{gcd(a,b)},y_0-\dfrac{ka}{gcd(a,b)})\)。

对于一般的\(ax+by=d\),我们可以解出\(ax+by=gcd(a,b)\),再把得到的解乘上\(\dfrac{d}{gcd(a,b)}\)倍,得到一组可行的解\((x_0,y_0)\)。而后,完全相同的,我们依然利用\(ap+bq=0\)把特解拓展到通解,这个过程中只用到了关于\(a,b\)的性质,与等式右侧是无关的,因为我们所作的修正是完全相同的。得到的通解依然是\((x_0+\dfrac{kb}{gcd(a,b)},y_0-\dfrac{ka}{gcd(a,b)})\)。

乘法逆元

我们想要定义模运算下的除法,本质上就是定义模运算下的倒数。如果\(ax \equiv 1 \pmod N\),就称\(x\)是\(a\)在模\(N\)意义下的乘法逆元(其中\(a,x\)都在\(0\)到\(N-1\)之间)。

并不是所有\(a\)都有乘法逆元的。比如\(a=0\)就找不到这样的\(x\);如果\(a=2,N=6\),也可以验证没有任何的\(x\)是满足的。由扩展欧几里得,我们发现乘法逆元存在当且仅当\(gcd(a,N)=1\):如果\(ax \equiv 1\pmod N\)有解,则\(ax+Ny=1\)有解,这要求\(gcd(a,N)=1\)。而只要\(gcd(a,N)=1\)满足,由EXGCD就可以给出\(ax+Ny=1\)的一组解。这样得到的解可能会使得\(x\)没有落在\(0\)到\(N-1\)内,但根据通解的表达式,\(x\)只能间隔\(\dfrac{N}{gcd(a,N)}=N\)地变化。所以我们把\(x\)加上若干个\(N\)使它落在\(0\)到\(N-1\)就可以了,这一定是能够做到的。

对于给定的\(a\),逆元\(x\)是唯一的。因为如果存在\(x' \neq x\)使得\(x'\)也是\(a\)的乘法逆元,那么就有\(ax \equiv ax' \equiv 1 \pmod N\)。于是\(a(x-x')=tN\)。由于\(gcd(a,N)=1\),所以\(a(x-x')\)一定是\(a,N\)的最小公倍数的倍数,即\(a(x-x')=kaN\),于是\(x-x'=kN\),这与\(0 \leq x,x'<N\)矛盾。

素数

费马小定理

费马小定理指出:如果\(p\)是素数,那么对任意满足\(1 \leq a < p\)的整数\(a\),一定有\(a^{p-1} \equiv 1 \pmod p\)

我们发现, 对于\(i=1,\cdots,p-1\),\(a \cdot i\)一定是两两不同且不为0的(刚好构成了一个permutation):首先,如果\(a \cdot i \equiv 0\pmod p\),根据\(gcd(a,p)=1\),\(a\)存在乘法逆元,两边同时乘上\(a\)的乘法逆元,就会得到\(i \equiv 0\),矛盾;如果对于\(i \neq j\)成立\(a \cdot i \equiv a \cdot j\),那么同样乘以乘法逆元会得到\(i \equiv j\),矛盾。因此,\(a \cdot i\)刚好取遍了\(1,\cdots,p-1\)。把他们全都乘起来,就得到\(a^{p-1}!(p-1)! \equiv (p-1)! \pmod p\)。而\(1\)到\(p-1\)的每个数都没有\(p\)这个质因子,因此\((p-1)!\)一定是与\(p\)互质的,因此也有乘法逆元。于是两边同时乘以这个逆元,得到\(a^{p-1} \equiv 1 \pmod p\)。这样就证明了费马小定理。

遗憾的是,费马小定理的逆命题并不成立——我们不能由\(a^{N-1} \equiv 1 \pmod N\)推出\(N\)是素数。反例:合数\(341=11\times 31\),却有\(2^{340} \equiv 1 \pmod {341}\)。甚至存在这样的合数\(N\),使得所有比\(N\)小的与它互质的数\(a\)全都满足\(a^{N-1} \equiv 1 \pmod N\)。(\(N=561=3\times 11 \times 17\))这样的数被称为Carmichael数。

\(a^{N-1} \equiv 1 \pmod N\)虽然无法推出\(N\)是素数,但可以推出\(gcd(a,N)=1\)。如果\(gcd(a,N)>1\),那么\(a^{N-1}-tN\)一定是\(gcd(a,N)\)的倍数,因此永远不可能等于1。

根据费马小定理,我们又得到了另一种求乘法逆元的方式:\(p\)是素数时,\(a\)在模\(p\)下的逆就是\(a^{p-2} \mod p\)。

素数判断

我们可以证明,如果\(N\)不是Carmichael数,即存在\(a\)(与\(N\)互质)使得\(a^{N-1} \not\equiv 1 \pmod N\),那么这样的\(a\)在\(1\)到\(N-1\)间至少有一半。我们先任意固定这样的一个\(a\),它是一个常数,满足\(N-1\)方不为1。假如我们能找到一个\(b\)满足\(b^{N-1} \equiv 1 \pmod{N}\),那么对于\(ab\)(模意义下)有\((ab)^{N-1} \equiv a^{N-1}b^{N-1} \equiv a^{N-1} \not\equiv 1 \pmod {N}\)。因为\(a\)与\(N\)互质,所以有乘法逆元,所以对于某个\(b'\)满足\(b \not\equiv b' \pmod{N}\),一定成立\(ab \not\equiv ab' \pmod N\)。因此对于每个这样的\(b\),我们一定可以对应地找到一个\(ab\)。前者满足\(N-1\)次方等于1的条件,后者满足不等于\(N-1\)次方等于1的条件。由于\(b^{N-1} \equiv 1 \pmod{N}\)成立,\(b\)一定存在乘法逆元,而显然\(a \not \equiv 1 \pmod N\),所以\(b \not\equiv ab \pmod N\)一定成立。我们一对对选这样的\(b,ab\),会形成两个独立的集合,一个全都是\(N-1\)次方等于1的,一个全都是\(N-1\)次方不等于1的。也就意味着,每找到一个等于1的就一定能找到一个不等于1的,所以不等于1的至少占了一半。

所以我们可以设计这样一个素数判定算法,对任意给定的\(N\),我们在\(2\)到\(N-1\)间随机一个\(a\),计算\(a^{N-1} \mod N\)。如果答案为1就返回yes否则返回no。如果\(N\)是素数,那么程序一定正确;如果\(N\)不是素数且\(N\)不是Carmichael数,此时返回yes的概率小于\(1/2\)。也就是说此时程序出错的概率小于\(1/2\)。如果我们连着做\(k\)次,出错的概率就被降到了\(1/2^k\)。

由于Carmichael数是罕见的,我们直接忽略这种情况。这样我们就有了一个高效的算法,复杂度就是计算快速幂的复杂度\(O(n^3)\)。考虑到\(k\)后的复杂度为\(O(kn^3)\)。

素数生成

Lagrange指出\(N\)以内的素数个数同阶于\(\dfrac{N}{\ln N}\)。这说明素数是非常密集的。我们可以不停地随机一个\(N\)以内的数,然后用上面的算法判定它是否是一个素数,直到它通过判定为止。这样我们就大概率能得到一个素数了。

每个随机到的数是素数的概率为\(\dfrac{1}{\ln N}\),因此期望\(\ln N\)次随机后算法结束,共\(O(n)\)次。每次判定的复杂度是\(O(n^3)\)。因此总复杂度为\(O(n^4)\)。

密码学

一切发送出去的(通信、电波……)信息都是可能被监听的。想要达成加密的目的就必须有从未发送出去过的“密钥”。如果Alice和Bob都有一个二进制串\(r\)(密钥),那么Alice可以先把要发送的串与\(r\)异或,Bob收到后再与\(r\)异或一次就能得到明文了。但这要求\(r\)很长,而且两人要在不被窃听的情况下交换\(r\)是不容易的。

RSA

RSA是公钥加密的典范。在这里,Bob手上有两个密钥,一个密钥全世界只有他知道,称为私钥;一个密钥向全世界公开,称为公钥。Alice用公钥加密自己的信息,Bob利用私钥解密,这样就永远都不可能有人破解了。

在这里,Bob随机生成两个质数\(p,q\)。定义\(N=pq\),再任意找一个与\((p-1)(q-1)\)互质的数\(e\)(这是容易办到的,比如取\(e=3\))。\(N,e\)构成了我们的公钥。求出\(e\)在模\((p-1)(q-1)\)意义下的乘法逆元\(d\),把\(d\)作为我们的私钥。

由于\(ed \equiv 1 \pmod{(p-1)(q-1)}\),可以写出\(ed = k(p-1)(q-1)+1\)。而\(x^{p-1} \equiv 1\pmod p\),因此\(x^{ed}\equiv x^{k(p-1)(q-1)+1} \equiv x \cdot (x^{p-1})^{k(q-1)} \equiv x \pmod p\)。同理,也有\(x^{ed} \equiv x \pmod q\)。因此\(x^{ed}-x=t_1p=t_2q\)。说明\(x^{ed}-x\)是\(p,q\)的公倍数,因此有\(x^{ed} \equiv x \pmod{pq}\),即\((x^{e})^d \equiv x \pmod N\)。这个式子对于任意的\(0\)到\(N-1\)的\(x\)都成立,如果把\(f(x)=x^e\)看作一个映射,\(g(x)=x^d\)看作一个映射,那么\(g(f(x))=x\)。这说明\(g(x)\)的象集至少有\(N\)个元素,那么其定义域也至少要有\(N\)个元素。而\(x\)最多只有\(N\)个元素。因此\(f,g\)都必须是双射。所以对于任意的\(x\),我们都可以用\(f\)加密(它依赖于\(e,N\),因此是公钥加密),接收者都可以用\(g\)来解密(由于涉及到\(d\),因此是私钥解密)。

RSA的核心在于,作为公钥外人只知道\(N\)而不知道\(p,q\)。在\(N\)很大时是难以找出\(p,q\)的。因此也就难以找出\((p-1)(q-1)\),也就难以求出逆元\(d\)。我们把要发送的信息以\(x^e\)的方式加密,不知道\(d\)就无法解密,只有Bob能够解密。