大整数乘法的Karatsuba算法实现

编程知识分享 11-10 218

前言：

此刻兄弟们对“高精度乘法程序设计”都比较注意，咱们都想要了解一些“高精度乘法程序设计”的相关文章。那么小编同时在网络上汇集了一些对于“高精度乘法程序设计””的相关资讯，希望兄弟们能喜欢，各位老铁们一起来学习一下吧！

两个整数相乘，使用基本算法，时间复杂度为O(n^2) ，这对于日趋庞大的数据来说是很慢的，目前比较常见的一种大整数的快速算法是 Karatsuba算法，当然他不是最快的，但是比基本算法要好的多，时间复杂度为O(n^1.59)，在密码运算中相差是很大的。

现在考虑分治算法。取m = (n+1)/2，把x写成10^m*a+b的形式，y写成10^m*c+d的形式，则a, b, c, d都是m位整数(如果不足m位，前面可以补0)。

递归方程为T(n) = 4T(n/2) + O(n)，其中系数4为四次乘法ac, bd, bc, ad，附加代价n为最后一个return语句的两次高精度加法。方程的解为T(n) = O(n^2)，和暴力乘法没有区别。

Anatolii Karatsuba在1962年提出一个改进方法（并由Knuth改进）：用ac和bd计算bc + ad，即：

bc + ad = ac + bd - (a - b) * (c - d)

这样一来，只需要进行三次递归乘法，即递归方程变为了T(n) = 3T(n/2)+O(n)，解为T(n) = O(nlog3) = O(n^1.585)，比暴力乘法快。

计算整数乘法的最快算法是基于FFT的，它的时间复杂度为O(n log n)。下面是Karatsuba's multiplication algorithm 大整数的快速乘法实现，使用递归。

理论上应该计算速度很快，但是在测试中，内存消耗巨大，速度也较慢，还希望大家帮帮忙，优化一下。和FFT_MULT相比，相差近100倍！

* 时间：2013年10月11日16:16:19

* 功能：Karatsuba's multiplication algorithm 大整数的快速乘法实现

* 输入：两个大整数f&g，位数不大于2048

* 输出：f*g

# include "miracl.h" //调用大整数库miracl.h

# include <time.h>

# include "math.h"

# define TIMES 1

big mult_k(big a, big b, long len); //快速乘法

long length(big a); //计算数字的长度

long max(long a, long b); //求较大值

long max_len(big a, big b); //求长度n，n满足是2的次幂

/*主函数*/

int main(void) {

int i=0, j=0, k=0;

big f, g, res;

long len = 0;

FILE *fp;

clock_t tBegin1, tEnd1;

clock_t tBegin2, tEnd2;

clock_t tBegin3, tEnd3;

miracl *mip = mirsys(4096, 16); //最大4096位,输入输出使用16进制

//初始化

f = mirvar(0);

g = mirvar(0);

res = mirvar(0);

fp = fopen("input.txt", "r+");

mip->IOBASE = 16;

cinnum(f, fp);

cinnum(g, fp);

fclose(fp);

printf("f:\n");

cotnum(f, stdout);

printf("g:\n");

cotnum(g, stdout);

//TIMES次普通大整数乘法

tBegin1 = clock();

for (i=0; i<TIMES; i++)

multiply(f, g, res);

cotnum(res, stdout);

tEnd1 = clock();

//TIMES次miracl FFT大整数乘法

tBegin2 = clock();

for (j=0; j<TIMES; j++)

fft_mult(f, g, res);

printf("FFT_MULT:\n");

cotnum(res, stdout);

tEnd2 = clock();

//TIMES次自写大整数乘法

tBegin3 = clock();

len = max_len(f, g);

printf("length: %ld\n", len);

for (j=0; j<TIMES; j++)

copy(mult_k(f,g, len), res);

tEnd3 = clock();

//输出

printf("XDC_MULT:\n");

cotnum(res, stdout);

//释放内存

mirkill(f);

mirkill(g);

mirkill(res);

mirexit();

//printf("\n\n进行%d次%ld比特的普通大整数乘法运算所消耗的时间为:%ld ms\n\n", TIMES, len, tEnd1 - tBegin1);

printf("\n\n进行%d次%ld比特的miracl FFT大整数乘法运算所消耗的时间为:%ld ms\n\n", TIMES, len, tEnd2 - tBegin2);

printf("\n\n进行%d次%ld比特的自写大整数乘法运算所消耗的时间为:%ld ms\n\n", TIMES, len, tEnd3 - tBegin3);

return (0);

}

//计算最大长度

long max_len(big a, big b) {

long len = 0, n = 1; //保存数字的长度,二进制

len = max(length(a), length(b)); //计算数字长度，取较大值

//len变为2的方幂

while ((pow(2, n) - len) < 0) {

n++;

}

len = (long)pow(2, n);

return (len);

}

/****************************************K氏乘法****************************************************/

