200字范文 > 密码学基础算法(一)基于整数的欧几里得算法和扩展欧几里得算法

密码学基础算法(一)基于整数的欧几里得算法和扩展欧几里得算法

时间：2018-06-22 13:11:31

图片来源: 随便谷歌的一个图片
图片地址: https://jason-chen-/uploads/1/0/8/5/108557741/euclidean_3_orig.jpg

密码学基础系列:

(一) 基于整数的欧几里得算法和扩展欧几里得算法(二) 扩展域 GF(2n)GF(2^n)GF(2n) 上的多项式运算(三) 扩展域 GF(2n)GF(2^n)GF(2n) 上的欧几里得算法和扩展欧几里得算法

文章目录

1. 欧几里得算法求最大公约数1. 递归版本2. 非递归版本2. 扩展欧几里得算法解线性方程 ax+by=d=gcd(a,b)ax + by = d = gcd(a, b)ax+by=d=gcd(a,b)1. 递归版本2. 非递归版本3. 求整数的乘法逆元4. 完整代码和测试4.1 源码4.2 基于gtest的单元测试5. 其它

求两个整数的最大公约数，以及求整数的乘法逆元是密码学中最基本的算法。

1. 欧几里得算法求最大公约数

欧几里得算法(Euclidean Algorithm)，即辗转相除算法，是求两个整数最大公约数(GCD, Greatest Common Devisor)最常用和高效的办法。

关于辗转相除的原理和证明这里不再赘述，推荐参考:

《数论概论》原书第3版，第5章，整除性与最大公因数《密码编码学与网络安全》第七版，第2.2节，欧几里得算法《深入浅出密码学》，第6.3.1节，欧几里得算法

以下是源码实现的递归和非递归版本。

1. 递归版本

/*** @description: 使用欧几里得算法(辗转相除)求最大公约数的递归版本* @param {int} a 整数a* @param {int} b 整数b* @return {*} 返回整数 (a, b) 的最大公约数*/int int_gcd(int a, int b){if (b == 0){return a;}else{return int_gcd(b, a % b);}}

2. 非递归版本

/*** @description: 使用欧几里得算法(辗转相除)求最大公约数的非递归版本* @param {int} a 整数a* @param {int} b 整数b* @return {*} 返回整数 (a, b) 的最大公约数*/int int_gcd(int a, int b){int t;while (b != 0){t = a % b;a = b;b = t;}return a;}

为了和后面基于多项式的欧几里得算法求最大公因式区分，这里将函数命名为int_gcd，前缀 int 表示 integer

2. 扩展欧几里得算法解线性方程 ax+by=d=gcd(a,b)ax + by = d = gcd(a, b)ax+by=d=gcd(a,b)

扩展欧几里得算法(Extend Euclidean Algorithm)的最大用途就是用于计算整数的乘法逆元。

关于扩展欧几里得算法(EEA)的原理和和推导这里不再赘述，推荐参考:

《数论概论》原书第3版，第6章，线性方程与最大公因数《密码编码学与网络安全》第七版，第2.3.6节，扩展的欧几里得算法《深入浅出密码学》，第6.3.2节，扩展的欧几里得算法

以下的源码基于扩展欧几里得算法实现求解线性方程:

ax+by=d=gcd(a,b),公式2.1ax + by = d = gcd(a, b) \text{, 公式2.1} ax+by=d=gcd(a,b),公式2.1

以下是源码实现的递归和非递归版本。

1. 递归版本

/*** @description: 使用扩展欧几里得算法(EEA)计算等式 a x ia + b x ib = gcd(a, b) 中的系数 ia, ib 和 gcd(a, b), 递归版本, 推导和证明参考: /lincifer/article/details/49391175* @param {int} a 整数 a* @param {int} b 整数 b* @param {int} *ia 整数 a 的系数 ia* @param {int} *ib 整数 b 的系数 ib* @return {*} 返回整数 (a, b) 的最大公约数 gcd(a, b)*/int int_gcd_ex(int a, int b, int *ia, int *ib){int x, y;int q, r;if (b == 0){*ia = 1;*ib = 0;return a;}r = int_gcd_ex(b, a % b, &x, &y);*ia = y;*ib = x - a / b * y;return r;}

一般的书讲述扩展欧几里得算法(Extend Euclidean Algorithm)都是基于递推原理进行的。
关于递归(非递推)版本的推导和证明，请参考参考:
/lincifer/article/details/49391175

