MbedTLS中的Montgomery算法实现解析(三)
MbedTLS中的Montgomery算法实现解析(三)
在前面的两篇博客中,我们详细讨论了Montgomery算法的数学原理以及MbedTLS中是如何快速求得-N-1(mod n)。在这一篇博客中,我们将探讨MbedTLS中的蒙哥马利约减(Montgomery reduction)和蒙哥马利模乘(Montgomery multiplication)是如何实现的。
在阅读本博客之前,建议先掌握MbedTLS的大数表示方式
蒙哥马利约减Montgomery reduction
在我们的数学推导中,蒙哥马利约减是蒙哥马利模乘的基础,也就是说我们首先有:
X ≡ a × R-1 (mod n)
然后才有:
X = a × R × b × R × R-1 (mod n)
但是MbedTLS将这两者结合了起来,在MbedTLS中定义的蒙哥马利模乘式为:
X = A × B × R-1 (mod n)
既包含了模乘(Montgomery multiplication),又包含了约减(Montgomery reduction),这样的好处就是约减可以写成模乘的一种特殊情况,即当B=1时:
X = A × B × R-1 (mod n) = A × R-1 (mod n)
因此在MbenTLS中,约减就是模乘的一种特殊情况。
下面我们来看代码。
// Montgomery reduction: A = A * R^-1 mod N
static int mpi_montred( mbedtls_mpi *A, const mbedtls_mpi *N, mbedtls_mpi_uint mm, const mbedtls_mpi *T )
{mbedtls_mpi_uint z = 1;mbedtls_mpi U;U.n = U.s = (int) z;U.p = &z;return( mpi_montmul( A, &U, N, mm, T ) );
}
可以看到这里初始化了一个值为1的大数变量U,随后调用了蒙哥马利模乘函数mpi_montmul
蒙哥马利模乘Montgomery multiplication
下面我们来看看这个蒙哥马利模乘函数mpi_montmul:
// Montgomery multiplication: A = A * B * R^-1 mod N
static int mpi_montmul( mbedtls_mpi *A, const mbedtls_mpi *B, const mbedtls_mpi *N, mbedtls_mpi_uint mm,const mbedtls_mpi *T )
{size_t i, n, m;mbedtls_mpi_uint u0, u1, *d;if( T->n < N->n + 1 || T->p == NULL )return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );memset( T->p, 0, T->n * ciL );d = T->p;n = N->n;m = ( B->n < n ) ? B->n : n;for( i = 0; i < n; i++ ){/** T = (T + u0*B + u1*N) / 2^biL*/u0 = A->p[i];u1 = ( d[0] + u0 * B->p[0] ) * mm;mpi_mul_hlp( m, B->p, d, u0 );mpi_mul_hlp( n, N->p, d, u1 );*d++ = u0; d[n + 1] = 0;}memcpy( A->p, d, ( n + 1 ) * ciL );if( mbedtls_mpi_cmp_abs( A, N ) >= 0 )mpi_sub_hlp( n, N->p, A->p );else/* prevent timing attacks */mpi_sub_hlp( n, A->p, T->p );return( 0 );
}
函数有五个形参:
mbedtls_mpi*A :第一个乘数
const mbedtls_mpi*B:第二个乘数
const mbedtls_mpi*N:模除的大数N
mbedtls_mpi_uint mm:预计算结果-N-1(mod n)
const mbedtls_mpi*T:用来保存中间计算结果的T
函数没有提到R的取值,在后面的实现中函数采用232×n作为R。而在每一轮的运算中,我们针对的是232,之所以这样实现是因为MbedTLS中的大整数采用的字长是32位,这样在涉及到模除运算的时候,我们只需要关注大整数的最后一个字即可。
我们一步一步来看代码,首先是:
if( T->n < N->n + 1 || T->p == NULL )return( MBEDTLS_ERR_MPI_BAD_INPUT_DATA );
首先T是用来保存两个大数相乘的中间变量,T的limbs数必须足够大,当不够大或者为空的时候要返回错误。
// 将T的内容全部置为0
memset( T->p, 0, T->n * ciL );d = T->p; //d指向T->p
n = N->n; //n记录N->n
m = ( B->n < n ) ? B->n : n; //m记录N->n和B->n中较小的一个
接下来是Montgomery multiplication的核心,一个循环:
for( i = 0; i < n; i++ ){/** T = (T + u0*B + u1*N) / 2^biL*/u0 = A->p[i];u1 = ( d[0] + u0 * B->p[0] ) * mm;mpi_mul_hlp( m, B->p, d, u0 );mpi_mul_hlp( n, N->p, d, u1 );*d++ = u0; d[n + 1] = 0;}
首先,有一点我们需要明确的是,函数接收的A和B其实是:
A = a × R (mod n)
B = b × R (mod n)
这也就是说明:
0 < A , B < n - 1
这意味着A和B只是低n×Bil上的内容是有意义的,在高于n×Bil位上的内容全部是0.
因此我们在处理的时候只需要考虑A和B的低n×Bil上的内容即可。
在这个循环中,u0=A->p[i],经过循环后u0会遍历A->p中的内容。
u1计算的内容是我们在博客一中提到的:s ≡ - a0 × n-1 (mod R)
d[0] + u0 × B->p[0] = A->p[i] × B->p[0] + T->p[0]
相当于每一轮计算的低32位再加上之前计算结果中的低32位。为什么要只计算低32位呢?这就是蒙哥马利算法很巧妙的地方,因为R的取值是232,这说明一个数在模232之后只用看他的低32位内容即可,也就是说我们只用关注第一个字即可(MbedTLS中低位放在前面的字中)。
再乘上mm,我们在博客二中提到了mm ≡ -n-1 (mod R)
于是我们现在就得到了u1 ≡ - T0 × n-1 (mod R)
读者可能会注意到,为什么u1是两个字相乘的结果,应该占2个字的存储空间,为什么我们只用了一个字来存放?这是因为高位字的内容模R等于0,我们只需要关注低位字就行了。
mpi_mul_hlp( m, B->p, d, u0 );
mpi_mul_hlp( n, N->p, d, u1 );
这里计算了:
d + = u0 × B + u1 × N
在这里我们讨论一下d现在有多少个字长?
首先u0和u1都是一个字长,而B的长度小于等于n个字,N的长度等于n个字。那么d的长度应该不超过n+1个字。
图解如下:
到这里我们就只差最后一步了,这个时候的d[0]的内容已经全是0(原因可以见博客一,这里就不解释了),下面我们需要进行右移操作,将整个d向右移动一个字,在这里MbedTLS采用了一种非常聪明的办法,并不是右移数字,而是左移指针。
*d++ = u0;
最后d[n+1]=0,是为了防止内存中出现其它的数据干扰计算结果。
因为我们通过图解可以看到d的数据不可能超过d[n],d[n+1]用来保存的是
d’[n]的内容,在下一个循环开始之前,我们将这个内容清空,方便保存d’[n]。
每左移一次指针,就相当于除以232,一共左移了n次指针,一共就相当于除以了232×n,就相当于除以R。
在这个循环结束后,我们就得到了一个长度为n+1个字的d,此时的
d终=d[n]。
此时的d[0]到d[n-1]也记录了A->p的n个字的内容。
在这个循环结束之后,我们就已经计算到了:
X = A × B × R-1 (mod n)
下面我们将这(n+1)个字的值copy给A:
memcpy( A->p, d, ( n + 1 ) * ciL );
就相当于完成了:
A = A × B × R-1 (mod n)
最后我们需要检查一下A是否大于N,如果大于N的话就直接 A=A-N就行,在数学上有严格的证明可以证明 A < 2×N 。
if( mbedtls_mpi_cmp_abs( A, N ) >= 0 )mpi_sub_hlp( n, N->p, A->p );else/* prevent timing attacks */mpi_sub_hlp( n, A->p, T->p );
最后将T->p中的低n个字(也就是A->p)的内容清除是为了防止时序攻击,与算法的实现原理关系不大。
写在最后
从刚开始研究蒙哥马利到写完这篇博客,已经断断续续地过了一个星期的时间,在这个时间里我慢慢了解到了蒙哥马利算法的精巧,同时在实现的过程中也有很多小技巧需要去琢磨。在师兄们的无私帮助之下,我总算写完了这三篇博客,博客内也可能有一些错误,欢迎大家的批评指正。
MbedTLS中的Montgomery算法实现解析(三)相关推荐
- 碰撞检测GJK算法论文解析三
碰撞检测GJK算法论文解析三 再探Appendix Ⅱ 内容详解 再探The Distance Subalgorithm 内容详解 过程1 过程2 过程3 这里要先纠正上篇文章的一些错误,就是上篇文章 ...
- 容迟网络中的路由算法笔记(三)
第三章 基于地理信息的路由算法 背景:为了提高消息成功投递的可能性,一种普遍受认可的方式是采用基于洪泛的多副本策略,通过引入更多的消息副本,增大消息与目的节点相遇的机会.为了实现消息的受控洪泛,减少消 ...
- STL中算法锦集(三)
STL中算法锦集(三) 文章目录 STL中算法锦集(三) 一.< algorithm > 1.std::find_if 2.std::find_if_not 3.std::for_each ...
- 二叉树的中序遍历非递归方法(算法导论第三版12.1-3)
二叉树的中序遍历非递归方法(算法导论第三版12.1-3) 1⃣️用栈实现 template<typename T> void inorder_tree_walk_non_recursion ...
- 二叉树的遍历(算法导论第三版12.1-4)(包含先序遍历,后序遍历和中序遍历)
二叉树的遍历(算法导论第三版12.1-4) 1⃣️先序遍历 template<typename T> void preorder_tree_wald(BinaryTreeNode<T ...
- 找出有序数组X和Y中所有元素的中位数(X,Y分别含n个元素)(算法导论第三版9.3-8)
找出有序数组X和Y中所有元素的中位数(X,Y分别含n个元素) (算法导论第三版9.3-8) 时间复杂度O(lgn) int find_median_two_ordered_arrays(int *ar ...
- 确定S中最接近中位数的k个元素(算法导论第三版9.3-7)
确定S中最接近中位数的k个元素 (算法导论第三版9.3-7题) 时间复杂度O(n) vector<int> k_elements_closest_to_median(int *array, ...
- 恺撒密码是古罗马恺撒大帝用来对军事情报进行加解密的算法,它采用了替换方法对信息中的每一个英文字符循环替换为字母表序列中该字符后面的第三个字符,即,字母表的对应关系如下:
题目: 恺撒密码是古罗马恺撒大帝用来对军事情报进行加解密的算法,它采用了替换方法对信息中的每一个英文字符循环替换为字母表序列中该字符后面的第三个字符,即,字母表的对应关系如下: 原文:A B C D ...
- c语言rr算法,[判断题] 在RR、PF、MAXC/I三种算法中,RR算法的用户公平性最好
[判断题] 在RR.PF.MAXC/I三种算法中,RR算法的用户公平性最好 更多相关问题 在等差数列中,有,则此数列的前13项之和为()A.24B.39C.52D.104 已知y=xlnx,则y(10 ...
最新文章
- 服务器温度3d显示,智能问答助手、3D可视化展示,腾讯医典“黑科技”助力科普更有温度...
- matlab中用于小数取整的函数的用法
- 2字节十六进制浮点数 qt_Qt中如何实现十六进制“41A4 0000”十六进制转为浮点数20.5呢?...
- ITK:分段线性曲线的数据结构
- 51nod 1575 Gcd and Lcm
- neo4jd3的使用流程(转载)
- 【Python】Scrapy的安装与使用
- 如何给table表格的tr行加border边框(解决篇)
- 数据随机丢失情况下多传感器多速率鲁棒融合估计
- 彻底告别加解密模块代码拷贝-JCE核心Cpiher详解
- paip.提升性能---- 网站并发数的总结.txt
- matlab 填补空洞,OpenCV空洞填充算法
- opengl简单模拟行星运转
- e4e反演框架:Designing an Encoder for StyleGAN Image Manipulation
- WinDbg非常简单的调试dmp文件
- 启动RocketMQ报错:Please set the JAVA_HOME variable in your environment, We need java 64
- 每个机器学习工程师都应该知道的机器学习算法
- 京东方劝退_正式确认,华为mate40正面照来了,国行价“劝退”!
- 星座运势,周公解梦流量主微信小程序源码下载
- 软件工程_东师站_第六周作业