2.3 大数乘法

774 阅读7分钟

大数乘法一般是给一个线性表,string格式的,让我们来求。

算法题我一般不写注释,大家如果不明白,就抄下代码对着过程一步步来。实在不会可以评论,2023上半年我还是有时间的。希望二月考研结果不错吧。如果考上了我想在b站给大家拍视频去讲,文字太生硬了。

简介

所谓大数相乘(Multiplication algorithm),就是指数字比较大,相乘的结果超出了基本类型的表示范围,所以这样的数不能够直接做乘法运算。

主要是有以下几种方法:

  • 小学生乘法:模拟小学生运算,最简单的乘法竖式手算的累加型。
  • 分治乘法:最简单的是Karatsuba乘法,一般化以后有Toom-Cook乘法
  • 快速傅里叶变换FFT:时间复杂度O(N lgN lglgN)。具体可参照Schönhage–Strassen algorithm
  • 中国剩余定理:把每个数分解到一些互素的模上,然后每个同余方程对应乘起来就行。
  • Furer’s algorithm:在渐进意义上FNTT还快的算法。

做参考的题目:43. 字符串相乘

小学生乘法

对比下图的手算步骤,我们的算法可总结成两步:

  1. 逐位相乘aibja_i * b_j
  2. 对应相加。

image.png 需要注意的事项:

  1. 符号单独判断。
  2. 需要一个int数组存储各个位的结果。
  3. 两重for循环算出aibja_i * b_j,并且进行移位操作(观察发现aibja_i * b_j经移位后的位置为i+ji+j
  4. 最后转成字符串输出。

代码:

class Solution {
public:
    string multiply(string num1, string num2) {
        string ans;
        vector<int> vec(num1.size() + num2.size());
        for(int i = num1.size() - 1;i >= 0;i--)
        for(int j = num2.size() - 1;j >= 0;j--){
            int num = (num1[i] - '0')  * (num2[j] - '0');
            num += vec[i + j + 1];
            vec[i + j + 1] = num % 10;
            vec[i + j] += num / 10;
        }
        for(int i = 0;i < vec.size();i++){
            if(!(ans.size() == 0 && vec[i] == 0)) ans += vec[i] + '0';
        }
         return ans.size() == 0? "0":ans;
    }
};

分治乘法

Karatsuba乘法

image.png

首先把数(x,y)给变成偶数位,不是的在高位补0。而x或者y可以分为为2个数,即前半部分a和后半部分b。 也就是: 10n2a+b=x10^{\frac{n}{2}} * a + b = x。同理10n2c+d=y10^{\frac{n}{2}} * c + d = y

所以算xyx * y我们算(10n2a+b)(10n2c+d)(10^{\frac{n}{2}} * a + b )* ( 10^{\frac{n}{2}} * c + d)即可,进一步化简为10nac+10n2(ad+bc)+bd10^n * ac + 10^{\frac{n}{2}}(ad + bc) + bd

从这里看,我们需要求ac,ad,bc,bd。实际上,我们不需要ad和bc需要的是(ad + bc),那么求:(a+b)(b+c)acbd=ad+bc(a + b) * (b + c) - ac - bd = ad + bc,那么求这一个式子就好了。

代码:

class Solution {
public:
    int to_int(string num){
        int ans = 0;
        for(int i = 0;i < num.size();i++) ans = ans*10 + (num[i] - '0');
        return ans;
    }
    string multiply(string num1, string num2) {
       while(num1.size() > num2.size()) num2 = '0' + num2;
       while(num1.size() < num2.size()) num1 = '0' + num1;
       int n = num1.size();
       if(n == 1) return to_string((num1[0] - '0') * (num2[0] - '0'));
       int mi = 1;
       for(int i = 0;i < n/2;i++) mi *= 10;
       int a = to_int(num1);
       int b = a % mi;
       int c = to_int(num2);
       int d = c % mi;
       a /= mi,c /= mi;
       int ac = to_int(multiply(to_string(a),to_string(c)));
       int bd = to_int(multiply(to_string(b),to_string(d)));
       int abcd = to_int(multiply(to_string(a + b),to_string(c + d))) - ac - bd;
       return to_string(mi*mi*ac + mi*abcd + bd);
    }
};

很遗憾,这个方法在leetcode题目上会出错。改成long long依旧不可

Toom-Cook乘法

Toom Cook也是基于分而治之的算法,Toom Cook-k算法就是指将乘数分别分为固定大小的k组进行计算的算法。Toom Cook算法可以当做Karatsuba算法的泛化版本,会使用并不难,但是想要理解为什么要这么操作是有难度的。leetcode的对应题目,大家写几个例子就行,别提交了,大概率过不了。

为了方便我们来探究Toom-Cook-3,同上把x,y分为3份(sizeof(x)/3 = t),三份分别为x1,x2,x3x_1,x_2,x_3,y1,y2,y3y_1,y_2,y_3:

x=x2t2+x1t+x0x = x_{2} * t^2 + x_{1} * t +x_{0}
y=y2t2+y1t+y0y = y_{2} * t^2 + y_{1} * t +y_{0}

同理:

xy=(x2t2+x1t+x0)(y2t2+y1t+y0)=x0y0+x0y1t1+x0y2t2+x1y0t1+x1y1t2+x1y2t3+x2y0t2+x2y1t3+x2y2t4=x0y0+x0y1+x1y0t1+(x0y2+x2y0+x1y1)t2+(x1y2+x2y1)t3+x2y2t4=w0+w1t1+w2t2+w3t3+w4t4x * y = (x_{2} * t^2 + x_{1} * t +x_{0}) * (y_{2} * t^2 + y_{1} * t +y_{0}) \newline = x_0 * y_0 + x_0 * y_1 * t^1 + x_0 * y_2 * t^2 + x_1 * y_0 * t^1 + x_1 * y_1 * t^2 \newline + x_1 * y_2 * t^3 + x_2 * y_0 * t^2 + x_2 * y_1 * t^3 + x_2 * y_2 * t^4 \newline = x_0 * y_0 +(x_0 * y_1 + x_1 * y_0)* t^1 + ( x_0 * y_2 + x_2 * y_0 + x_1 * y_1) * t^2 \newline + (x_1 * y_2 + x_2 * y_1) * t^3 + x_2 * y_2 *t^4 \newline = w_0 + w_1 * t^1 + w_2 * t^2 + w_3 * t^3 + w_4 * t^4

那么我们令:

a=x0y0a = x_0 * y_0
b=(x0+x1+x2)(y0+y1+y2)b = (x_0 + x_1 + x_2) * (y_0 + y_1 + y_2)
c=(x0x1+x2)(y0y1+y2)c = (x_0 - x_1 + x_2) * (y_0 - y_1 + y_2)
d=(x0+2x1+4x2)(y0+2y1+4y2)d = (x_0 + 2x_1 + 4x_2) * (y_0 + 2y_1 + 4y_2)
e=(x02x1+4x2)(y02y1+4y2)e = (x_0 - 2x_1 + 4x_2) * (y_0 - 2y_1 + 4y_2) \newline

可以得到:

w0=aw_0 = a
w1=(8b8cd+e)/12w_1 = (8b - 8c - d + e)/12
w2=(30a+16b+16cde)/24w_2 = (-30a + 16b + 16c - d - e)/24
w3=(2b2c+de)/24w_3 = (-2b - 2c + d - e)/24
w4=(6a4b4c+d+e)/24w_4 = (6a - 4b - 4c + d + e)/24

就可以得到结果:

xy=w0+w1t1+w2t2+w3t3+w4t4x * y = w_0 + w_1 * t^1 + w_2 * t^2 + w_3 * t^3 + w_4 * t^4
class Solution {
public:
    int to_int(string num){
        int ans = 0;
        for(int i = 0;i < num.size();i++) ans = ans*10 + (num[i] - '0');
        return ans;
    }
    string multiply(string num1, string num2) {
       while(num1.size() > num2.size()) num2 = '0' + num2;
       while(num1.size() < num2.size()) num1 = '0' + num1;
       int n = num1.size();
       if(n < 3) return to_string(to_int(num1) * to_int(num2));
       int mi = 1;
       for(int i = 0;i < n/3;i++) mi *= 10;
       int x0 = to_int(num1) / mi;
       int x1 = (to_int(num1) % mi) / mi;
       int x2 = (to_int(num1) / mi) / mi;
       int y0 = to_int(num2) / mi;
       int y1 = (to_int(num2) % mi) / mi;
       int y2 = (to_int(num2) / mi) / mi;
       int a = to_int(multiply(to_string(x0),to_string(y0)));
       int b = to_int(multiply(to_string(x0 + x1 + x2),to_string(y0 + y1 +y2)));
       int c = to_int(multiply(to_string(x0 - x1 + x2),to_string(y0 - y1 +y2)));
       int d = to_int(multiply(to_string(x0 + 2 * x1 + 4 * x2),to_string(y0 + 2 * y1 + 4 * y2)));
       int e = to_int(multiply(to_string(x0 - 2 * x1 + 4 * x2),to_string(y0 - 2 * y1 + 4 * y2)));
       return to_string(a + ((8 * b - 8 * c - d + e) / 12 * mi) + ((-30 * a + 16 * b + 16 * c - d - e) / 24 * mi * mi) + (-2 * b - 2 * c +d - e) / 24 * mi * mi * mi + (6 * a - 4 * b - 4 * c + d + e) / 24 * mi * mi * mi * mi);
    }
};

快速傅里叶变换FFT

快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。

FFT(Fast Fourier Transformation ) 是离散傅氏变换(DFT)的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。

FFT其实是一个用O(nlog2n)O(nlog_2 n)的时间将一个用系数表示的多项式转换成它的点值表示的算法。多项式的系数表示和点值表示可以互相转换。

系数表示法:

一个n-1次n项多项式f(x)可以表示为: i=0n1aixi\sum_{i=0}^{n-1} a_i x^i 也可以用每一项的系数f(x)来表示即: f(x)={a0,a1,a2....,an1}f(x) = \{ a_0,a_1,a_2....,a_{n-1} \} 这就是系数表示法。

点值表示法:

  1. 把多项式放到平面直角坐标系里,把n个x都给带入,每一个x对应一个y。
  2. 把n个点连起来,如果n各不相同,那么n个点可以解出一个n个待定系数的方程组,如:f(x)=ax2+bx+c f(x) = ax^2 + bx + c 需要三个点就可以求出结果。
  3. 所有的点表示出来可以用:f(x)={{x0,f(x0)},{x1,f(x1)},.....,{xn1,f(xn1)}}f(x) = \{ \{ x_0,f(x_0) \},\{ x_1,f(x_1) \},.....,\{ x_{n-1},f(x_{n-1}) \} \}来表示

高精度乘法下两种多项式表示法的区别:

对于两个用系数表示的多项式,我们把它们相乘,类似上面的分治乘法,就是这样直接乘起来的,复杂度为O(n2)O(n^2)

点值表示的多项式相乘,只需要O(n)O(n)的时间。但是朴素的系数表示法转点值表示法还是O(n2)O(n^2)。快速傅里叶变换能达到O(nlog2n)O(nlog_2 n),快速傅里叶逆变换的复杂度O(n)O(n)

设两个点值多项式分别为:

f(x)={{x0,f(x0)},{x1,f(x1)},.....,{xn1,f(xn1)}}f(x) = \{ \{ x_0,f(x_0) \},\{ x_1,f(x_1) \},.....,\{ x_{n-1},f(x_{n-1}) \} \}

g(x)={{x0,g(x0)},{x1,g(x1)},.....,{xn1,g(xn1)}}g(x) = \{ \{ x_0,g(x_0) \},\{ x_1,g(x_1) \},.....,\{ x_{n-1},g(x_{n-1}) \} \}

设它们的乘积是h(x),那么:

h(x)={{x0,f(x0)g(x0)},{x1,f(x1)g(x1)},.....,{xn1,f(xn1)g(xn1)}}h(x) = \{ \{ x_0,f(x_0) * g(x_0) \},\{ x_1,f(x_1) * g(x_1) \},.....,\{ x_{n-1},f(x_{n-1}) * g(x_{n-1}) \} \}

从这里看时间复杂度只有O(n),但是就如同上面说的事实不是这样的。

离散傅里叶变换(DFT):

复数:我们把形如a+bi(a,b均为实数)的数称为复数,其中a称为实部,b称为虚部,i称为虚数单位。当虚部等于零时,这个复数可以视为实数;当z的虚部不等于零时,实部等于零时,常称z为纯虚数。复数域是实数域的代数闭包,也即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入,经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作,此概念逐渐为数学家所接受。(中学知识就不多说了)

复数的运算:复数不像点或向量,它和实数一样可以进行四则运算。设两个复数分别为:z1=a+biz_1 = a +bi,z2=c+diz_2 = c + di那么:

z1+z2=(a+c)+(c+d)i z_1 + z_2 = (a + c) + (c + d)i

z1z2=(acbd)+(ad+bc)iz_1 z_2 = (ac - bd) + (ad + bc)i

复数相加也满足平行四边形法则:

image.png

即AB + AD = AC,令αk(k=1,2..n) \alpha_k(k=1,2..n)为复数相对应向量的模长,θk(k=1,2..n)\theta_k(k=1,2..n)为负数对应的极角。那么复数相乘是:(α1,θ1)(α2,θ2)=(α1α2,θ1+θ2)(\alpha_1,\theta_1) * (\alpha_2,\theta_2) = (\alpha_1\alpha_2,\theta_1+ \theta_2)

离散傅里叶变换(DFT):对于任意系数多项式转点值,当然可以随便取任意n个x值代入计算。但是暴力计算需要O(n2)O(n^2)的时间。(这里n一定是2的整数次幂) 其实可以代入一组神奇的x,代入以后不用做那么多的次方运算 这些x当然不是乱取的,而且取这些x值应该就是傅里叶的主意了

考虑一下,如果我们代入一些x,使每个x的若干次方等于1,我们就不用做全部的次方运算了,±1和±i都是可以的,但只有这四个数远远不够。

这个圆圈上面的点都可以做到 image.png

以原点为圆心,画一个半径为1的单位圆,那么单位圆上所有的点都可以经过若干次次方得到1,傅里叶说还要把它给n等分了,比如此时n = 8。

image.png

橙色的点即为n = 8时要取的点,从(1,0)开始,逆时针从0开始编号,标到7号,记编号为k的点代表的复数的值为ωnk ω_{n}^k,那么由模长相乘,极角相加可知(ωn1)k=ωnk(ω_{n}^1)^k = ω_{n}^k 其中ωn1ω_{n}^1称为n次单位根,而且每一个ω都可以求出。

ωnk=coskn2π+sinkn2πi ω_{n}^k = \cos \frac{k}{n} 2\pi + \sin \frac{k}{n}2\pi i

也就是说cos2kn2π+sin2kn2π=1 \cos^2 \frac{k}{n} 2\pi + \sin^2 \frac{k}{n}2\pi = 1,具体见下图:

image.png

其中sin2θ+cos2θ=1\sin^2 \theta + \cos^2 \theta = 1

那么ωn0ω_{n}^0,ωn1ω_{n}^1...ωnn1ω_{n}^{n-1},就是我们要带入的x0x_0,x1x_1....xn1x_{n-1}

推FFT的过程中需要用到ω\omega的一些性质:

ωnk=ω2n2k\omega_{n}^k = \omega_{2n}^2k,它们表示的点(或向量)表示的复数是相同的。

证明: ωnk=coskn2π+sinkn2πi=cos2k2n2π+sin2k2n2πi=ω2n2k\omega_{n}^{k} = \cos \frac{k}{n} 2\pi + \sin \frac{k}{n}2\pi i = \cos \frac{2k}{2n} 2\pi + \sin \frac{2k}{2n}2\pi i = \omega_{2n}^{2k}

ωnk+n2=ωnk \omega_{n}^{k + \frac{n}{2}} = - \omega_{n}^{k},它们表示的点关于原点对称,所表示的复数实部相反,所表示的向量等大反向。

证明:ωnk+n2=cosk+n2n2π+sink+n2n2πi=coskn2π+sinkn2πi=ωnk ω_{n}^{k + \frac{n}{2}} = \cos \frac{k + \frac{n}{2}}{n} 2\pi + \sin \frac{k + \frac{n}{2}}{n}2\pi i = - \cos \frac{k}{n} 2\pi + \sin \frac{k}{n}2\pi i = -ω_{n}^k

ωn0=ωnn\omega_{n}^{0} = \omega_{n}^{n}

都等于1或1 + 0i。

快速傅里叶变换(FFT):

虽然DFTωω\omegaω作为代入多项式的x值,代入这类特殊x值法的变换叫做DFT,但是DFT还是要代入单位根暴力计算,复杂度依旧是O(n2)O(n^2)

A(x)=i=0n1aixiA(x) = \sum_{i=0}^{n-1} a_i x^i,然后按照ii的奇偶的不同,分成两部分A(x)=(a0+a2x2+....+an2xn2)+(a1+a3x2+....+an1xn2)xA(x)= (a_0 +a_2 x^2+....+a_{n-2} x^{n-2}) + (a_1 +a_3 x^2+....+a_{n-1} x^{n-2})x。这样可以设出来A1(x)A_1(x),A2(x)A_2(x),其中A1(x)=a0+a2x1+....+an2xn21A_1(x) = a_0 +a_2 x^1+....+a_{n-2} x^{\frac{n}{2} - 1}A2(x)=a1+a3x1+....+an1xn21A_2(x) = a_1 +a_3 x^1+....+a_{n-1} x^{\frac{n}{2} - 1},明显可知 A(x)=A1(x2)+A2(x2)x A(x) = A_1(x^2) + A_2(x^2)x

再设k<n2k < \frac{n}{2},把ωnk\omega_n^{k}作为x代入A(x),也就是A(ωnk)=A1((ωnk)2)+A2((ωnk)2)ωnkA(\omega_n^{k}) = A_1((\omega_n^{k})^2) + A_2((\omega_n^{k})^2)\omega_n^{k},进一步化简为A(ωnk)=A1(ωn2k)+A2(ωn2k)ωnk=A1(ωn2k)+A2(ωn2k)ωnkA(\omega_n^{k}) = A_1(\omega_{n}^{2k}) + A_2(\omega_{n}^{2k})\omega_n^{k} = A_1(\omega_{\frac{n}{2}}^{k}) + A_2(\omega_{\frac{n}{2}}^{k})\omega_n^{k}

同理对A(ωnk+n2)A(\omega_n^{k+\frac{n}{2}})化简:A(ωnk+n2)=A1(ωn2k+n)+A2(ωn2k+n)ωnk=A1(ωn2k)A2(ωn2k)ωnkA(\omega_n^{k+\frac{n}{2}}) = A_1(\omega_n^{2k+n}) + A_2(\omega_n^{2k + n})\omega_n^{k} = A_1(\omega_{\frac{n}{2}}^{k}) - A_2(\omega_{\frac{n}{2}}^{k})\omega_n^{k}

通过上面公式,我们可以知道A(ωnk+n2)A(\omega_n^{k+\frac{n}{2}})ωnk\omega_n^{k}最后只有尾部的A2(ωn2k)ωnk A_2(\omega_{\frac{n}{2}}^{k})\omega_n^{k}这一部分前面的符号是不同的,那么也就是说我们得到A2(ωn2k)A_2(\omega_{\frac{n}{2}}^{k})A1(ωn2k)A_1(\omega_{\frac{n}{2}}^{k})就可以求出A(ωnk+n2)A(\omega_n^{k+\frac{n}{2}})ωnk\omega_n^{k},这样就可以堪称一个分治的问题,也就是,每次都把一整个分成两部分把AA分成A1,A2A_1,A_2。然后再基于A1,A2A_1,A_2再继续去分。直到只有一个数,也就是长度等于1.就返回退到上一步。如果不能理解,可以去学一下递归。其实这里不会也没事,看程序,或者看博客后续内容。这里时间复杂度为O(nlog2n) O(nlog_2 n)

快速傅里叶逆变换(IFFT):

我们把两个多项式相乘 (也叫求卷积),做完两遍FFT也知道了积的多项式的点值表示,算完以后得到结果可我们平时用系数表示的多项式,点值表示没有意义,需要变回去,变回去的过程就是傅里叶逆变换。

刚考完数一DNA动了,范德蒙行列式卷积定理我可太会了,但是这是第二章,可以讲明白但是没必要,以后慢慢来。

看到路人黑的纸巾最后就写了个结论,我觉得现在这个阶段直接知道结论就行: 一个多项式在分治的过程中乘上单位根的共轭复数,分治完的每一项除以n即为原多项式的每一项系数

很多比如蝴蝶操作霍纳法则等,以后专门讨论多项式的时候再写博客。会参考算导的内容,但是现在刚刚开始,我觉得有点过了。

代码:

class Solution {
public:
const double PI = acos(-1.0);
//建立复数的数据结构
    struct complex{
        double re,mi;
        complex(double re = 0.0,double mi = 0.0){
            this->re = re;
            this->mi = mi;
        }
        complex operator-(const complex& other) {
            return complex(re - other.re,mi - other.mi);
        }
        complex operator+(const complex& other){
            return complex(re + other.re,mi + other.mi);
        }
        complex operator*(const complex& other){
            return complex(re * other.re - mi * other.mi,re * other.mi + mi * other.re);
        }
    };
//复现书上的fft函数
    vector<complex> fft(vector<complex> &a,bool invert){
        int n = a.size();
        if(n==1) return a;
        vector<complex> pe(n/2),po(n/2);
        for(int i = 0;2 * i < n; i++)
        {
            pe[i] = a[2*i];
            po[i] = a[2 * i + 1];
        }
        pe = fft(pe,invert);
        po = fft(po,invert);
        vector<complex> snowball(n);
        double aug = 2 * PI / n * (invert?-1:1);
        complex wn(cos(aug),sin(aug));
        complex w(1,0);
        for(int i = 0;i < n/2; i++){
            snowball[i] = pe[i] + w*po[i];
            snowball[n/2 + i] = pe[i] - w*po[i];
            w = wn * w;
        }
        return snowball;
    }
  string multiply(string num1, string num2) {
        if (num1 == "0" || num2 == "0") {
            return "0";
        }

        int len1 = num1.size();
        int len2 = num2.size();

        int n = 1;
        while (n < len1 + len2) {
            n = n << 1;
        }

        vector<complex> a(n);
        vector<complex> b(n);

        for (int i = len1 - 1; i >= 0; i--) {
            a[i] = complex(num1[len1 - 1 - i] - '0', 0);
        }

        for (int i = len2 - 1; i >= 0; i--) {
            b[i] = complex(num2[len2 - 1 - i] - '0', 0);
        }

        a = fft(a, false);
        b = fft(b, false);

        for (int i = 0; i < n; i++) {
            a[i] = a[i] * b[i];
        }

        a = fft(a, true);

        string snowball;
        int carrsonwball = 0;
        for (int i = 0; i < n; i++) {
            int sum = round(round(a[i].re) / n) + carrsonwball;
            carrsonwball = sum / 10;
            snowball += sum % 10 + '0';
        }

        if (carrsonwball > 0) {
            snowball += carrsonwball % 10 + '0';
        }

        int idx = snowball.size() - 1;
        while (snowball[idx] == '0' && idx > 0) {
            idx--;
        }

        snowball = snowball.substr(0, idx + 1);
        reverse(snowball.begin(), snowball.end());
        return snowball;
    }
};

大概记住过程即可,只打比赛记住就够了。这个代码可以跑通,直接抄就行。

中国剩余定理

简介:

在《孙子算经》中有这样一个问题:“今有物不知其数,三三数之剩二(除以3余2),五五数之剩三(除以5余3),七七数之剩二(除以7余2),问物几何?”这个问题称为“孙子问题”,该问题的一般解法国际上称为“中国剩余定理”。

解法:

  1. 找出三个数:从3和5的公倍数中找出被7除余1的最小数15,从3和7的公倍数中找出被5除余1 的最小数21,最后从5和7的公倍数中找出除3余1的最小数70。
  2. 用15乘以2(2为最终结果除以7的余数),用21乘以3(3为最终结果除以5的余数),同理,用70乘以2(2为最终结果除以3的余数),然后把三个乘积相加15∗2+21∗3+70∗215∗2+21∗3+70∗2得到和233。
  3. 用233除以3,5,7三个数的最小公倍数105,得到余数23,即233%105=23233%105=23。这个余数23就是符合条件的最小数。

我们将“孙子问题”拆分成几个简单的小问题,从零开始,试图揣测古人是如何推导出这个解法的。

  首先,我们假设n1n_1是满足除以3余2的一个数,比如2,5,8等等,也就是满足3k+2k>=03∗k+2(k>=0)的一个任意数。同样,我们假设n2n_2是满足除以5余3的一个数,n3n_3是满足除以7余2的一个数。

  有了前面的假设,我们先从n1n_1这个角度出发,已知n1n_1满足除以3余2,能不能使得n1+n2n_1+n_2的和仍然满足除以3余2?进而使得n1+n2+n3n_1+n_2+n_3的和仍然满足除以3余2?

  这就牵涉到一个最基本数学定理,如果有a%b=ca\%b=c,则有(a+kb)%b=c(kZ,k0)(a+k∗b)\%b = c (k∈Z,k≠0),换句话说,如果一个除法运算的余数为c,那么被除数与k倍的除数相加(或相减)的和(差)再与除数相除,余数不变。这个是很好证明的。

以此定理为依据,如果n2n_2是3的倍数,n1+n2n_1+n_2就依然满足除以3余2。同理,如果n3n_3也是3的倍数,那么n1+n2+n3n_1+n_2+n_3的和就满足除以3余2。这是从n1n_1的角度考虑的,再从n2n_2n3n_3的角度出发,我们可推导出以下三点:

  1. 为使n1+n2+n3n_1+n_2+n_3的和满足除以3余2,n2n_2n3n_3必须是3的倍数。
  2. 为使n1+n2+n3n_1+n_2+n_3的和满足除以5余3,n1n_1n3n_3必须是5的倍数。
  3. 为使n1+n2+n3n_1+n_2+n_3的和满足除以7余2,n1n_1n2n_2必须是7的倍数。

  因此,为使n1+n2+n3n_1+n_2+n_3的和作为“孙子问题”的一个最终解,需满足:

  1. n1n_1除以3余2,且是5和7的公倍数。
  2. n2n_2除以5余3,且是3和7的公倍数。
  3. n3n_3除以7余2,且是3和5的公倍数。

所以,孙子问题解法的本质是从5和7的公倍数中找一个除以3余2的数n1n_1,从3和7的公倍数中找一个除以5余3的数n2n_2,从3和5的公倍数中找一个除以7余2的数n3n_3,再将三个数相加得到解。在求n1n_1n2n_2n3n_3时又用了一个小技巧,以n1n_1为例,并非从5和7的公倍数中直接找一个除以3余2的数,而是先找一个除以3余1的数,再乘以2。也就是先求出5和7的公倍数模3下的逆元,再用逆元去乘余数:

a%b=ca \% b = c
(ak)%b=c+c...+c=ck(a * k)\% b = c + c...+c = c * k

最后,我们还要清楚一点,n1+n2+n3n_1+n_2+n_3只是问题的一个解,并不是最小的解。如何得到最小解?我们只需要从中最大限度的减掉掉3,5,7的公倍数105即可。道理就是前面讲过的定理“如果a%b=ca\%b=c,则有(akb)%b=c(a−k∗b)\%b=c。所以例子里n1+n2+n3%105(n_1+n_2+n_3)\%105就是最终的最小解。

这样就都得到了中国剩余定理的公式:

有一个数x模3的结果是2,模5的结果是3,模7的结果是2。这是问题给的条件。 设正整数两两互素,则同余方程组

                             

有整数解。并且在模下的解是唯一的,解为

                               

其中,而的逆元。 普通的中国剩余定理要求所有的互素,那么如果不互素呢,怎么求解同余方程组?

  这种情况就采用两两合并的思想,假设要合并如下两个方程:

  那么得到:

  我们需要求出一个最小的x使它满足:

  那么x1x_1x2x_2就要尽可能的小,于是我们用扩展欧几里得算法求出x1x_1的最小正整数解,将它代回a1+m1x1a_1+m_1x_1,得到x的一个特解x′,当然也是最小正整数解。

  所以x的通解一定是x′加上lcm(m1,m2)klcm(m_1,m_2)∗k,这样才能保证x模m1m_1m2m_2的余数是a1a_1a2a_2。由此,我们把这个x′当做新的方程的余数,把lcm(m1,m2)lcm(m_1,m_2)当做新的方程的模数。(这一段是关键

  合并完成:

具体代码以后再写,有点难。我找了找合适的题,是noi省选,没法做。等我能力强大了再给大家写题解吧。我太菜了。呜呜呜~

Furer’s algorith

2007年,宾夕法尼亚州立大学数学家Martin Fürer提出逼近O(nlog2n)O(nlog_2n)Fürer算法,打破没有进展的僵局,过去十年乘法算法不断改善,无限接近O(nlog2n)O(nlog_2n)

都是很牛的算法了,感兴趣可以看维基百科,兄弟们我crt还写不出来,这个就先咕咕了。

参考

[1] 大数乘法的深入讨论一(模拟小学手算法),楠子小先生

[2] <<算法详解(卷1)—算法基础>>,[美]蒂姆·拉夫加登

[3] 十分简明易懂的FFT(快速傅里叶变换),路人黑的纸巾

[4] 中国剩余定理学习笔记,MashiroSky