大整数算法[10] Comba乘法（实现）
★ 引子
上一篇文章讲了 Comba 乘法的原理，这次来讲讲如何实现。为了方便移植和充分发挥不同平台下的性能，暂时用了三种不同的实现方式：
1、单双精度变量都有的情况。
2、只有单精度变量的情况。
3、可以使用内联汇编的情况。
前面已经介绍了 Comba 乘法的原理和实现思路，为了方便，再把它贴到这里：
计算 c = a * b，c0，c1，c2 为单精度变量。
1. 增加 c 到所需要的精度，并且令 c = 0，c->used = a->used + b->used。
2. 令 c2 = c1 = c0 = 0。
3. 对于 i 从 0 到 c->used - 1 循环，执行如下操作：
3.1 ty = BN_MIN(i, b->used - 1)
3.2 tx = i - ty
3.3 k = BN_MIN(a->used - tx, ty + 1)
3.4 三精度变量右移一个数位：(c2 || c1 || c0) = (0 || c2 || c1)
3.5 对于 j 从 0 到 k - 1 之间执行循环，计算：
(c2 || c1 || c0) = (c2 || c1 || c0) + a(tx + j) * b(ty - j)
3.6 c(i) = c0
4. 压缩多余位，返回 c。
上面所说的三种不同实现方式，主要体现在 3.5 中的循环。 Comba 乘法的实现代码如下：（暂时不考虑 x = x * y 这种输入和输出是同一个变量的情况）
static void bn_mul_comba(bignum *z, const bignum *x, const bignum *y)
{bn_digit c0, c1, c2;bn_digit *px, *py, *pz;size_t nc, i, j, k, tx, ty;pz = z->dp;nc = z->used;c0 = c1 = c2 = 0;for(i = 0; i < nc; i++){ty = BN_MIN(i, y->used - 1);tx = i - ty;k = BN_MIN(x->used - tx, ty + 1);px = x->dp + tx;py = y->dp + ty;c0 = c1;c1 = c2;c2 = 0U;//Comba 32for(j = k; j >= 32; j -= 32){COMBA_INITCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_STOP}//Comba 16for(; j >= 16; j -= 16){COMBA_INITCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_STOP}//Comba 8for(; j >= 8; j -= 8){COMBA_INITCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_STOP}//Comba 4for(; j >= 4; j -= 4){COMBA_INITCOMBA_MULADDC COMBA_MULADDCCOMBA_MULADDC COMBA_MULADDCCOMBA_STOP}//Comba 1for(; j > 0; j--){COMBA_INITCOMBA_MULADDCCOMBA_STOP}*pz++ = c0;} //for ibn_clamp(z);
}
在 3.5 步中的内循环，乘法的计算按照 1，4，8，16，32 的步进展开，展开的过程需要比较大的空间，但是减少了循环控制的开销，所以可以节省大量时间。执行乘法的关键代码都定义在 COMBA_INIT，COMBA_MULADDC 和 COMBA_STOP 这三个宏内。COMBA_INIT 用于变量定义或者初始化，COMBA_MULADDC 用于计算单精度乘法和累加，COMBA_STOP 用于存储当前计算结果。另外，这里假设 z 的精度已经足够。
★ 单双精度变量都有的情况
这种情况下，因为 bn_digit 和 bn_udbl 同时有定义，所以是最容易实现的一种。三个宏的定义如下：
#define COMBA_INIT \
{ \bn_udbl r; \#define COMBA_MULADDC \\r = (bn_udbl)(*px++) * (*py--) + c0; \c0 = (bn_digit)r; \r = (bn_udbl)c1 + (r >> biL); \c1 = (bn_digit)r; \c2 += (bn_digit)(r >> biL); \#define COMBA_STOP \
}
在乘法器中，定义双精度变量 r 存储单精度乘法结果。在计算中，先计算单精度乘法，然后在与三精度变量 c2||c1||c0 相加。因为 c0，c1，c2 都是单精度数，为了避免溢出，先把累加计算的结果放到 r 上，然后提取低半部分作为当前数位的值。这里暂时不需要存储计算结果，所以 COMBA_STOP 没有什么操作。
★ 只有单精度变量的情况
这种情况下，bn_udbl 没有定义，每个数被拆成高半部分和低半部分，要执行 4 次单精度乘法，所以计算会比较复杂。
#define COMBA_INIT \
{ \bn_digit a0, a1, b0, b1; \bn_digit t0, t1, r0, r1; \#define COMBA_MULADDC \\a0 = (*px << biLH) >> biLH; \b0 = (*py << biLH) >> biLH; \a1 = *px++ >> biLH; \b1 = *py-- >> biLH; \r0 = a0 * b0; \r1 = a1 * b1; \t0 = a1 * b0; \t1 = a0 * b1; \r1 += (t0 >> biLH); \r1 += (t1 >> biLH); \t0 <<= biLH; \t1 <<= biLH; \r0 += t0; \r1 += (r0 < t0); \r0 += t1; \r1 += (r0 < t1); \c0 += r0; \c1 += (c0 < r0); \c1 += r1; \c2 += (c1 < r1); \#define COMBA_STOP \
}
a1，a0 分别存储 x 的高半部分和低半部分，b1，b0 分别存储 y 的高半部分和低半部分，这两个操作可以通过将 x 和 y 向左向右移动半个数位得到。biLH 前面提到过，这里在说一下：biLH 是每个数位的比特大小的一半。上面这个乘法计算比较复杂，所以我先讲讲原理。
设 f(x)，g(x) 分别表示单精度数 a 和 b，要计算 c = a * b。
a 和 b 用 f(x) 和 g(x) 表示：f(x) = a1 * x + a0，g(x) = b1 * x + b0，其中这里的 x 就是 2 ^ (n / 2)。
那么 h(x) = f(x) * g(x) = a1 * b1 * (x ^ 2) + (a1 * b0 + a0 * b1) * x + a0 * b0。这里默认 a0，a1，b0，b1 都是半精度的变量，所以他们的乘积可以用一个单精度变量存储。
令 t0 = a1 * b0 = p * x + q，t1 = a0 * b1 = r * x + s，r0 = a0 * b0，r1 = a1 * b1，其中 p，q，r，s 是半精度数，那么 h(x) 可以改写成：
h(x) = a1 * b1 * (x ^ 2) + (p * x + q + r * x + s) * x + a0 * b0
= a1 * b1 * (x ^ 2) + ((p + r) * x + (q + s)) * x + a0 * b0
= a1 * b1 * (x ^ 2) + (p + r) * (x ^ 2) + (q + s) * x + a0 * b0
= (a1 * b1 + p + r) * (x ^ 2) + (q + s) * x + a0 * b0
= (r1 + p + r) * (x ^ 2) + (q + s) * x + r0
所以，要计算 (r1 + p + r) * (x ^ 2)，需要提取 t0，t1 的高半部分和 r1 相加（因为 r1 和 p，r 都是同类项，所以直接计算 x^2 的系数即可）。计算剩余项 (q + s) * x 和 r0 的和，需要将 t0 和 t1 的低半部分提取出来，左移半个数位后和 r0 相加（r0 是单精度数，而 q 和 s 是半个精度的，乘以 x 要左移半个数位）。计算 (q + s) * x + r0 后结果可能会溢出，表明产生了进位，所以还要把进位要传递到 (r1 + p + r) * (x ^ 2)，最终乘积 c 就可以用 r0，r1 组成的双精度变量 r1||r0 表示。最后将乘积 c 和三精度变量 c2||c1||c0 累加即可。
★ 使用内联汇编的情况
注：本小节涉及到 x86 的汇编，如果对汇编不了解的话，可以跳过，免得浪费心情。对于 C 的内联汇编，细节暂时不说了，参考下面两篇文章：
【GCC 和 VC 的内联汇编】 https://github.com/1184893257/simplelinux/blob/master/inlineasm.md
【Linux 中 x86 的内联汇编】 http://www.ibm.com/developerworks/cn/linux/sdk/assemble/inline/
如果编译环境可以使用内联汇编，则使用汇编指令执行乘法会加快计算速度。由于我只接触过 x86 平台的汇编，所以暂时只考虑 x86 平台，ARM，PowerPC 或者 MIPS 之类的架构暂时不管。 在 x86 环境下，需要考虑到 GCC 和 VC 环境下的内联汇编，毕竟他们的语法差别很大。
A. VC x86 环境下的内联汇编：（和 GCC 的比起来要简单得多）
#define COMBA_INIT \
{ \__asm mov esi, px \__asm mov edi, py \#define COMBA_MULADDC \\__asm lodsd \__asm mov ebx, [edi] \__asm mul ebx \__asm add c0, eax \__asm adc c1, edx \__asm adc c2, 0 \__asm sub edi, 4 \#define COMBA_STOP \\__asm mov px, esi \__asm mov py, edi \
}
1. 首先把指针 px 和 py 的地址用 MOV 指令分别送到源变址寄存器 ESI 和目的变址寄存器 EDI。
2. 执行单精度乘法和累加：
A. 执行 LODSD 指令，把 ESI 寄存器指向的数据段中某单元的内容（也就是 x 的某个数位的值）送到 EAX 寄存器中，并根据方向标志和数据类型修改 ESI 寄存器的内容，具体的操作是将 ESI 中 的地址增加 4，相当于 C 中的 px++。一般情况下，方向标志 DF = 0， ESI 中的地址增加。如果 DF = 1，要执行 CLD 指令把 DF 置为 0，不过一般情况下不会出现这种情况的。
B. 用 MOV 指令把 EDI 寄存器所指向的数据段中某单元的内容（也就是 y 的某个数位的值）送到 EBX 寄存器中。
C. 执行 MUL 指令执行单精度乘法：（EDX，EAX）= EAX * EBX，计算 EAX 和 EBX 的乘积，结果的高位放在 EDX 寄存器中，结果的低位放在 EAX 寄存器中。
D. 执行 ADD 加法指令，将乘积的低位累加到三精度变量的最低位 c0 上。
E. 执行 ADC 带进位的加法指令，将乘积的高位和上一步加法产生的进位累加到三精度变量的第二个数位 c1 上。
F. 执行 ADC 带进位的加法指令，将剩余的进位加到三精度变量的最高位 c2 上。
G. 执行 SUB 减法指令，将 EDI 中的地址往前挪 4 字节，相当于 C 中的 py--。
3. 存储计算结果：用 MOV 指令把 ESI 和 EDI 中的地址送回到指针 px 和 py 中，之所以要做这一步，是因为在循环中，循环控制可能会修改 CPU 寄存器的值，如果在计算结束后不存储 ESI 和 EDI 的地址值，那么地址可能就会因为寄存器被修改而丢失。
B. GCC x86 环境下的内联汇编：（比较复杂）
#define COMBA_INIT \
{ \asm \( \"movl %5, %%esi \n\t" \"movl %6, %%edi \n\t" \#define COMBA_MULADDC \\"lodsl \n\t" \"movl (%%edi), %%ebx \n\t" \"mull %%ebx \n\t" \"addl %%eax, %2 \n\t" \"adcl %%edx, %3 \n\t" \"adcl $0, %4 \n\t" \"subl $4, %%edi \n\t" \#define COMBA_STOP \\"movl %%esi, %0 \n\t" \"movl %%edi, %1 \n\t" \:"=m"(px),"=m"(py),"=m"(c0),"=m"(c1),"=m"(c2) \:"m"(px),"m"(py) \:"%eax","%ebx","%ecx","%edx","%esi","%edi" \); \
}
这段 GCC 的内联汇编指令和上面的 VC 内联汇编指令操作是一样的，只是语法不同而已，具体的语法请自行参考上面的资料或者 Google 搜索 :)
C. GCC x64 环境下的内联汇编：（和上面的 GCC x86 类似） 【补充】
考虑到在 64 bit 的 GCC 编译器下，可以使用 x64 的内联汇编，所以又加了这一部分代码。GCC x64 环境下的内联汇编和 GCC x86 的比较类似，只是寄存器和指令有些不同。细节可以参考这里：http://www.dutor.net/index.php/2013/11/att-syntax-and-gcc-inlined-assembly/
#define COMBA_INIT \
{ \asm \( \"movq %5, %%rsi \n\t" \"movq %6, %%rdi \n\t" \#define COMBA_MULADDC \\"lodsq \n\t" \"movq (%%rdi), %%rbx \n\t" \"mulq %%rbx \n\t" \"addq %%rax, %2 \n\t" \"adcq %%rdx, %3 \n\t" \"adcq $0, %4 \n\t" \"subq $8, %%rdi \n\t" \#define COMBA_STOP \\"movq %%rsi, %0 \n\t" \"movq %%rdi, %1 \n\t" \:"=m"(px),"=m"(py),"=m"(c0),"=m"(c1),"=m"(c2) \:"m"(px),"m"(py) \:"%rax","%rbx","%rcx","%rdx","%rsi","%rdi" \); \
}
要注意的是：在 GCC x64 的环境下，bn_digit 是 64 bit的，也就是 8 个字节，所以执行完乘法后，py 的地址应该往前挪动 8 个字节，因此 COMBA_MULADDC 中最后的减法指令是把地址减去 8，而不是减 4。
★ 总结
Comba 乘法大概就讲这么多吧。由于 Comba 乘法的时间复杂度仍然是 O(n^2)，所以当输入的规模 n 越大时，所需的时间仍然会急剧增加。下一篇文章将讲讲如何使用分治的方式降低乘法的时间复杂度。
【回到本系列目录】
版权声明
原创博文，转载必须包含本声明，保持本文完整，并以超链接形式注明作者Starrybird和本文原始地址：http://www.cnblogs.com/starrybird/p/4441022.html
转载于:https://www.cnblogs.com/starrybird/p/4441022.html
大整数算法[10] Comba乘法（实现）相关推荐
- 大整数算法[09] Comba乘法（原理）
★ 引子 原本打算一篇文章讲完,后来发现篇幅会很大,所以拆成两部分,先讲原理,再讲实现.实现的话相对复杂,要用到内联汇编,要考虑不同平台等等. 在大整数计算中,乘法是非常重要的,因为 ...
- 大整数乘法（Comba 乘法 （Comba Multiplication）原理）
Comba 乘法以(在密码学方面)不太出名的 Paul G. Comba 得名.上面的笔算乘法,虽然比较简单, 但是有个很大的问题:在 O(n^2) 的复杂度上进行计算和向上传递进位,看看前面的那个竖 ...
- javascript乘法和加法_js 大整数加法、乘法、除法
有一定的编程经验的人都知道,当我们对数据操作的时候,若数据过大,计算机将这个大数从十进制的转为二进制的,是没有那个类型的放的了的,因此,我们经常将一个大数转化为字符串来操作.接下来我将演示大整数的加法 ...
- 山东大学软件工程应用与实践——GMSSL开源库（十）——重要的大整数
2021SC@SDUSC 文章标题 BIGNUM的结构 大整数的加减法(绝对值加减) 绝对值加减 带符号加减 大整数的乘法 经典乘法 递归和comba乘法 小结 无论是SM2.SM9等国密算法还是以R ...
- 数据挖掘十大经典算法（详解）
数据挖掘十大经典算法 一. C4.5 C4.5算法是机器学习算法中的一种分类决策树算法 ...
- [转]数据挖掘十大经典算法
国际权威的学术组织the IEEE International Conference on Data Mining (ICDM) 2006年12月评选出了数据挖掘领域的十大经典算法:C4.5, k-M ...
- DataMining学习2_数据挖掘十大经典算法
数据挖掘十大经典算法 一. C4.5 C4.5算法是机器学习算法中的一种分类决策树算法,其核心算法是ID3 算法. C4.5算法继承了ID3算法的优点,并在以下几方面对ID3算法进行了改进: ...
- 【转】十大数据挖掘算法
数据挖掘十大经典算法 转自:https://www.cnblogs.com/lcchuguo/p/4874288.html 一. C4.5 C4.5算法是机器学习算法中的一种分类决策树算法,其核心 ...
- python两数相乘代码_Python 实现大整数乘法算法的示例代码
我们平时接触的长乘法,按位相乘,是一种时间复杂度为 O(n ^ 2) 的算法.今天,我们来介绍一种时间复杂度为 O (n ^ log 3) 的大整数乘法(log 表示以 2 为底的对数). 介绍原理 ...
- 计算机算法设计与分析 大整数乘法
大整数乘法 问题描述 求两个不超过200位的非负整数的积. 输入形式 有两行,每行是一个不超过200位的非负整数,没有多余的前导0. 输出形式 一行,即相乘后的结果.结果里不能有多余的前导0,即如果结 ...
最新文章
- 数字图像缩放之最近邻插值与双线性插值处理效果对比
- CV：基于Keras利用训练好的hdf5模型进行目标检测实现输出模型中的脸部表情或性别的gradcam(可视化)
- [Medical Image Processing] 2. Image Binary -【OTSU Algorithm Entropy Method】
- 堆排序（基于完全二叉树）
- php 递归中的全局变量,PHP中递归的实现实例详解
- 微信开发必看，使用.Net Core 开发微信跨平台应用
- java登陆session用法_Java web 登录 使用shiro和基于session的方式有何不同？
- python numpy数组画图_python – Numpy和Pyplot的条件绘图
- 2018程序员必备的4大网站推荐（代码、编程、SQL等）
- centos7.0KVM虚拟化
- Flutter的SnackBar
- Qt编写onvif搜索和云台控制工具源码
- php 公众号授权登录,微信公众号授权登录
- YS-M8C自助终端主板RK3288详解和说明
- 放弃75W年薪，回老家当公务员，提离职被领导教育，网友：leader嫉妒了
- 时统ptp_IEEE1588对时系统,PTP校时模块,PTP时钟服务器
- 互联网医院系统软件开发|互联网医院管理系统开发的好处
- Mezzanine Markdown的问题
- 调查问卷的反向计分题数据如何录入？
- c语言程序设计分值,计算机考试题分值分布
热门文章
- docker-compose部署EFKL，存储，过滤laravel日志文件
- java.lang.UnsatisfiedLinkError: Library not loaded: @rpath/libo.dylib Reason: image not found
- 全网首发：跟踪分析This parser does not support specification “null“ version “null“
- 文字处理技术：终于明白了压缩要怎么做
- 全网首发：神奇的WORD文字渲染效果，18个字符一组，每个字符渲染效果不同
- 安卓GLSurfaceView使用简单范例
- 解决办法：在查找预编译头使用时跳过
- 管理感悟：如何处理不干活的主管
- 求最大公约数欧几里得算法
- access html导出,AccessToFile