MAPM代码：任意有效位数的除法实现解析

先看看MAPM代码中包含的ARTICLE.PDF，对除法设计的描述是这样的：

使用了两个除法函数。对于少于250位的数字，我使用Knuth在The Art of Computer Programming, Volume 2一书中写的算法，稍有一些改动。

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_divide(M_APM rr, int places, M_APM aa, M_APM bb)
{
M_APM   tmp0, tmp1;
int     sn, nexp, dplaces;

sn = aa->m_apm_sign * bb->m_apm_sign;

if (sn == 0)                  /* one number is zero, result is zero */
  {
   if (bb->m_apm_sign == 0)
     {
      M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_divide\', Divide by 0");
     }

   M_set_to_zero(rr);
   return;
  }

/*
 *    Use the original 'Knuth' method for smaller divides. On the
 *    author's system, this was the *approx* break even point before
 *    the reciprocal method used below became faster.
 */

if (places < 250)
  {
   M_apm_sdivide(rr, places, aa, bb);
   return;
  }

/* mimic the decimal place behavior of the original divide */

nexp = aa->m_apm_exponent - bb->m_apm_exponent;

if (nexp > 0)
  dplaces = nexp + places;
else
  dplaces = places;

tmp0 = M_get_stack_var();
tmp1 = M_get_stack_var();

m_apm_reciprocal(tmp0, (dplaces + 8), bb);
m_apm_multiply(tmp1, tmp0, aa);
m_apm_round(rr, dplaces, tmp1);

M_restore_stack(2);
}
/****************************************************************************/
void    m_apm_reciprocal(M_APM rr, int places, M_APM aa)
{
M_APM   last_x, guess, tmpN, tmp1, tmp2;
char    sbuf[32];
int    ii, bflag, dplaces, nexp, tolerance;

if (aa->m_apm_sign == 0)
  {
   M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_reciprocal\', Input = 0");

   M_set_to_zero(rr);
   return;
  }

last_x = M_get_stack_var();
guess  = M_get_stack_var();
tmpN   = M_get_stack_var();
tmp1   = M_get_stack_var();
tmp2   = M_get_stack_var();

m_apm_absolute_value(tmpN, aa);

/*
    normalize the input number (make the exponent 0) so
    the 'guess' below will not over/under flow on large
    magnitude exponents.
*/

nexp = aa->m_apm_exponent;
tmpN->m_apm_exponent -= nexp;

m_apm_to_string(sbuf, 15, tmpN);
m_apm_set_double(guess, (1.0 / atof(sbuf)));

tolerance = places + 4;
dplaces   = places + 16;
bflag     = FALSE;

m_apm_negate(last_x, MM_Ten);

/*   Use the following iteration to calculate the reciprocal :


         X     =  X  *  [ 2 - N * X ]
          n+1
*/

ii = 0;

while (TRUE)
  {
   m_apm_multiply(tmp1, tmpN, guess);
   m_apm_subtract(tmp2, MM_Two, tmp1);
   m_apm_multiply(tmp1, tmp2, guess);

   if (bflag)
     break;

   m_apm_round(guess, dplaces, tmp1);

   /* force at least 2 iterations so 'last_x' has valid data */

   if (ii != 0)
     {
      m_apm_subtract(tmp2, guess, last_x);

      if (tmp2->m_apm_sign == 0)
        break;

      /*
       *   if we are within a factor of 4 on the error term,
       *   we will be accurate enough after the *next* iteration
       *   is complete.
       */

      if ((-4 * tmp2->m_apm_exponent) > tolerance)
        bflag = TRUE;
     }

   m_apm_copy(last_x, guess);
   ii++;
  }

m_apm_round(rr, places, tmp1);
rr->m_apm_exponent -= nexp;
rr->m_apm_sign = aa->m_apm_sign;
M_restore_stack(5);
}
/****************************************************************************/

其中调用了m_apm_reciprocal(), m_apm_multiply()和m_apm_round()等函数。其中m_apm_reciprocal()用于取倒数，m_apm_round()用于

void    m_apm_round(M_APM btmp, int places, M_APM atmp)
{
int    ii;

if (M_util_firsttime)
  {
   M_util_firsttime = FALSE;

   M_work_0_5 = m_apm_init();
   m_apm_set_string(M_work_0_5, "5");
  }

ii = places + 1;

if (atmp->m_apm_datalength <= ii)
  {
   m_apm_copy(btmp,atmp);
   return;
  }

M_work_0_5->m_apm_exponent = atmp->m_apm_exponent - ii;

if (atmp->m_apm_sign > 0)
  m_apm_add(btmp, atmp, M_work_0_5);
else
  m_apm_subtract(btmp, atmp, M_work_0_5);

btmp->m_apm_datalength = ii;
M_apm_normalize(btmp);
}