big mult_k(big a, big b, long len) {

big res, tmp_0; //保存结果

big a1, b1, a0, b0, power_2; //保存分解的因子

big a0_b0, a1_b1;

long m = 0; //保存数字的长度

//初始化

res = mirvar(0);

a1 = mirvar(0);

b1 = mirvar(0);

a0 = mirvar(0);

b0 = mirvar(0);

power_2 = mirvar(0);

tmp_0 = mirvar(0);

a0_b0 = mirvar(0);

a1_b1 = mirvar(0);

if (len == 1) {

multiply(a, b, res);

//释放内存，内存泄露

mirkill(b1);

mirkill(b0);

mirkill(a1);

mirkill(a0);

mirkill(tmp_0);

mirkill(power_2);

mirkill(a0_b0);

mirkill(a1_b1);

return (res);

}

else

{

m = len / 2;

//优化,减少中间变量可以减少内存消耗

/*分解 a & b*/

sftbit(a, (-1)*m, a1); //移位，相当于除2的m次方

sftbit(a1, m, power_2);

negify(power_2, power_2);

add(a, power_2, a0);

sftbit(b, (-1)*m, b1); //移位，相当于除2的m次方

sftbit(b1, m, power_2);

negify(power_2, power_2);

add(b, power_2, b0);

copy(mult_k(a1, b1, m), a1_b1);

copy(mult_k(a0, b0, m), a0_b0);

add(a0, a1, a0); //计算a0+a1

add(b0, b1, b0); //计算b0+b1

copy(mult_k(a0, b0, m), a0);

//计算返回值

sftbit(a1_b1, len, tmp_0); //移位，相当于乘2的len次方

add(tmp_0, a0_b0, tmp_0);

//求负值

negify(a0_b0, a0_b0);

negify(a1_b1, a1_b1); //取负值

add(a0_b0, a1_b1, a0_b0); //负值相加

add(a0_b0, a0, a0_b0);

sftbit(a0_b0, m, a0_b0);

add(a0_b0, tmp_0, res);

}

//释放内存

mirkill(b1);

mirkill(b0);

mirkill(a1);

mirkill(a0);

mirkill(tmp_0);

mirkill(power_2);

mirkill(a0_b0);

mirkill(a1_b1);

return (res);

}

/***********************计算整数的长度**************************************************/

long length(big a) {

long i = 0;

big tmp;

tmp = mirvar(0);

expb2(i, tmp);

while (compare(tmp, a) == -1) {

i++;

expb2(i, tmp);

}

mirkill(tmp);

return (i);

}

//max

long max(long a, long b){

return (a >= b ? a : b);

}

测试环境：

---------------------------------------------------------------------------------

操作系统：windows7， x86

硬件：2G内存，主频2.67GHz

编译平台：VC6.0企业版

---------------------------------------------------------------------------------

在VC6.0中输出测试结果为：

---------------------------------------------------------------------------------

F05085869EF4BA2514D08635E180E138DCD2AAAF1B04C69A4C3A9B612A6FAF9784393B5B49026FEA

2F0E244D84506A7A1D44B8745CE4B9B0C83668FD83BADEFC2A6EEC3D80BA5A3CEB1CB538C25199B0

5E3E3535F3276020F53C8E9D2B518465BD2F6322C1751A00C6EF5186614D9EC955841B2CCFD59882

853E4131233BC2E3D98E5FC36267464CE6947FEEE0EC8BC7AA611AD15D68F234BAC62C18C9DEF38B

A135550D54EBCD179EA40F377A01066B13E61FF8C9639B2D3A19EC7B8CC58877F7266FDDDC776C56

3D277DB0204C9CE7213D87E76750478531E3B09685629B1B9FEB06E118A5F3E978F8AED1D0C202A5

728021831A5012D43DE53C9CAFFF4E1D

D98E5FC36267464CE6947FEEE0EC8BC7AA611AD15D68F234BAC62C18C9DEF38BA135550D54EBCD17

9EA40F377A01066B13E61FF8C9639B2D3A19EC7B8CC58877F7266FDDDC776C563D277DB0204C9CE7

213D87E76750478531E3B09685629B1B9FEB06E118A5F3E978F8AED1D0C202A5728021831A5012D4

3DE53C9CAFFF4E1DD98E5FC36267464CE6947FEEE0EC8BC7AA611AD15D68F234BAC62C18C9DEF38B

A135550D54EBCD179EA40F377A01066B13E61FF8C9639B2D3A19EC7B8CC58877F7266FDDDC776C56

3D277DB0204C9CE7213D87E76750478531E3B09685629B1B9FEB06E118A5F3E978F8AED1D0C202A5

728021831A5012D43DE53C9CAFFF4E1D

FFT_MULT:

CC39E7BE78AC0D920B95B029E2005B1DB44E62400D8C455AE03E02FAD84143091F245DDBAD74DCDF

32B5C19991E06B04E8A06F716943966FB2C53F62129DADCF841A24094C453D407FB65C6C7B2140AC

494BCE1DE2C574274D6BE136D8D3E0B6AA7F0E625B50C0663EA1EB84D68D16C9F7C695AE1BF5C545

52B4D8446E0D212DAFB66B4EE69F5BE5E1208D2B655D9CB54F68B239A88C9EBBB51C292C2D967BEB

936FB96D29FFC64CE50A0BEDF66020C2D3B6982D1DED645C6603EA7C8E2477C8CE76B1E3B8669D54

C9A29B2F50A98D598F0CF44166E0C538817B37177FE46171AF5094D4424DE3EEE0E99FC3BE237320

91623AA160D5BB56F7B641549845911438F2963CBDD7FAFC854E3DC24F88E4EC04D93E59150476F2

0CA4B384A4E2042A1FA261A4EF3C7D10A51976C090B3624259A9F42DBCCC1975C9278C66FAC6857E

01BA4D7FA09F91FB795ECD0D2EF2FBA3B9B086A51F490FE4282898F426CBF4FF11C86F84FD5F4A8D

E2514778C625189500E0647B8B61C67BDBA73FD085D41F30557612AC4FE4ACA8AFC360C0CC2BA354

69BEEE5F7A041D9137C68D534F8CCB47AB57061372B193A2F2C52C6C2C33AC846E93CB7208224B89

15E8E14C7F3FBB84B75DBFA5347E31E72F728E4A596AAEF673EF60819B2DBED2F41943137FBB7444

0CF6E9131662270540099339DE8EBC3E6744BF884681D06A36A5D6C05B9BAF49

length: 2048

xdc_mult: