Miracle密码算法开源库（零） 项目综述

2021SC@SDUSC 山东大学软件学院软件工程应用与实践

一．MIRACL简介

MIRACL(Multiprecision Integer and RationalArithmetic C/c++ Library)是一套由Shamus Software Ltd.所开发的一套关于大数运算函数库，用来设计与大数运算相关的密码学之应用，包含了RSA 公开密码学、Diffie-Hellman密钥交换(Key Exchange)、AES、DSA数字签名，还包含了较新的椭圆曲线密码学(Elliptic CurveCryptography)等等。运算速度快，并提供源代码。MIARCL是当前使用比较广泛的基于公钥加密算法保护实现的大数库之一，据说要使用该库用于商业软件，需要交纳一笔昂贵的授权费——1000$。

二、下载

下载地址：https://github.com/miracl/MIRACL

三、 编译miracl

将从github网站上下载的安装包解压，如下图：

2.新建文件夹miracle，将解压出来的所有文件复制到miracl文件夹中，如下图：
3.打开vs 2017，新建 Visual C++–空项目CompileMiracl，把miracl文件夹里的所有文件复制到CompileMiracl的工程目录，和工程源文件放在一个文件夹，如下图：
4.右键单击项目–添加–现有项，如下图：

5.添加如下头文件和源文件：

miracl.h
mirdef.h

mraes.c
mralloc.c
mrarth0.c
mrarth1.c
mrarth2.c
mrarth3.c
mrbits.c
mrbrick.c
mrbuild.c
mrcore.c
mrcrt.c
mrcurve.c
mrdouble.c
mrebrick.c
mrec2m.c
mrgf2m.c
mrfast.c
mrflash.c
mrflsh1.c
mrflsh2.c
mrflsh3.c
mrflsh4.c
mrfrnd.c
mrgcd.c
mrgcm.c
mrio1.c
mrio2.c
mrjack.c
mrlucas.c
mrmonty.c
mrmuldv.c
mrpi.c
mrpower.c
mrprime.c
mrrand.c
mrround.c
mrscrt.c
mrshs.c
mrshs256.c
mrshs512.c
mrsmall.c
mrsroot.c
mrstrong.c
mrxgcd.c
mrzzn2.c
mrzzn2b.c
mrzzn3.c
mrzzn4.c
mrecn2.c

6.右键单击项目–属性，如下图：

7.更改目标文件扩展名和配置类型，如下图：

8.右键单击项目–生成/重新生成，如下图：

编译完成
想要使用miracle库的函数时，进入之前新建的CompileMiracl项目的debug以及源代码文件夹，拷贝 “mircal.h”、”mirdef.h”、”CompileMiracl.lib”、”CompileMiracl.pdb”到新建空项目的源代码路径，并且把”CompileMiracl.lib”、”CompileMiracl.pdb”更名为 “miracl.lib” , “miracl.pdb”即可。

四、miracl底层算法

1.大整数的乘法
1.1 基线乘法 （Baseline Multiplication）
基线乘法是 Bignum Math 一书中的一个术语，这里直接拿来用。实际上基线乘法就是通常所说的笔算乘法。先来看看计算 123 * 456 用笔算乘法应怎么做：

1
2
3
x
4
5
6
-----------------------------------------------
7
3
8
6
1
5
4
9
2
------------------------------------------------
5
6
0
8
8

在笔算算法中，第二个数的每一位分别与第一个数的每一个位相乘，并且将进位传递到下一位，比如 3 * 6 = 18 保留本位 8，进位 1 传递到下一位。按照笔算的原理，很容易得到基线乘法的算法思路：
设大整数 x 和 y 分别有 x->used 和 y->used 个数位，进位 c，相乘结果放在 z 中，算法的步骤是：

1. 对于下标 i，从 0 到 x->used + y->used - 1 循环，将 z 中每一个数位清 0。其实就是把 z 设为 0。
2. 对于下标 i，从 0 到 y->used - 1 循环，执行如下几个操作：
   2.1 c = 0
   2.2 对于下标 j，从 0 到 x->used - 1循环计算：（注：i，j 都是从 0 开始的，r 是一个双精度变量）
   r = z(i + j) + x(j) * y(i) + c。其中 z(i + j)，x(j)，y(i) 分别表示 z 的第 i + j + 1个数位，x 的第 j + 1 个 数位，y 的第 i +1 个数位。
   z(i + j) = r & BN_MASK //取 r 的低半部分作为 z(i + j) 的本位。
   c = r >> biL //取 r 的高半部分作为进位 c。
   2.3 z(i + x->used) = c
3. 返回结果 z。

2.1 Comba 乘法 （Comba Multiplication）原理
Comba 乘法以（在密码学方面）不太出名的 Paul G. Comba 得名。上面的笔算乘法，虽然比较简单， 但是有个很大的问题：在 O(n^2) 的复杂度上进行计算和向上传递进位，看看前面的那个竖式，每计算一次单精度乘法都要计算和传递进位，这样的话就使得嵌套循环的顺序性很强，难以并行展开和实现。Comba 乘法则无需进行嵌套进位来计算乘法，所以虽然其时间复杂度和基线乘法一样，但是速度会快很多。还是以计算 123 * 456 为例：
1 2 3

