MAPM代码：任意有效位数的乘法实现解析（1）基本乘法

先看看MAPM代码中附带的ARTICLE.PDF中的描述:

乘法通常是一个O(N^2)次运算，对于10000位数字的乘法要进行一亿次的运算。但是利用分部运算法（divide-and-conquer algorithm）可以达到更快的乘法，对于N位数字的时间复杂度大约是O(N^1.585)，也就是说对于10000位数字的乘法只需要2,188,000次运算。而真正快速的乘法（Really Fast Multiplication）试用了FFT(快速傅里叶变换Fast Fourier Transform)。FFT乘法的时间复杂度仅仅是O(N*log2(N))，这比一般的N^2和分布运算的N^1.585要快得多。

FFT的基本思路是这样的。首先，对两个数都做正向傅立叶变换，数字的位数作为输入FFT的样本。这样在频域（frequency domain）生成两套系数，每一套系数都包含和原始数字中的位数一样多的样本。（同样，每个系数是一个复数）第二步，将这两套系数相乘，即假设两个数a和b相乘，那么执行Ca0×Cb0，Ca1×Cb1，等等，Can、Cbn是a和b的频域系数。这个在频域逐对（pair-wise）相乘等同于“时间”或采样域的卷积（convolution）。这样生成的是正确的结果，因为当将N位数数字相乘的时候，你事实上就是在它们的位上逐个进行卷积。第三步，将结果进行反变换。最后，将逆FFT的实数部分转换成整数，释放计算过程中使用的数。MAPM中使用的傅立叶变换来自Takuya Ooura。这个FFT是快速、可移植的，并可自由分发。

通过一个实际情况比较一下三种算法，计算1,000,000位数字乘法所消耗的时间：

FFT Method : 40 seconds
Divide-and-Conquer : 1 hr, 50 min
Ordinary N2 : 23.9 days*

*projected估计的!

在MAPM中，这三种算法都被用到了。对于输入的普通大小的数，如43×87，使用O(N^2)算法。然后就是使用FFT了。不过基于FFT的算法会在某个点上出现浮点运算溢出，此时先使用分部运算法解决这个问题，然后再用FFT法完成计算。

我的这个帖子讲的是O(N^2)算法，分部运算和基于FFT的算法将写在另一篇帖子里。

MAPM的基本数据结构定义如下：

typedef struct  {
    UCHAR    *m_apm_data;
    long    m_apm_id;
    int     m_apm_refcount;       /* <- used only by C++ MAPM class */
    int    m_apm_malloclength;
    int    m_apm_datalength;
    int    m_apm_exponent;
    int    m_apm_sign;
} M_APM_struct;

typedef M_APM_struct *M_APM;

其中的成员定义如下：

名称	定义
m_apm_data	以10进制的方式按字节存储数值。例如1750（0x6D6）的存储方式是0x11 0x32，而不是0x06 0xD6。
m_apm_id	常数，定义为：#define    M_APM_IDENT 0x6BCC9AE5
m_apm_refcount	for MAPM C++ class
m_apm_malloclength	第一次初始化为80字节。
m_apm_datalength	指significant digits有效数
m_apm_exponent	指数，又称阶码。
m_apm_sign	正数为0x00000001，负数为0xFFFFFFFF，零数时为0x00000000。

在Windows x86(little endian)系统中，典型的内存结构是：

0x003D3260  b8 32 3d 00 e5 9a cc 6b 01 00 00 00 50 00 00 00  .2=....k....P...
0x003D3270  09 00 00 00 05 00 00 00 01 00 00 00 fd fd fd fd  ................

其中m_apm_data指向的内存信息是，此时通过m_apm_set_string(a, "99831.5355")给变量a赋值：

0x003D32B8  63 53 0f 23 32 00 cd cd cd cd cd cd cd cd cd cd  cS.#2...........
0x003D32C8  cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd  ................
0x003D32D8  cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd  ................
0x003D32E8  cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd  ................
0x003D32F8  cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd cd  ................
0x003D3308  cd cd cd cd

从中我们可以看出，a的阶码是5，在内存中的表示为：0x63 53 0f 23 32，对应的10进制值为：99 83 15 35 50 00。由此我们可以看出，MAPM将有效位全部作为小数部分，与阶码一起定义数值。这样在计算乘法的时候，只需要把小数部分当作整型数值来对待就可以了。

下面是乘法函数：

