[转]大整数算法[11] Karatsuba乘法

★ 引子

        前面两篇介绍了 Comba 乘法,最后提到当输入的规模很大时,所需的计算时间会急剧增长,因为 Comba 乘法的时间复杂度仍然是 O(n^2)。想要打破乘法中 O(n^2) 的限制,需要从一个完全不同的角度来看待乘法。在下面的乘法算法中,需要使用 x 和 y 这两个大整数的多项式基表达式 f(x) 和 g(x) 来表示。

令 f(x) = a * x + b,g(x) = c * x + d,h(x) = f(x) * g(x)。这里的 x 相当于一个基,比如十进制下,123456 可以表示成 123 * 100 + 456,这时 x 就是 100 了。

★ Karatsuba 乘法原理

         既然当输入规模 n 增加时计算量会以 n^2 增加,那么考虑使用分治的方式把一个规模较大的问题分解成若干个规模较小的问题,这样解决这几个规模较小的问题就会比较容易了。

对于 x 和 y 的乘积 z = x * y,可以把 x 和 y 都拆成两段:x 的左半段和右半段分别是 a 和 b,y 的左半段和右半段分别是 c 和 d。

现在考虑乘积 h(x) = f(x) * g(x):

h(x) = (a * x + b) * (c * x + d)

= a * c * (x^2) + (a * d + b * c) * x + b * d

可以看到,计算 h(x),需要计算 4 个一半大小的乘法,3 次加法,乘以 x 或者乘以 x^2 可以通过移位实现,所有的加法和移位共用 O(n) 次计算。

设 T(n) 是两个 n 位整数相乘所需的运算总数,则递归方程有:

T(n) =  4 * T(n / 2) + O(n)          当 n > 1

T(n) = O(1)                                       当 n = 1

解这个递归方程,得到 T(n) = O(n^2),即时间复杂度和 Comba 乘法是一样的,没有什么改进。要想降低计算的复杂度,必须减少乘法的计算次数。

注意到 a * d + b * c 可以用 a * c 和 b * d 表示: (a + b) * (c + d) - a * c - b * d。

原式 = (a + b) * (c + d) - a * c - b * d

= a * c + a * d + b * c + b * d - a *c - b *d

= a * d + b * c

上边的换算,说明计算 a * d + b * c 可以通过两次加法一次乘法搞定,这样总的乘法次数就减少到 3 次,列出新的递归方程有:

T(n) = 3 * T(n / 2) + O(n)         当 n > 1

T(n) = O(1)                                     当 n = 1

解递归方程得到 T(n) = O(n^log3),注意这里的 log 是以 2 为底,近似计算,时间复杂度为:T(n) = O(n^1.585),比 Comba 乘法的 O(n^2) 要小。

以上就是使用分治的方式计算乘法的原理。上面这个算法,由 Anatolii Alexeevitch Karatsuba 于1960年提出并于1962年发表,所以也被称为 Karatsuba 乘法。

★ 实现思路

原理弄明白了,现在整理一下思路:

1. 拆分输入:

计算分割的基:B = MIN(x->used, y->used) / 2

x0,y0 分别存储 x 和 y 的低半部分,x1,y1 分别存储 x 和 y 的高半部分。

2. 计算三个乘积:

x0y0 = x0 * y0           //递归调用乘法 bn_mul_bn

x1y1 = x1 * y1

t1 = x0 + x1

x0 = y0 + y1

t1 = t1 * x0

3. 计算中间项:

x0 = x0y0 + x1y1

t1 = t1 - x0

4. 计算最终乘积:

t1 = t1 * (2^(n * B))                    //左移 B个数位

x1y1 = x1y1 * (2^(2 * B))         //左移 2 * B 个数位

t1 = x0y0 + t1

z = t1 + x1y1                        //最终结果

★ 实现

         根据上面的思路,Karatsuba乘法的实现代码如下:

static int bn_mul_karatsuba(bignum *z, const bignum *x, const bignum *y)
{
    int ret;
    size_t i, B;
    register bn_digit *pa, *pb, *px, *py;
    bignum x0[1], x1[1], y0[1], y1[1], t1[1], x0y0[1], x1y1[1];

    B = BN_MIN(x->used, y->used);
    B >>= 1;

    BN_CHECK(bn_init_size(x0, B));
    BN_CHECK(bn_init_size(x1, x->used - B));
    BN_CHECK(bn_init_size(y0, B));
    BN_CHECK(bn_init_size(y1, y->used - B));
    BN_CHECK(bn_init_size(t1, B << 1));
    BN_CHECK(bn_init_size(x0y0, B << 1));
    BN_CHECK(bn_init_size(x1y1, B << 1));

    x0->used = y0->used = B;
    x1->used = x->used - B;
    y1->used = y->used - B;

    px = x->dp;
    py = y->dp;
    pa = x0->dp;
    pb = y0->dp;

    for(i = 0; i < B; i++)
    {
        *pa++ = *px++;
        *pb++ = *py++;
    }

    pa = x1->dp;
    pb = y1->dp;

    for(i = B; i < x->used; i++)
        *pa++ = *px++;

    for(i = B; i < y->used; i++)
        *pb++ = *py++;

    bn_clamp(x0);
    bn_clamp(y0);

    BN_CHECK(bn_mul_bn(x0y0, x0, y0));
    BN_CHECK(bn_mul_bn(x1y1, x1, y1));

    BN_CHECK(bn_add_abs(t1, x0, x1));
    BN_CHECK(bn_add_abs(x0, y0, y1));
    BN_CHECK(bn_mul_bn(t1, x0, t1));

    BN_CHECK(bn_add_abs(x0, x0y0, x1y1));
    BN_CHECK(bn_sub_abs(t1, t1, x0));

    BN_CHECK(bn_lshd(t1, B));
    BN_CHECK(bn_lshd(x1y1, B << 1));

    BN_CHECK(bn_add_abs(t1, x0y0, t1));
    BN_CHECK(bn_add_abs(z, t1, x1y1));

clean:
    bn_free(x0);
    bn_free(x1);
    bn_free(y0);
    bn_free(y1);
    bn_free(t1);
    bn_free(x0y0);
    bn_free(x1y1);

    return ret;
}

  

上面的代码,需要很多临时的 bignum 变量,但由于一开始就知道各个 bignum 的大小,所以使用 bn_init_size 函数初始化并且分配指定的数位,避免后面再进行内存的重新分配了,节约了时间。

算法一开始将输入的 x 和 y 拆分成两半,使用三个循环搞定。为了提高效率,这里在指针的前边加上了 register 关键字来暗示在执行过程中尽量把这几个指针变量放到 CPU 的寄存器中,以此来加快变量的访问速度。

乘法的操作是递归调用 bn_mul_bn 函数,这个是有符号数乘法的计算函数,后面会讲,当递归调用到某一个临界点后,乘法的计算会直接调用 Comba 方法进行计算,而不是一直使用 Karatsuba 递归下去。 bn_mul_bn 的函数原型是:int bn_mul_bn(bignum *z, const bignum *x, const bignum *y);

需要注意的是,每个计算操作(加法,减法,乘法和移位)在执行过程中都有可能出错,所以必须加上 BN_CHECK 宏进行错误检查,一旦函数调用出错,调到 clean 后面执行内存清理操作。

★ 分割点

虽然 Karatsuba 乘法执行时所需的单精度乘法比 Comba 方法少,但是也多了一项 O(n) 级别的开销来解一个方程组,用于计算中间项以及合并最后的结果,这就使得 Karatsuba 乘法在应对输入比较小的数字时所需的计算时间会更多。因此在实际操作中,递归计算到一定大小后,就应该改用 Comba 方法计算了。在 bn_mul_bn 函数中,分割点大小是 80(bn_digit 字长为 32 bit 时)或 64(bn_digit 字长为 64 bit 时),当输入两个数的规模有一个小于分割点时,就应该改用 Comba 方法计算乘法,只有当两个数的规模大于或等于分割点才使用 Karatsuba 方法递归计算。

★ 总结

Karatsuba 算法是比较简单的递归乘法,把输入拆分成 2 部分,不过对于更大的数,可以把输入拆分成 3 部分甚至 4 部分。拆分为 3 部分时,可以使用 Toom-Cook 3-way 乘法,复杂度降低到 O(n^1.465)。拆分为 4 部分时,使用 Toom-Cook 4-way 乘法,复杂度进一步下降到 O(n^1.404)。对于更大的数字,可以拆成 100 段,使用快速傅里叶变换,复杂度接近线性,大约是 O(n^1.149)。可以看出,分割越大,时间复杂度就越低,但是所要计算的中间项以及合并最终结果的过程就会越复杂,开销会增加,因此分割点上升,对于公钥加密,暂时用不到太大的整数,所以使用 Karatsuba 就合适了,不用再去弄更复杂的递归乘法。

下一篇将使用 Comba 方法和 Karatsuba 方法构建有符号数的乘法。

回到本系列目录

版权声明
原创博文,转载必须包含本声明,保持本文完整,并以超链接形式注明作者Starrybird和本文原始地址:http://www.cnblogs.com/starrybird/p/4445566.html

原文地址:https://www.cnblogs.com/wlzy/p/8697994.html

时间: 2024-10-10 17:23:47

[转]大整数算法[11] Karatsuba乘法的相关文章

大整数算法[11] Karatsuba乘法

★ 引子         前面两篇介绍了 Comba 乘法,最后提到当输入的规模很大时,所需的计算时间会急剧增长,因为 Comba 乘法的时间复杂度仍然是 O(n^2).想要打破乘法中 O(n^2) 的限制,需要从一个完全不同的角度来看待乘法.在下面的乘法算法中,需要使用 x 和 y 这两个大整数的多项式基表达式 f(x) 和 g(x) 来表示. 令 f(x) = a * x + b,g(x) = c * x + d,h(x) = f(x) * g(x).这里的 x 相当于一个基,比如十进制下,