2. 非递归版本

/*** @description: 使用扩展欧几里得算法: Extend Euclidean Algorithm (EEA) 计算等式 a x ia + b x ib = gcd(a, b) 中的 ia, ib 和 gcd(a, b)* @param {int} a 整数 a* @param {int} b 整数 b* @param {int} *ia 整数 a 的系数 ia* @param {int} *ib 整数 b 的系数 ib* @return {*} 返回整数 (a, b) 的最大公约数 gcd(a, b)*/int int_gcd_ex(int a, int b, int *ia, int *ib){int x, y, x0, y0, x1, y1;int q, r;/* 消除警告: "warning: ‘x’/‘y’ may be used uninitialized in this function" */x = y = 0;/* 初始化最初两项系数 */x0 = 1; y0 = 0;x1 = 0; y1 = 1;q = a / b;r = a % b;while (r != 0){/* 计算当前项 x/y */x = x0 - q * x1;y = y0 - q * y1;/* 依次保存前两项到 x0/y0, x1/y1 */x0 = x1; x1 = x;y0 = y1; y1 = y;a = b;b = r; /* 前一次的余数 */q = a / b;r = a % b;}*ia = x;*ib = y;return b;}

函数int_gcd_ex的返回值和int_gcd一样，都返回其参数a和b的最大公约数，不同的是，除了返回最大公约数之外，还通过指针参数返回线性方程 2.1 中的x和y值。

二者的函数原型如下:

int int_gcd(int a, int b);int int_gcd_ex(int a, int b, int *x, int *y);

3. 求整数的乘法逆元

在第2节的公式 2.1 中，当 gcd(a,b)=1gcd(a, b) = 1gcd(a,b)=1 时，有:

ax+by=1=gcd(a,b)ax + by = 1 = gcd(a, b) ax+by=1=gcd(a,b)

等式两边同时 modbmod\ bmodb，有: axmodb+bymodb≡1modbax\ mod\ b + by\ mod\ b \equiv 1\ mod\ baxmodb+bymodb≡1modb，而 bymodb=0by\ mod\ b = 0bymodb=0。

所以:

ax≡1modbax \equiv 1\ mod\ bax≡1modb

因此，a 和 x 在模b 的前提下互为逆元。

第2节的函数int_gcd_ex在求整数 a 和整数 b 最大公约数的同事，可以顺带求出 a 和 b 的乘法逆元。

以下代码复用int_gcd_ex，求整数 a 基于整数 b 的乘法逆元:

/*** @description: 基于扩展欧几里得算法计算整数 a 关于整数 b 的乘法逆元* @param {int} a 整数 a* @param {int} b 整数 b* @return {*} 返回整数 a 关于整数 b 的最小正整数乘法逆元, 乘法逆元不存在(gcd(a, b) != 1)则返回 0;*/int int_inv(int a, int b){int res, ia, ib;res = int_gcd_ex(a, b, &ia, &ib);/* gcd(a, b) != 1, 没有乘法逆元 */if (res != 1){return 0;}/* 调整小于 0 的情况 */if (ia < 0){ia += b;}return ia;}

这里的函数名是integer inverse的简写，如果整数 a 和整数 b 的最大公约数为1(互素)，则返回整数 a 关于整数 b 的最小正整数逆元，如果不能满足这个条件，则函数返回 0。

上面的代码中，在求乘法逆元时，通过扩展欧几里得算法int_gcd_ex同时得到了整数 a 和整数 b 的逆元 ia 和 ib，但整数 b 的逆元 ib 并没有被使用，所以可以在int_inv函数中将int_gcd_ex的代码展开，去掉所有整数 b 逆元计算相关代码，得到的新的代码如下:

/*** @description: 基于扩展欧几里得算法计算整数 a 关于整数 b 的乘法逆元(展开 int_gcd_ex 后的优化版本)* @param {int} a 整数 a* @param {int} b 整数 b* @return {*} 返回整数 a 关于整数 b 的最小正整数乘法逆元, 乘法逆元不存在(gcd(a, b) != 1)则返回 0;*/int int_inv(int a, int b){int x, x0, x1;int q, r, t;/* 消除警告: "warning: ‘x’/‘y’ may be used uninitialized in this function" */x = 0;/* 初始化最初两项系数 */x0 = 1;x1 = 0;/* 临时保存 b */t = b;q = a / b;r = a % b;while (r != 0){/* 计算当前项 x */x = x0 - q * x1;/* 依次保存前两项到 x0, x1 */x0 = x1; x1 = x;a = b;b = r; /* 前一次的余数 */q = a / b;r = a % b;}/* gcd(a, b) != 1, 没有乘法逆元 */if (b != 1){return 0;}/* 调整小于 0 的情况 */if (x < 0){x += t;}return x;}

优化的内容主要是减少整数 b 计算乘法逆元所使用相关变量 y, y0, y1 的赋值，减少这 3 个变量所占用的空间。实际加解密中，求逆元的代码并不会非常频繁被调用，因此这个优化带来的效率提升我认为比较有限。

相比之下，我还是更喜欢调用int_gcd_ex求乘法逆元的版本，因为这个版本代码更简洁，思路更清晰。

4. 完整代码和测试

4.1 源码

头文件gcd.h:

/** @ file: gcd.h* @ description: 整数和扩展域上多项式的最大公约数(gcd), 乘法逆元(multiple reverse)的计算* @author: Yongqiang Gu* @ blog: /guyongqiangx*/#ifndef __GCD_H__#define __GCD_H__#ifdef __cplusplusextern "C"{#endif/*** @description: 使用欧几里得算法(辗转相除)求最大公约数* @param {int} a 整数a* @param {int} b 整数b* @return {*} 返回整数 (a, b) 的最大公约数*/int int_gcd(int a, int b);/*** @description: 使用扩展欧几里得算法: Extend Euclidean Algorithm (EEA) 计算等式 a x ia + b x ib = gcd(a, b) 中的 ia, ib 和 gcd(a, b)* @param {int} a 整数 a* @param {int} b 整数 b* @param {int} *ia 整数 a 的系数 ia* @param {int} *ib 整数 b 的系数 ib* @return {*} 返回整数 (a, b) 的最大公约数 gcd(a, b)*/int int_gcd_ex(int a, int b, int *ia, int *ib);/*** @description: 基于扩展欧几里得算法计算整数 a 关于整数 b 的乘法逆元* @param {int} a 整数 a* @param {int} b 整数 b* @return {*} 返回整数 a 关于整数 b 的最小正整数乘法逆元, 乘法逆元不存在(gcd(a, b) != 1)则返回 0;*/int int_inv(int a, int b);#ifdef __cplusplus}#endif#endif

源码文件int_gcd.c:

/** @ file: int_gcd.c* @ description: 整数的最大公约数(欧几里得算法 Euclidean Algorithm)，和扩展欧几里得算法(Extend Euclidean Algorithm)求乘法逆元的实现* @author: Yongqiang Gu* @ blog: /guyongqiangx*/#include "gcd.h"///**// * @description: 使用欧几里得算法(辗转相除)求最大公约数的递归版本// * @param {int} a 整数a// * @param {int} b 整数b// * @return {*} 返回整数 (a, b) 的最大公约数// *///int int_gcd(int a, int b)//{// if (b == 0)// {// return a;// }// else// {// return int_gcd(b, a % b);// }//}/*** @description: 使用欧几里得算法(辗转相除)求最大公约数的非递归版本* @param {int} a 整数a* @param {int} b 整数b* @return {*} 返回整数 (a, b) 的最大公约数*/int int_gcd(int a, int b){int t;while (b != 0){t = a % b;a = b;b = t;}return a;}///**// * @description: 使用扩展欧几里得算法(EEA)计算等式 a x ia + b x ib = gcd(a, b) 中的系数 ia, ib 和 gcd(a, b), 递归版本, 推导和证明参考: /lincifer/article/details/49391175// * @param {int} a 整数 a// * @param {int} b 整数 b// * @param {int} *ia 整数 a 的系数 ia// * @param {int} *ib 整数 b 的系数 ib// * @return {*} 返回整数 (a, b) 的最大公约数 gcd(a, b)// *///int int_gcd_ex(int a, int b, int *ia, int *ib)//{// int x, y;// int q, r;//// if (b == 0)// {// *ia = 1;// *ib = 0;//// return a;// }//// r = int_gcd_ex(b, a % b, &x, &y);//// *ia = y;// *ib = x - a / b * y;//// return r;//}/*** @description: 使用扩展欧几里得算法(EEA) 计算等式 a x ia + b x ib = gcd(a, b) 中的系数 ia, ib 和 gcd(a, b), 非递归版本* @param {int} a 整数 a* @param {int} b 整数 b* @param {int} *ia 整数 a 的系数 ia* @param {int} *ib 整数 b 的系数 ib* @return {*} 返回整数 (a, b) 的最大公约数 gcd(a, b)*/int int_gcd_ex(int a, int b, int *ia, int *ib){int x, y, x0, y0, x1, y1;int q, r;/* 消除警告: "warning: ‘x’/‘y’ may be used uninitialized in this function" */x = y = 0;/* 初始化最初两项系数 */x0 = 1; y0 = 0;x1 = 0; y1 = 1;q = a / b;r = a % b;while (r != 0){/* 计算当前项 x/y */x = x0 - q * x1;y = y0 - q * y1;/* 依次保存前两项到 x0/y0, x1/y1 */x0 = x1; x1 = x;y0 = y1; y1 = y;a = b;b = r; /* 前一次的余数 */q = a / b;r = a % b;}*ia = x;*ib = y;return b;}///**// * @description: 基于扩展欧几里得算法计算整数 a 关于整数 b 的乘法逆元// * @param {int} a 整数 a// * @param {int} b 整数 b// * @return {*} 返回整数 a 关于整数 b 的最小正整数乘法逆元, 乘法逆元不存在(gcd(a, b) != 1)则返回 0;// *///int int_inv(int a, int b)//{// int res, ia, ib;//// res = int_gcd_ex(a, b, &ia, &ib);//// /* gcd(a, b) != 1, 没有乘法逆元 */// if (res != 1)// {// return 0;// }//// /* 调整小于 0 的情况 */// if (ia < 0)// {// ia += b;// }//// return ia;//}/*** @description: 基于扩展欧几里得算法计算整数 a 关于整数 b 的乘法逆元(展开 int_gcd_ex 后的优化版本)* @param {int} a 整数 a* @param {int} b 整数 b* @return {*} 返回整数 a 关于整数 b 的最小正整数乘法逆元, 乘法逆元不存在(gcd(a, b) != 1)则返回 0;*/int int_inv(int a, int b){int x, x0, x1;int q, r, t;/* 消除警告: "warning: ‘x’ may be used uninitialized in this function" */x = 0;/* 初始化最初两项系数 */x0 = 1;x1 = 0;/* 临时保存 b */t = b;q = a / b;r = a % b;while (r != 0){/* 计算当前项 x */x = x0 - q * x1;/* 依次保存前两项到 x0, x1 */x0 = x1; x1 = x;a = b;b = r; /* 前一次的余数 */q = a / b;r = a % b;}/* gcd(a, b) != 1, 没有乘法逆元 */if (b != 1){return 0;}/* 调整小于 0 的情况 */if (x < 0){x += t;}return x;}

4.2 基于gtest的单元测试

测试代码

基于gtest的单元测试:

// g++ int_gcd.c -o int_gcd -Igtest/include -Lgtest/lib -lgtest_main -lgtest -lpthread#include "gcd.h"#include <gtest/gtest.h>TEST(IntegerGCDTest, CompositeTest){/* gcd(60, 0) = 60, gcd(0, 20) = 20 */EXPECT_EQ(60, int_gcd(60, 0));EXPECT_EQ(20, int_gcd( 0, 20));EXPECT_EQ(12, int_gcd(60, 24));EXPECT_EQ(11, int_gcd(55, 22));EXPECT_EQ( 3, int_gcd(33, 18));EXPECT_EQ( 3, int_gcd(27, 21));EXPECT_EQ( 6, int_gcd(84, 30));EXPECT_EQ(10, int_gcd(50, 30));EXPECT_EQ( 7, int_gcd(973, 301));EXPECT_EQ(77, int_gcd(7469, 2464));EXPECT_EQ(34, int_gcd(24140, 16762));EXPECT_EQ(35, int_gcd(4655, 12075));EXPECT_EQ(1078, int_gcd(1160718174, 316258250));}TEST(IntegerGCDTest, PrimeTest){EXPECT_EQ(1, int_gcd(67, 12));EXPECT_EQ(1, int_gcd(29, 100));EXPECT_EQ(1, int_gcd(1759, 550));EXPECT_EQ(1, int_gcd(2689, 4001));}TEST(IntegerExtEuclideanTest, GCDTest){int ia, ib;EXPECT_EQ(3, int_gcd_ex(33, 18, &ia, &ib));EXPECT_EQ(3, int_gcd_ex(18, 33, &ia, &ib));EXPECT_EQ(1, int_gcd_ex(100, 29, &ia, &ib));EXPECT_EQ(1, int_gcd_ex(29, 100, &ia, &ib));EXPECT_EQ(10, int_gcd_ex(50, 30, &ia, &ib));EXPECT_EQ(1, int_gcd_ex(1759, 550, &ia, &ib));}TEST(IntegerExtEuclideanTest, InverseTest){int res, ia, ib;/* 198 x (-11) + 243 x 9 = 9 = gcd(198, 243) */res = int_gcd_ex(198, 243, &ia, &ib);EXPECT_EQ( 9, res);EXPECT_EQ(-11, ia);EXPECT_EQ( 9, ib);/* 1819 x 71 + 3587 x (-36) = 17 = gcd(1819, 3587) */res = int_gcd_ex(1819, 3587, &ia, &ib);EXPECT_EQ( 17, res);EXPECT_EQ( 71, ia);EXPECT_EQ(-36, ib);res = int_gcd_ex(100, 29, &ia, &ib);EXPECT_EQ( 1, res);EXPECT_EQ( 9, ia);EXPECT_EQ(-31, ib);res = int_gcd_ex(12, 67, &ia, &ib);EXPECT_EQ( 1, res);EXPECT_EQ(28, ia);EXPECT_EQ(-5, ib);res = int_gcd_ex(973, 301, &ia, &ib);EXPECT_EQ( 7, res);EXPECT_EQ( 13, ia);EXPECT_EQ(-42, ib);res = int_gcd_ex(1234, 4321, &ia, &ib);EXPECT_EQ( 1, res);EXPECT_EQ(-1082, ia);EXPECT_EQ( 309, ib);res = int_gcd_ex(24140, 40902, &ia, &ib);EXPECT_EQ( 34, res);EXPECT_EQ(-571, ia);EXPECT_EQ(337, ib);res = int_gcd_ex(1759, 550, &ia, &ib);EXPECT_EQ( 1, res);EXPECT_EQ(-111, ia);EXPECT_EQ( 355, ib);res = int_gcd_ex(550, 1769, &ia, &ib);EXPECT_EQ( 1, res);EXPECT_EQ( 550, ia);EXPECT_EQ(-171, ib);}TEST(IntegerInverseTest, InverseTest){/* gcd(33, 18) = 3, no inverse */EXPECT_EQ(0, int_inv(33, 18));/* gcd(30, 50) = 10, no inverse */EXPECT_EQ(0, int_inv(30, 50));/* gcd(12, 67) = 1, 12 x 28 mod 67 = 1 */EXPECT_EQ(28, int_inv(12, 67));/* gcd(1759, 550) = 1, 1759 x 439 mod 550 = 1 */EXPECT_EQ(439, int_inv(1759, 550));EXPECT_EQ(355, int_inv(550, 1759));EXPECT_EQ(3239, int_inv(1234, 4321));EXPECT_EQ(0, int_inv(24140, 40902));EXPECT_EQ(550, int_inv(550, 1769));EXPECT_EQ(69, int_inv(29, 100));EXPECT_EQ(9, int_inv(100, 29));}

测试结果

[==========] Running 5 tests from 3 test suites.[----------] Global test environment set-up.[----------] 2 tests from IntegerGCDTest[ RUN] positeTest[ OK ] positeTest (0 ms)[ RUN] IntegerGCDTest.PrimeTest[ OK ] IntegerGCDTest.PrimeTest (0 ms)[----------] 2 tests from IntegerGCDTest (0 ms total)[----------] 2 tests from IntegerExtEuclideanTest[ RUN] IntegerExtEuclideanTest.GCDTest[ OK ] IntegerExtEuclideanTest.GCDTest (0 ms)[ RUN] IntegerExtEuclideanTest.InverseTest[ OK ] IntegerExtEuclideanTest.InverseTest (0 ms)[----------] 2 tests from IntegerExtEuclideanTest (0 ms total)[----------] 1 test from IntegerInverseTest[ RUN] IntegerInverseTest.InverseTest[ OK ] IntegerInverseTest.InverseTest (0 ms)[----------] 1 test from IntegerInverseTest (0 ms total)[----------] Global test environment tear-down[==========] 5 tests from 3 test suites ran. (0 ms total)[ PASSED ] 5 tests.