大数模运算相关资料

 RSA算法中用到的大数运算

C. 大数的运算
1. 大数的运算原理
RSA算法依赖于大数的运算，目前主流RSA算法都建立在512位到1024位的大数运算之

上，所以我们首先需要掌握大数（比如1024位）的运算原理。
大多数的编译器只能支持到32位（或64位）的整数运算，即我们在运算中所使用的

整数必须小于等于32位，即0xFFFFFFFF，这远远达不到RSA的需要，于是需要专门建

立大数运算库，来解决这一问题。
最简单的办法是将大数当作字符串进行处理，也就是将大数用10进制字符数组进行

表示，然后模拟人们手工进行“竖式计算”的过程编写其加减乘除函数。但是这样

做效率很低。当然其优点是算法符合人们的日常习惯，易于理解。
另一种思路是将大数当作一个二进制流进行处理，使用各种移位和逻辑操作来进行

加减乘除运算，但是这样做代码设计非常复杂，可读性很低，难以理解也难以调试

。
这里我们采用了一种介于两者之间的思路：将大数看作一个N进制数组，对于目前的

32位系统而言，N可以取2的32次方，即0x100000000，假如将一个1024位的大数转化

成0x10000000进制，它就变成了32位，而每一位的取值范围是0－0xFFFFFFFF。我们

正好可以用一个无符号长整数来表示这一数值。所以1024位的大数就是一个有32个

元素的unsigned long数组。而且0x100000000进制的数组排列与2进制流对于计算机

来说，实际上是一回事，但是我们完全可以针对unsigned long数组进行“竖式计算

”，而循环规模被降低到了32次之内，并且算法很容易理解。
但考虑到乘法和除法，都要进行扩展才能进行快速的计算（如果把除法变减法而不

扩展，速度将慢的无法忍受）。所以我们将N取2的16次方的，即0xFFFF。每一位用

unsigned short来表示，当进行乘除运算时，将short扩展成long，这是编译器所支

持的，所以运算起来，比较快。
参考资料：http://hi.baidu.com/gaojinshan/blog/item/abf80613cc713b26dd540110.html

RSA算法的介绍


A. 加密解密
1. 密钥的产生
1) 找出两个相异的大素数P和Q，令N＝P×Q，M＝（P－1）（Q－1）。
2) 找出与M互素的大数E，用欧氏算法计算出大数D，使D×E≡1 MOD M。
3) 丢弃P和Q，公开E，D和N。E和N即加密密钥，D和N即解密密钥。
2. 加密的步骤
1) 计算N的有效位数tn（以字节数计），将最高位的零忽略掉，令tn1＝tn－1

。比如N＝0x012A05，其有效位数tn＝5，tn1＝4。
2) 将明文数据A分割成tn1位（以字节数计）的块，每块看成一个大数，块数

记为bn。从而，保证了每块都小于N。
3) 对A的每一块Ai进行Bi＝Ai^E MOD N运算。Bi就是密文数据的一块，将所有

密文块合并起来，就得到了密文数据B。
3. 解密的步骤
1) 同加密的第一步。
2) 将密文数据B分割成tn位（以字节数计）的块，每块看成一个大数，块数记

为bn。
3) 对B的每一块Bi进行Ci＝Bi^D MOD N运算。Ci就是密文数据的一块，将所有

密文块合并起来，就得到了密文数据C。
4. 定理及证明
<定理>
费马小定理：P是任意一个素数，Q是任意一个整数，则P^Q≡P MOD Q。
换句话说，如果P和Q互质，则P^（Q－1）≡1 MOD Q。
<证明>
运用一些基本的群论知识，可以很容易地证出来，请参考群论的相关书籍。
<定理>
若P，Q是相异素数，N＝P×Q，M＝（P－1）（Q－1）。D×E≡1 MOD M。A是任意一

个正整数，B ≡A^E MOD N，C ≡B^D MOD N。
则 C≡A MOD N。
<证明>
因为D×E≡1 MOD M，所以D×E=kM＋1，其中k是整数。所以，C≡B^D≡（A^E）^D≡