void    m_apm_multiply(M_APM r, M_APM a, M_APM b)
{
int    ai, itmp, sign, nexp, ii, jj, indexa, indexb, index0, numdigits;
UCHAR   *cp, *cpr, *cp_div, *cp_rem;
void    *vp;

sign = a->m_apm_sign * b->m_apm_sign;
nexp = a->m_apm_exponent + b->m_apm_exponent;

if (sign == 0)      /* one number is zero, result is zero */
  {
   M_set_to_zero(r);
   return;
  }

numdigits = a->m_apm_datalength + b->m_apm_datalength;
indexa = (a->m_apm_datalength + 1) >> 1;
indexb = (b->m_apm_datalength + 1) >> 1;

/*
 *    If we are multiplying 2 'big' numbers, use the fast algorithm.
 *
 *    This is a **very** approx break even point between this algorithm
 *      and the FFT multiply. Note that different CPU's, operating systems,
 *      and compiler's may yield a different break even point. This point
 *      (~96 decimal digits) is how the test came out on the author's system.
 */

if (indexa >= 48 && indexb >= 48)
  {
   M_fast_multiply(r, a, b);
   return;
  }

ii = (numdigits + 1) >> 1;     /* required size of result, in bytes */

if (ii > r->m_apm_malloclength)
  {
   if ((vp = MAPM_REALLOC(r->m_apm_data, (ii + 32))) == NULL)
     {
      /* fatal, this does not return */

      M_apm_log_error_msg(M_APM_FATAL, "\'m_apm_multiply\', Out of memory");
     }

   r->m_apm_malloclength = ii + 28;
   r->m_apm_data = (UCHAR *)vp;
  }

M_get_div_rem_addr(&cp_div, &cp_rem);

index0 = indexa + indexb;
cp = r->m_apm_data;
memset(cp, 0, index0);
ii = indexa;

while (TRUE)
  {
   index0--;
   cpr = cp + index0;
   jj  = indexb;
   ai  = (int)a->m_apm_data[--ii];

   while (TRUE)
     {
      itmp = ai * b->m_apm_data[--jj];

      *(cpr-1) += cp_div[itmp];
      *cpr     += cp_rem[itmp];

      if (*cpr >= 100)
        {
         *cpr     -= 100;
         *(cpr-1) += 1;
    }

      cpr--;

      if (*cpr >= 100)
        {
         *cpr     -= 100;
         *(cpr-1) += 1;
    }

      if (jj == 0)
        break;
     }

   if (ii == 0)
     break;
  }

r->m_apm_sign       = sign;
r->m_apm_exponent   = nexp;
r->m_apm_datalength = numdigits;

M_apm_normalize(r);
}

第17、18行中计算出的indexa、indexb分别是乘法左右操作数的有效位数的一半个数（注意，超出部分进1，即从9位得出5位），也就是存储数值所需字节数，因为一个字节存储两个字符。在MAPM中，使用普通字符（CHAR类型）存储每一位数字。这用于判断是否是大数相乘，此处在代码编写者的系统中用96来判断大数位数，在不同的系统上这个数字是有很大差别的，所以仅仅是一个经验数值。

第35行中计算出计算结果所需的字节数，算法同indexa、indexb的计算方法。

第50行中获取两个变量div、rem的地址，分别指向静态变量M_mul_div和M_mul_rem，分别存储10000个字节数值。另外还有静态变量M_mul_div_10和M_mul_rem_10，分别存储100个字节的数值。这四个变量都是计算过程中的工具数，M_mul_div头100个字节数值为0，以后每递增100个字节数值增1；M_mul_rem的每一百个字节均递增存储0～99的数值；M_mul_div_10头10个字节数值为0，以后每递增10个字节数值增1；M_mul_rem_10的每十个字节均递增存储0～9的数值。这其中的作用是，因为左右计算数和最终结果计算数中每个字节存储的是0～99的数值，也就是对于百位数相乘，M_uml_div存储的是百位以上的结果数值，例如35×50=1750，在M_mul_div的第1700～1799字节的数值也正好是17，而在M_mul_rem的第1700～1799的数值是0～99（第1750字节的数值正好是1750）。这样，MAPM将乘法替换为每个字节与每个字节的乘法，依次按字节移动内存位置，通过加法得到字节相乘的结果。注意，这些计算都相当于是人工10进制的计算过程，参见成员m_apm_data的定义。另外，如果不使用工具数，而是从计算结果中通过计算（需要用除法和求余，例如1750/100=17，1750%100=50）得到累加数，会耗费额外的时间。这就是一个以空间换效率的方法。

在乘法计算过程中，通过分别对两个操作数逐字节的运算结果的累加，得到最终的数值结果。