大整数乘法的Karatsuba算法实现

最新推荐文章于 2024-03-28 10:38:45 发布

原创最新推荐文章于 2024-03-28 10:38:45 发布 · 6.7k 阅读

5 ·

CC 4.0 BY-SA版权

文章标签：

#大整数乘法 #快速乘法 #Karatsuba #算法

信息安全专栏收录该内容

2 篇文章

订阅专栏

本文介绍了Karatsuba算法作为大整数乘法的快速算法，其时间复杂度为O(n^1.59)，优于基本算法的O(n^2)。通过将大整数分解并进行三次递归乘法，实现效率提升。然而，与基于FFT的更快算法相比，Karatsuba算法在实际应用中可能存在内存消耗和速度的问题。

两个整数相乘，使用基本算法，时间复杂度为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)。【参考：http://blog.youkuaiyun.com/jiyanfeng1/article/details/8543846】

下面是Karatsuba's multiplication algorithm 大整数的快速乘法实现，使用递归。

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

/*
* 时间：2013年10月11日16:16:19
* 作者：xdc
* 功能：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中输出测试结果为：
---------------------------------------------------------------------------------
f:
F05085869EF4BA2514D08635E180E138DCD2AAAF1B04C69A4C3A9B612A6FAF9784393B5B49026FEA
2F0E244D84506A7A1D44B8745CE4B9B0C83668FD83BADEFC2A6EEC3D80BA5A3CEB1CB538C25199B0
5E3E3535F3276020F53C8E9D2B518465BD2F6322C1751A00C6EF5186614D9EC955841B2CCFD59882
853E4131233BC2E3D98E5FC36267464CE6947FEEE0EC8BC7AA611AD15D68F234BAC62C18C9DEF38B
A135550D54EBCD179EA40F377A01066B13E61FF8C9639B2D3A19EC7B8CC58877F7266FDDDC776C56
3D277DB0204C9CE7213D87E76750478531E3B09685629B1B9FEB06E118A5F3E978F8AED1D0C202A5
728021831A5012D43DE53C9CAFFF4E1D
g:
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:
CC39E7BE78AC0D920B95B029E2005B1DB44E62400D8C455AE03E02FAD84143091F245DDBAD74DCDF
32B5C19991E06B04E8A06F716943966FB2C53F62129DADCF841A24094C453D407FB65C6C7B2140AC
494BCE1DE2C574274D6BE136D8D3E0B6AA7F0E625B50C0663EA1EB84D68D16C9F7C695AE1BF5C545
52B4D8446E0D212DAFB66B4EE69F5BE5E1208D2B655D9CB54F68B239A88C9EBBB51C292C2D967BEB
936FB96D29FFC64CE50A0BEDF66020C2D3B6982D1DED645C6603EA7C8E2477C8CE76B1E3B8669D54
C9A29B2F50A98D598F0CF44166E0C538817B37177FE46171AF5094D4424DE3EEE0E99FC3BE237320
91623AA160D5BB56F7B641549845911438F2963CBDD7FAFC854E3DC24F88E4EC04D93E59150476F2
0CA4B384A4E2042A1FA261A4EF3C7D10A51976C090B3624259A9F42DBCCC1975C9278C66FAC6857E
01BA4D7FA09F91FB795ECD0D2EF2FBA3B9B086A51F490FE4282898F426CBF4FF11C86F84FD5F4A8D
E2514778C625189500E0647B8B61C67BDBA73FD085D41F30557612AC4FE4ACA8AFC360C0CC2BA354
69BEEE5F7A041D9137C68D534F8CCB47AB57061372B193A2F2C52C6C2C33AC846E93CB7208224B89
15E8E14C7F3FBB84B75DBFA5347E31E72F728E4A596AAEF673EF60819B2DBED2F41943137FBB7444
0CF6E9131662270540099339DE8EBC3E6744BF884681D06A36A5D6C05B9BAF49


进行1次2048比特的miracl FFT大整数乘法运算所消耗的时间为:177 ms



进行1次2048比特的自写大整数乘法运算所消耗的时间为:16001 ms

Press any key to continue
---------------------------------------------------------------------------------
*/

这是自己第一次使用大数库写代码，请大家多指教。