A^（E×D）≡A^（kM＋1）MOD N。
1) 如果A不是P的倍数，也不是Q的倍数。
则A^（P－1）≡1 MOD P =>A^（k（P－1）（Q－1））≡1 MOD P
A^（Q－1）≡1 MOD Q =>A^（k（P－1）（Q－1））≡1 MOD Q
所以P，Q均能整除A^（k（P－1）（Q－1））－1
=>A^（k（P－1）（Q－1））≡1 MOD PQ
=>C ≡A^（k（P－1）（Q－1）＋1）≡A MOD N
2) 如果A是P的倍数，但不是Q的倍数。
则A^（Q－1）≡1 MOD Q =>A^（k（P－1）（Q－1））≡1 MOD Q
=>C ≡A^（k（P－1）（Q－1）＋1）≡A MOD Q => Q | C－A
因 P | A => C ≡A^（k（P－1）（Q－1）＋1）≡0 MOD P
=> P | C－A，所以，PQ | C－A =>C≡A MOD N
3) 如果A是Q的倍数，但不是P的倍数（同上）。
4) 如果A同时是P和Q的倍数。
则PQ | A=> C ≡A^（k（P－1）（Q－1）＋1）≡0 MOD PQ
=> PQ | C－A => C≡A MOD N


B. 素数的产生
1. 实际考虑
1) 产生一个n位的随机大数p。
2) 设p的最高位和最低位都为1，以确保达到要求的长度，和确保该大数是奇

数。
3) 检查以确保p不能被任何小素数整除，如3，5，7，11等等。有效的方法是

测试小于2000的素数。使用字轮方法更快。
4) 对某随机数a运行Rabin-Miller检测，如果P通过，则另外产生一个随机数

a,在测试。选取较小的a值，以保证速度。做5次Rabin-Miller测试，如果有一次失

败，返回第一步，从新产生p，再测试。
2. Rabin-Miller素数检测算法
对一个待测的随机大数p，计算p的有效位n，比如p＝0x5A，则n＝7。
1) 选择一个小于p的随机数a。
2) 令z＝1，i＝n。
3) x＝z，z＝d×d MOD N。
4) 如果z＝1并且x<>1并且x<>p－1，那么测试通不过，p不是素数。
5) 如果p－1的第i位为1，令z＝z×a MOD p，i＝i－1。如果i>=0，则转到第

三步。
6) 如果z<>1，则同不过检测，p不是素数，否则通过检测，p可能为素数。
这个测试算法速度快，但是有很小的概率会出错。如果用5次该算法进行检测，p是

素数的概率可达到99.9％（1－（0.25）^5）。这里给出的该算法的过程可能和有些

资料上不太一样，本人从经验上觉得这里给出的更可靠些。
3. 其它素数检测算法
除了常用的Rabin-Miller检测，还有Solovag-Strassen检测方法，Lehmann检测方法

等。这里不作详细介绍了，请参考应用密码学的相关内容。


C. 大数的运算
1. 大数的运算原理
RSA算法依赖于大数的运算，目前主流RSA算法都建立在512位到1024位的大数运算之

上，所以我们首先需要掌握大数（比如1024位）的运算原理。
大多数的编译器只能支持到32位（或64位）的整数运算，即我们在运算中所使用的

整数必须小于等于32位，即0xFFFFFFFF，这远远达不到RSA的需要，于是需要专门建

立大数运算库，来解决这一问题。
最简单的办法是将大数当作字符串进行处理，也就是将大数用10进制字符数组进行

表示，然后模拟人们手工进行“竖式计算”的过程编写其加减乘除函数。但是这样

做效率很低。当然其优点是算法符合人们的日常习惯，易于理解。
另一种思路是将大数当作一个二进制流进行处理，使用各种移位和逻辑操作来进行

加减乘除运算，但是这样做代码设计非常复杂，可读性很低，难以理解也难以调试