大整数算法[08] Comba乘法(原理)

★ 引子          原本打算一篇文章讲完,后来发现篇幅会很大,所以拆成两部分,先讲原理,再讲实现.实现的话相对复杂,要用到内联汇编,要考虑不同平台等等. 在大整数计算中,乘法是非常重要的,因为在公钥密码学中模幂运算要频繁使用乘法,所以乘法的性能会直接影响到模幂运算的效率.下面将会介绍两种乘法:基线乘法和 Comba 乘法,尽管他们的原理和计算看起来十分类似,而且算法的时间复杂度都是 O(n^2),但是他们的效率差别是很大的. ★ 基线乘法 (Baseline Multiplication

大整数算法[10] Comba乘法(实现)

★ 引子 上一篇文章讲了 Comba 乘法的原理,这次来讲讲如何实现.为了方便移植和充分发挥不同平台下的性能,暂时用了三种不同的实现方式: 1.单双精度变量都有的情况. 2.只有单精度变量的情况. 3.可以使用内联汇编的情况. 前面已经介绍了 Comba 乘法的原理和实现思路,为了方便,再把它贴到这里: 计算 c = a * b,c0,c1,c2 为单精度变量. 1. 增加 c 到所需要的精度,并且令 c = 0,c->used = a->used + b->used. 2. 令 c2

大整数算法[00] 概述

★ 为啥要做这个 早在大一的时候,我便对密码学产生兴趣.那时在计算机导论后面看到RSA加密的计算原理,觉得十分有趣,于是就很想自己实现一个RSA加密,不过我很快就放弃了,因为实在搞不定那超长的整数计算.C里面最长的整数类型也就64位,对于动辄就1024位的RSA密钥,这连个零头都没有.为了完成这个目标,我便开始琢磨着弄一个用来计算大整数的库.原本我也打算使用别人已经写好的大数库,不过最终还是决定自己搞一个,因为凡是效率高速度快的大数库 (OpenSSL, Crypto++, GMP),要么使用的

大整数算法[01] 大整数的表示和相关定义

★ 相关的数据类型定义 在干正事之前,先定义好各种数据类型还是很有必要的,避免在以后的编码中引起混乱. uintX   X位无符号整形,如uint32表示32位无符号整形 intX    X位有符号整形,如int32表示32位有符号整形 基本数据类型定义: #ifdef _MSC_VER            typedef __int8              int8;            typedef __int16             int16;            typ

大整数算法[12] 有符号乘法

★ 引子 前面三篇文章讲了 Comba 乘法和 Karatsuba 乘法,有了这两个算法,就可以很轻松的构造有符号数乘法. 顺便提一下:讲Comba 乘法的实现的时候,给出了 x86 环境下的内联汇编实现,最近添加了 GCC x64 环境的内联汇编,已经补充到帖子当中. ★ 实现         有符号数的乘法,基本实现是这样:大的整数用 Karatsuba 乘法搞定,小的整数用 Comba 乘法搞定.对于大的整数,Karatsuba 乘法会不断递归计算,直到输入的整数小到一定规模,就改用 Co

大整数算法[13] 单数位乘法

★ 引子 最近在折腾 wxWidgets,同时拖延症又犯了,所以中断了好久.这次来讲讲单数位乘法,前面讲到 Comba 和 Karatsuba 乘法,这两个算法适合用来处理比较大的整数,但是对于一个大整数和一个单精度数相乘,其效果反而会不好,因为计算量过多.实际上单数位乘法只是基线乘法的一个特例,不存在嵌套循环进位,因此可以通过优化减少计算量.另外与完整的乘法不同的是,单数位乘法不需要什么临时变量存储和内存分配(目标精度增加除外). ★ 算法思路         单数位乘法类似于计算 12345

大整数算法

本文主要整理了几个常用的大整数的算法:大整数加法大整数乘法大整数阶乘大整数幂其实大体的思路都差不多,都是用数组来存储大整数.以下的代码仅仅实现功能,并没有充分详细的参数判断,在实际运用中,肯定是需要考虑的. 大整数相加 1 #include <stdio.h> 2 #include <string.h> 3 #define N 1000 4 void get_num(int *array) 5 { 6 char str[N] = {'\0'}; 7 int loop = 0; 8

大整数算法[08] 有符号加法和减法

★ 引子 前面几篇文章介绍了比较操作,绝对值加法和绝对值减法,现在就可以利用这几个算法构建有符号数的加减算法. ★ 有符号数加法            有符号数的加法分成两种情况:同号和异号. 1.  如果两个数同号,则执行绝对值加法,如果两个数为非负数,则结果为非负数:如果两个数都是负数,则结果也为负数. 2.  如果两个数异号,则要执行绝对值减法,用绝对值较大的数去减绝对值较小的数.最终结果 z 的符号由 x 和 y 的绝对值大小决定:如果 x 的绝对值大于或等于 y,则 z 的符号与 x