x
4
5
6
-----------------------------------------------
6
12
18
5
10
15
4
8
12
------------------------------------------------
4
13
28
27
18
4
13
28
28
8
4
13
30
8
4
16
0
5
6
0
5
------------------------------------------------------
5
6
0
8
8

和普通的笔算乘法很类似，只是每一次单精度乘法只是单纯计算乘法，不计算进位，进位留到每一列累加后进行。所以原来需要 n * n 次进位，现在最多只需要 2n 次即可。
以上就是 Comba 乘法的原理，不过这里有个比较严重的问题：如何保证累加后结果不溢出。上面的例子，假设单精度数 1 位数，双精度是两位数，那万一累加后的结果超过两位数则么办？那没办法，只能用三精度变量了。在大整数算法中，单精度能表示的最大整数是 2^n - 1（n 是单精度变量的比特数），用三个单精度变量 c2，c1，c0 连在一起作为一个三精度变量（高位在左，低位在右），则 c2 || c1 || c0 能表示的最大整数是 2^(3n) - 1，最多能存放 (2^(3n) - 1) / ((2^n - 1)^2) 个单精度乘积结果。当 n = 32 时，能够存放 4294967298 个单精度乘积结果；当 n = 64 时，能够存放约 1.845 * 10^19 个单精度乘积结果，而我一开始规定 bignum 不能超过 25600 个数位，这样使用三精度变量就可以保证累加结果不会溢出了。
有了上面的铺垫，下面就把 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。
   BN_MIN 是宏定义：#define BN_MIN(x, y) (((x) < (y)) ? (x) : (y))
   以上就是 Comba 乘法的思路，先计算每一列，然后求和累加，将累加结果求余数得到本位，进位传递到下一列。
   第三步中，分别计算 tx， ty 和 k 的值，用于控制列的计算，光这么说比较抽象，所以还是举个例子比较直观。假设要计算 12345 * 678。

    Index i:
   6
   5
   4
   3
   2
   1
   0
   -----------------------------------------------------------------------
   1
   2
   3
   4
   5
   x
   6
   7
   8
   -----------------------------------------------------------------------
   8
   16
   24
   32
   40
   7
   14
   21
   28
   35
   6
   12
   18
   24
   30
   2
   4
   6
   8
   7
   4
   0
   ------------------------------------------------------------------------
   8
   3
   6
   9
   9
   1
   0
   

过程：a->used = 5， b->used = 3
i = 0，ty = 0，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 1，k = 1，循环计算第一列：5 * 8 = 40
i = 1，ty = 1，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 2，k = 2，循环计算第二列：5 * 7 = 35，4 * 8 = 32
i = 2，ty = 2，tx = i - ty = 0，a->used - tx = 5，ty + 1 = 3，k = 3，循环计算第三列：5 * 6 = 30，4 * 7 = 28，3 * 8 = 24
i = 3，ty = 2，tx = i - ty = 1，a->used - tx = 4，ty + 1 = 3，k = 3，循环计算第四列：4 * 6 = 24，3 * 7 = 21，2 * 8 = 16
i = 4，ty = 2，tx = i - ty = 2，a->used - tx = 3，ty + 1 = 3，k = 3，循环计算第五列：3 * 6 = 18，2 * 7 = 14，1 * 8 = 8
i = 5，ty = 2，tx = i - ty = 3，a->used - tx = 2，ty + 1 = 3，k = 2，循环计算第六列：2 * 6 = 12，1 * 7 = 7
i = 6，ty = 2，tx = i - ty = 4，a->used - tx = 1，ty + 1 = 3，k = 1，循环计算第七列：1 * 6 = 6
可以看出，Comba 乘法每一次都是计算一列，所以如果想进行并行计算的话，那会方便很多：同时计算所有列并累加，最后进行一次进位传递。
前面提到基线乘法和 Comba 乘法的时间复杂度是一样的，这是对于单精度乘法而言，Comba 乘法之所以要比基线乘法快，主要是因为减少了进位传递的次数，所以减少了加法的计算量。
2.模乘
先定义：若N，R为正整数，且gcd（N，R）=1， R>N。
设定：R‘是R的模N逆元， （-N‘）是N的模R逆元。
对于整数a，b；模乘就是M：“a*b mod N”
模乘基本运算法则：
(a + b) % p = (a % p + b % p) % p （1）

(a – b) % p = (a % p – b % p) % p （2）

(a * b) % p = (a % p * b % p) % p （3）

(a^b) % p = ((a % p)^b) % p （4）

结合律：

((a+b) % p + c) % p = (a + (b+c) % p) % p （5）

((ab) % p * c)% p = (a * (bc) % p) % p （6）

交换律：

(a + b) % p = (b+a) % p （7）

(a * b) % p = (b * a) % p （8）

分配律：

((a +b)% p * c) % p = ((a * c) % p + (b * c) % p) % p （9）

重要定理：

若a≡b (% p)，则对于任意的c，都有(a + c) ≡ (b + c) (%p)；（10）

若a≡b (% p)，则对于任意的c，都有(a * c) ≡ (b * c) (%p)；（11）

若a≡b (% p)，c≡d (% p)，则 (a + c) ≡ (b + d) (%p)，(a – c) ≡ (b – d) (%p)，

(a * c) ≡ (b * d) (%p)，(a / c) ≡ (b / d) (%p)；

3.Montgomery模乘
Montgomery乘法是公钥算法实现中的一个核心算法，其主要作用是为模乘运算加速。

在公钥算法实现中，通常需要计算a mod M、a*b mod M、a^b mod M等，一般看见mod M，最直接想到的当然是除法，可是除法运算慢且实现难，于是，就有人发明了一种不需要计算除法的计算模的方法，这就是所谓的Montgomery乘法。

Montgomery乘法的数学表达是ABR^(-1) mod M，其中，A、B是与M同位长的大数，R = 2^{bitlen(bitlen指M位长)，R}(-1)是指R相对于M的模逆，即R^{(-1)是满足如下条件的数：R*R}(-1) mod M = 1；这个条件成立的充要条件是R与M互素，这一点只需要M为奇数即可，所以Montgomery乘法通常适用于对奇数求模。
蒙哥马利模乘流程：

4.扩展欧几里得算法
为了介绍扩展欧几里得，我们先介绍一下贝祖定理：

即如果a、b是整数，那么一定存在整数x、y使得ax+by=gcd(a,b)。

换句话说，如果ax+by=m有解，那么m一定是gcd(a,b)的若干倍。（可以来判断一个这样的式子有没有解）

有一个直接的应用就是 如果ax+by=1有解，那么gcd(a,b)=1；

要求出这个最大公因数gcd(a,b)，我们最容易想到的就是古老悠久而又相当强大的辗转相除法：

1. int gcd(int a,int b)
2. {
3. return b==0?a:gcd(b,a%b);
4. }

但是，对于上面的式子ax+by=m来说，我们并不仅仅想要知道有没有解，而是想要知道在有解的情况下这个解到底是多少。

所以，扩展欧几里得

当到达递归边界的时候，b==0，a=gcd(a,b) 这时可以观察出来这个式子的一个解：a*1+b*0=gcd(a,b)，x=1,y=0，注意这时的a和b已经不是最开始的那个a和b了，所以我们如果想要求出解x和y，就要回到最开始的模样。
初步想法：由于是递归的算法，如果我们知道了这一层和上一层的关系，一层一层推下去，就可以推到最开始的。类似数学上的数学归纳法。
假设当前我们在求的时a和b的最大公约数，而我们已经求出了下一个状态：b和a%b的最大公因数，并且求出了一组x1和y1使得 b*x1+(a%b)*y1=gcd

（注意在递归算法中，永远都是先得到下面一个状态的值）

这时我们可以试着去寻找这两个相邻状态的关系：

首先我们知道：a%b=a-(a/b)*b；带入：

b*x1 + (a-(a/b)*b)*y1

= bx1 + ay1 – (a/b)by1

= ay1 + b(x1 – a/by1) = gcd 发现 x = y1 , y = x1 – a/by1

这样我们就得到了每两个相邻状态的x和y的转化，就可以在求gcd的同时对x和y进行求值了hiahia
5.模幂算法
模幂是一种对模进行的冪运算，在计算机科学，尤其是公开密钥加密方面有一定用途。 模幂运算是指求整数b的e次方b 被正整数m所除得到的余数c的过程，可用数学符号表示为c = b mod m。由c的定义可得0 ≤ c < m。 例如，给定b = 5，e = 3和m = 13，5 = 125被13除得的余数c = 8。 指数e为负数时可使用扩展欧几里得算法找到b模除m的模逆元d来执行模幂运算，即： c = b mod m = d mod m，其中 e < 0且b ⋅ d ≡ 1 (mod m)。 即使在整数很大的情况下，上述模幂运算其实也是易于执行的。然而，计算模的离散对数（即在已知b，c和m时求出指数e）则较为困难。这种类似單向函數的表现使模幂运算可用于加密算法。

五、关键代码

通过阅读源代码中的说明文档，可以了解到源代码分为示例程序和函数，其中的示例程序都是对不同函数的应用，如下图：
示例程序

函数说明

仔细观察这些函数所属的模块，可以发现这些模块都是上文中编译miracl时拷贝到CompileMiracl工程目录中的模块。所以，miracl开源库的关键代码就在这些模块之中，如下图：

六、团队分工

项目成员共两人，由于一个模块中有多个函数，且每个人不可以分析相同代码，我们决定将上述关键代码模块平均分，一人分析一半。

最后

以上就是温暖自行车为你收集整理的2021-09-30Miracle密码算法开源库（零） 项目综述一．MIRACL简介的全部内容，希望文章能够帮你解决2021-09-30Miracle密码算法开源库（零） 项目综述一．MIRACL简介所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。