。
这里我们采用了一种介于两者之间的思路：将大数看作一个N进制数组，对于目前的

32位系统而言，N可以取2的32次方，即0x100000000，假如将一个1024位的大数转化

成0x10000000进制，它就变成了32位，而每一位的取值范围是0－0xFFFFFFFF。我们

正好可以用一个无符号长整数来表示这一数值。所以1024位的大数就是一个有32个

元素的unsigned long数组。而且0x100000000进制的数组排列与2进制流对于计算机

来说，实际上是一回事，但是我们完全可以针对unsigned long数组进行“竖式计算

”，而循环规模被降低到了32次之内，并且算法很容易理解。
但考虑到乘法和除法，都要进行扩展才能进行快速的计算（如果把除法变减法而不

扩展，速度将慢的无法忍受）。所以我们将N取2的16次方的，即0xFFFF。每一位用

unsigned short来表示，当进行乘除运算时，将short扩展成long，这是编译器所支

持的，所以运算起来，比较快。
2. 大数的各种运算
这些运算都是常见的，同时也是常用的。要实现RSA算法，就必须先实现大数的这些

运算。
1) 大数的比较。两个无符号或有符号的大数进行大小比较。大数和一般整数

进行比较。大于，等于，小于，返回值各异，以区别比较结果。
2) 大数的赋值。将一个大数的值，符号等，逐位赋给另一个大数。将一般整

数的值，符号等赋给一个大数。
3) 大数的加法。两个大数从低位到高位逐位相加，要考虑到进位的问题。或

大数与一般整数的相加。
4) 大数的减法。两个大数从低位到高位逐位相减，要考虑到借位的问题。或

大数与一般整数的相减。
5) 大数的乘法。两个大数的乘法，从一个大数的低位到高位，逐位与另一个

大数相乘，然后将结果低位对齐相加，要考虑进位，类似于普通的竖式乘法。或大

数与一般整数的乘法。
6) 大数的除法。两个大数的除法，从一个大数的高位到低位，逐步与另一个

大数相除，要考虑借位，类似于普通的竖式除法。或大数与一般整数的乘法。
7) 大数的取余。两个大数的取余，类似于大数的除法，只是当除到底时，返

回的是余数而已，也存在借位的问题。或大数于一般整数的取余。
8) 大数的欧氏算法。它是已知大数A、B，求满足AX≡1 MOD B的X，是最大公

约数算法的扩展，同样用辗转相除法。再递归的过程中，两个参数都要用到，到要

变化的。具体算法请参考源代码。
9) 大数的蒙氏算法。它是已知大数A、B和C，求X=A^B MOD C的值。要对指数

进行逐渐降阶，直到变成2次方，也就是转换成乘法和取余运算。降阶的方法和算法

的具体过程，请参考相关书籍和源代码。
10) 大数的最大公约数。求两个大数的最大公约数，用辗转相除法。

RSA算法的实现

A. 生成密钥函数
gChar gGenerateKey(gBigInt *n,gBigInt *e,gBigInt *d);
功能：该函数实现了产生密钥的功能。先产生两个随机的大素数p和q，然后计算n＝p×q，
随机产生（或固定）一个大数e，计算d，使得ed≡1 MOD (p-1)(q-1)。
参数：
n：两个大数的乘积，和e或d联合构成加密密钥或解密密钥。
e：一个大数，作为加密密钥。
d：一个和e互异的大数，作为解密密钥。
返回值：1－成功，0－失败。
B. 数据加密函数
char gEncRSA(unsigned char* indata, unsigned long inlen, gBigInt* e, gBigInt* n,/
        unsigned char *outdata, unsigned long* outlen, int flg);
功能：把待加密的明文数据indata，用加密密钥e和n进行分段加密。
加密后的数据放到提前开辟好的内存outdata中，其长度outlen不得小于((inlen+k-12)/(k-11))*k，
这里k为n的位数。
参数：
indata：指向明文数据的指针。
Inlen：传入数据的长度，即明文数据的长度。
e和n：两个大数，作为加密密钥。
Outdata：存放加密后密文数据的指针。
Outlen：传入outdata的长度，传出数据的长度，即密文数据的长度。
Flg：公钥加密或私钥加密的标志，g_PUBLIC－公钥，g_PRIVATE－私钥。
返回值：1－成功，0－失败。
C. 数据解密函数
char gDecRSA(unsigned char* indata, unsigned long inlen, gBigInt* d, gBigInt* n,/
        unsigned char *outdata, unsigned long* outlen, int flg);
功能：把待解密的密文数据indata，用解密密钥e和n进行分段解密。
解密后的数据放到提前开辟好的内存outdata中，其长度outlen不得小于(inlen)*k/(k-11)，
这里k为n的位数。
参数：
indata：指向密文数据的指针。
Inlen：传入数据的长度，即密文数据的长度。
d和n：两个大数，作为加密密钥。
Outdata：存放解密后明文数据的指针。
Outlen：传入outdata的长度，传出数据的长度，即明文数据的长度。
Flg：公钥加密或私钥加密的标志，g_PUBLIC－公钥，g_PRIVATE－私钥。
返回值：1－成功，0－失败。
D. 用小素数检测
gChar gIsNotPrime(gBigInt *p);
功能：检测随机大数p能否被小素数整除。
参数：
p：待检测的随机大数。
返回值：0－可能是素数，其它－为能整除p的小素数或1。
E. Rabin-Miller检测
gChar gRabinMillerTest(gBigInt *p,int tNum);
功能：用Rabin-Miller算法，对大数p进行tNum次检测，判断p是否可能为素数。
参数：
p：待检测的随机大数。
tNum：检测的次数，它可以决定该检测的可信度。
返回值：0－失败或通不过，1－成功并通过，说明p可能是素数。
F. 随机产生大素数
gChar gGeneratePrime(gBigInt *p);
功能：随机生成一个大素数。
参数：
p：随机生成的大素数。
返回值：0－失败，1－成功。
G. 顺序搜索大素数
gChar gGeneratePrimeEx(gBigInt *p, int flg);
功能：根据标志flg，递增或递减地搜索p附近的大素数。
参数：
p：搜索到的大素数。
Flg：方向标志，g_INCREASE－递增，g_DECREASE－递减。
返回值：0－失败，1－成功。
H. 字符串与大数的转换
gChar gStr2BigInt(gBigInt *gbi,char *str,int flg);
功能：根据标志flg，将表示成10进制或16进制形式的字符串str转换成大数gbi。
参数：
gbi：转换成的大数。
Str：待转换的字符串。
Flg：标志，OCT_FLG－8进制，DEC_FLG－10进制，HEX_FLG－16进制。
返回值：0－失败，1－成功。
gChar gBigInt2Str(char *str,gBigInt *gbi,int flg);
功能：根据标志flg，将大数gbi转换成10进制或16进制形式表示的字符串str。
参数：
Str：转换成的字符串。
gbi：待转换的大数。
Flg：标志，OCT_FLG－8进制，DEC_FLG－10进制，HEX_FLG－16进制。
返回值：0－失败，1－成功。
I. 一般整数的乘除
gUInt gIntMul(gUInt *a, gUInt b);
功能：一般整数的乘法。在函数内部，对整数a和b进行扩展，然后与b相乘，把结果

放在a中返回，进位作为返回值。
参数：
a：传入一个整数，传出相乘的结果。
b：传入一个整数，作乘数。
返回值：两个整数相乘的进位。失败时，返回0。
gUInt gIntDiv(gUInt *ah, gUInt *al, gUInt b);
功能：一般整数的除法。在函数内部，将整数ah和al合并成一个扩展了的整数，
然后除以整数b，把结果的高位放到ah中，低位放到al中返回，余数作为返回值。
参数：
ah：传入一个整数，作为被除数的高位；传出相除后的结果的高位。
al：传入一个整数，作为被除数的低位；传出相除后的结果的低位。
b：传入一个整数，作除数。
返回值：相除的余数。失败时，返回0。
J. 大数的比较
gChar gCmp(gBigInt *a, gBigInt *b);
功能：对大数a和大数b进行比较，不考虑它们的符号，只比较数值。
参数：
a：传入一个大数。
b：传入另一个大数。
返回值：-1－大数a小于大数b，0－大数a等于大数b，+1－大数a大于大数b。
gChar gCmpEx(gBigInt *a, gBigInt *b);
功能：对大数a和大数b进行比较，考虑它们的符号，整数大于负数。
参数：
a：传入一个大数。
b：传入另一个大数。
返回值：-1－大数a小于大数b，0－大数a等于大数b，+1－大数a大于大数b。
gChar gCmpInt(gBigInt *a, gUInt b);
功能：大数与整数的比较。对大数a和一般整数b进行比较，不考虑它们的符号。
参数：
a：传入一个大数。
b：传入一个一般整数。
返回值：-1－大数a小于b，0－大数a等于b，+1－大数a大于b。
K. 大数的赋值
gChar gMov(gBigInt *a, gBigInt *b);
功能：将大数b的值，赋给大数a，考虑它们的符号。
参数：
a：传出赋值的结果。
b：传入另一个大数。
返回值：0－失败，１－成功。
gChar gMovInt(gBigInt *a, gUInt b);
功能：将整数b的值，赋给大数a，考虑它们的符号。
参数：
a：传出赋值的结果。
b：传入一个一般整数。
返回值：0－失败，１－成功。
L. 大数的加法
gChar gAdd(gBigInt *a, gBigInt *b);//加
功能：对大数a和大数b进行加法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相加的结果。
b：传入另一个大数。
返回值：0－失败，１－成功。
gChar gAddInt(gBigInt *a, gSInt b);
功能：对大数a和整数b进行加法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相加的结果。
b：传入一个一般整数。
返回值：0－失败，１－成功。
M. 大数的减法
gChar gSub(gBigInt *a, gBigInt *b);//减
功能：对大数a和大数b进行减法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相减的结果。
b：传入另一个大数。
返回值：0－失败，１－成功。
gChar gSubInt(gBigInt *a, gSInt b);
功能：对大数a和整数b进行减法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相减的结果。
b：传入一个一般整数。
返回值：0－失败，１－成功。
N. 大数的乘法
gChar gMul(gBigInt *a, gBigInt *b);//乘
功能：对大数a和大数b进行乘法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相乘的结果。
b：传入另一个大数。
返回值：0－失败，１－成功。
gChar gMulInt(gBigInt *a, gSInt b);
功能：对大数a和整数b进行乘法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相乘的结果。
b：传入一个一般整数。
返回值：0－失败，１－成功。
O. 大数的除法
gChar gDiv(gBigInt *a, gBigInt *b);//除
功能：对大数a和大数b进行除法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相除的结果。
b：传入另一个大数。
返回值：0－失败，１－成功。
gChar gDivInt(gBigInt *a, gSInt b);
功能：对大数a和整数b进行除法运算，考虑它们的符号。
参数：
a：传入一个大数，传出相除的结果。
b：传入一个一般整数。
返回值：0－失败，１－成功。
P. 大数的取余
gChar gMod(gBigInt *a, gBigInt *b);//取余
功能：对大数a和大数b进行取余运算，考虑它们的符号。
参数：
a：传入一个大数，传出取余的结果。
b：传入另一个大数。
返回值：0－失败，１－成功。
gChar gModInt(gBigInt *a, gSInt b);
功能：对大数a和整数b进行取余运算，考虑它们的符号。
参数：
a：传入一个大数，传出取余的结果。
b：传入一个一般整数。
返回值：0－失败，１－成功。
Q. 欧几里德算法
gChar gEuc(gBigInt *a,gBigInt *b);//欧几里德算法
功能：对大数a和大数b执行欧几里德算法，结果放到a中。即用辗转相除的思想，
求满足aX≡1 MOD b的X的值，也即求满足aX－bY＝1的较小的大数X和Y的值。
其中用到了递归的思想，运算后，大数a和大数b，都改变了。
参数：
a：传入一个大数，传出满足aX－bY＝1的X的值。
b：传入另一个大数，传出满足aX－bY＝1的Y的值
返回值：0－失败，１－成功。
R. 蒙格马利算法
gChar gMon(gBigInt *a, gBigInt *b, gBigInt *c);//蒙格马利算法
功能：对大数a，大数b和大数c执行蒙格马利算法，结果放到a中。
即用指数降阶的思想，求算式a^b MOD c的值。至于降阶的原理，请参考相关书籍。
参数：
a：传入一个大数，传出运算的结果。
b：传入另一个大数。
c：传入第三个大数。
返回值：0－失败，１－成功。
S. 最大公约数
gChar gGcd(gBigInt *a, gBigInt *b);//最大公约数
功能：用辗转相除法，求大数a和大数b的最大公约数。
参数：
a：传入一个大数，传出a和b的最大公约数
b：传入另一个大数。
返回值：0－失败，１－成功。
卡列宁的微笑
是时候了
一个丑丑的模运算算法实现
#最近我发现了一个问题，就是关于浮点数求模的问题。该问题描述如下：给定两个浮点数#####（double）a,b,求a modb.浮点数的模运算是仿照整数的模运算得来的。它的结果的取值#范围应该是在0到b之间，可以无限接近b。它是基于形如如下的应用：f(x)是一个周期为T的##函数，T是一个浮点数，给定一个数x求f(x).如果x是一个很大的数，直接求f(x)的值是比较##麻烦的，往往得不到正确的结果，如果有了浮点数求模运算，问题就可以得到很好解决。现##在的关键就是如何解决浮点数的求模问题。


最好想的做法帖子下面已经有人以C++实现鸟，我用伪代码表示如下：
while a is not in (b(n-1),bn]
n++;
道理如此的简单，两行就够说明了，但是，在大数模运算中，这个算法的效率是极其低下的，基本上要做a/b次乘法
改进的办法我查到的虽然有几种，但是思路上都是以空间换时间，我自己用python给出了一个丑陋的实现，用python是因为正在学，而且python表达起来比较清晰一点，下面就把源码bia出来，我的基本思路是：
    通过一定方式构造一个b的整倍数数组，这个数组可以以极快的速度逼近a，然后让a和这个数组进行比较，递归到数组起始元素的时候（也就是b），就可以得到模运算的结果
另外，有朋友说TAOPC中有关于这个东东的内容，烦请哪位朋友赐教，我没有那本书。
 1 # antmod.py
 2 def antmod (a,b,ANTLIMIT):
 3     ant=0.0
 4     if a<b:
 5         ant=a          
 6         a=b
 7         b=ant              #确保a总是比b大
 8     li=[]                  
 9     temp=1                 
10     for i in range(1,ANTLIMIT):     #构造比较数组
11         temp=temp*b
12         li=li+[temp]     
13     else:
14         k=0
15     while k!=ANTLIMIT-2:            #当a还比b大的时候进行循环
16         if a>li[-1-k]:         
17             a-=li[-1-k]      #a和当前数组元素进行比较，大则减之，若a比当前元素小，
18         else:                    #将当前元素前移一位
19             k+=1
20     else:
21         return a                  #返回结果     
22
23 print antmod(10000,3,5)                   #小测试一下
24
25  


       
       
    一上面的数据为例，若使用群组讨论中提供的算法，要做10000/3次乘法运算，而使用上面的代码，需要5次乘法和几十次减法运算，效率上的区别是巨大的。
   不过，上面的代码只是一个大致的实现，结构上是非常不合理的，比方生产数组的接口应当提供给用户，这个貌似用C++的模板机制比较好一点，python也未必不能实现。这里只是稍做说明。