大整数乘法算法

一 转换为二进制求，推导出的公式适合十进制计算

设X和Y都是n位的二进制整数，现在要计算它们的乘积XY。我们可以用小学所学的方法来设计一个计算乘积XY的算法，但是这样做计算步骤太多，显得效率较低。如果将每2个1位数的乘法或加法看作一步运算，那么这种方法要作O(n²)步运算才能求出乘积XY。下面我们用分治法来设计一个更有效的大整数乘积算法。

　

图6-3 大整数X和Y的分段　

我们将n位的二进制整数X和Y各分为2段，每段的长为n/2位(为简单起见，假设n是2的幂)，如图6-3所示。

由此，X=A2^n/2+B ，Y=C2^n/2+D。这样，X和Y的乘积为：

XY=(A2^n/2+B)(C2^n/2+D)=AC2ⁿ+(AD+CB)2^n/2+BD             (1)

如果按式(1)计算XY，则我们必须进行4次n/2位整数的乘法(AC，AD，BC和BD)，
以及3次不超过n位的整数加法(分别对应于式(1)中的加号)，此外还要做2次移位(分别对应于式(1)中乘2ⁿ和乘2^n/2)。所有这些加法和移位共用O(n)步运算。设T(n)是2个n位整数相乘所需的运算总数，则由式(1)，我们有:

                             (2)

由此可得T(n)=O(n²)。因此，用(1)式来计算X和Y的乘积并不比小学生的方法更有效。要想改进算法的计算复杂性，必须减少乘法次数。为此我们把XY写成另一种形式:

XY=AC2ⁿ+[(A-B)(D-C)+AC+BD]2^n/2+BD                     (3)

虽然，式(3)看起来比式(1)复杂些，但它仅需做3次n/2位整数的乘法(AC，BD和(A-B)(D-C))，6次加、减法和2次移位。由此可得:

                                 (4)

用解递归方程的套用公式法马上可得其解为T(n)=O(n^log3)=O(n^1.59)。利用式(3)，并考虑到X和Y的符号对结果的影响，我们给出大整数相乘的完整算法MULT如下:

function MULT(X，Y，n); {X和Y为2个小于2n的整数，返回结果为X和Y的乘积XY}
begin
  S:=SIGN(X)*SIGN(Y); {S为X和Y的符号乘积}
  X:=ABS(X);
  Y:=ABS(Y); {X和Y分别取绝对值}
  if n=1 then
     if (X=1)and(Y=1) then return(S)
                      else return(0)
         else begin
                A:=X的左边n/2位;
                B:=X的右边n/2位;
                C:=Y的左边n/2位;
                D:=Y的右边n/2位;
                ml:=MULT(A,C,n/2);
                m2:=MULT(A-B,D-C,n/2);
                m3:=MULT(B,D,n/2);
                S:=S*(m1*2n+(m1+m2+m3)*2n/2+m3);
                return(S);
              end;
end;

上述二进制大整数乘法同样可应用于十进制大整数的乘法以提高乘法的效率减少乘法次数。下面的例子演示了算法的计算过程。

设X=314l,Y=5327，用上述算法计算XY的计算过程可列表如下，其中带^'号的数值是在计算完成AC，BD，和(A-B)(D-C)之后才填入的。

---

X=3141        A=31       B=41        A-B=-10

Y=5327        C=53       D=27        D-C=-26

           AC=(1643)^'

           BD=(1107)^'

          (A-B)(D-C)=(260)^'

XY=(1643)^'10⁴+[(1643)^'+(260)^'+(1107)^']10²+(1107)^'

  =(16732107)^'

---

A=31        A₁=3       B₁=1        A₁-B₁=2

C=53        C₁=5       D₁=3        D₁-C₁=-2

           A₁C₁=15     B₁D₁=3     (A₁-B₁)(D₁-C₁)=-4

AC=1500+(15+3-4)10+3=1643

---

B=41        A₂=4       B₂=1        A₂-B₂=3

D=27        C₂=2       D₂=7        D₂-C₂=5

           A₂C₂=8     B₂D₂=7     (A₂-B₂)(D₂-C₂)=15

BD=800+(8+7+15)10+7=1107

---

|A-B|=10        A₃=1       B₃=0        A₃-B₃=1

|D-C|=26        C₃=2       D₃=6        D₃-C₃=4

           A₃C₃=2     B₃D₃=0     (A₃-B₃)(D₃-C₃)=4

(A-B)(D-C)=200+(2+0+4)10+0=260

二、另外的思路：

        比较容易想到的是做多位数乘法时列竖式进行计算的方法，只要写出模拟这一过程的程序，就能实现任意大整数的乘法运算。经过查阅资料，找到一种更易于编程的方法，即“列表法”。 

下面先介绍“列表法”： 
例如当计算8765 x 234时，把乘数与被乘数照如下列出，见表1：                         

                        
把表1中的数按图示斜线分组（横纵坐标和相等的数分为一组），把每组数的累加起来所得的和记在表格下方，见表 2：

        从最低位的 20 开始，保留个位数字“0”，把个位以外的数“2”进到前一位；把次低位的 39 加上低位进上来的 2 得 41，保留个位数字“1”，把“4”进到前一位；以此类推，直至最高位的 16，16 加上低位进上来的4得 20，保留“0”，把2进到最高位，得乘积答数 2051010。                      

根据以上思路就可以编写C 程序了，再经分析可得： 
1、一个m 位的整数与一个 n 位的整数相乘，乘积为m+n-1 位或m+n 位。 
2、程序中，用三个字符数组分别存储乘数、被乘数与乘积。由第 1 点分析知，存放乘积的字符数组
的长度应不小于存放乘数与被乘数的两个数组的长度之和。 
3、可以把第二步“计算填表”与第三四步“累加进位”放在一起完成，可以节省存储表格 2所需的空间。 
4、程序关键部分是两层循环，内层循环累计一组数的和，外层循环处理保留的数字与进位。 

编写的程序如下： 
#define MAXLENGTH 1000 
#include <stdio.h> 
#include <string.h> 

void compute(char *a, char *b, char *c); 

void main(void) 
{ 
    char a[MAXLENGTH], b[MAXLENGTH], c[MAXLENGTH * 2]; 

    puts("Input multiplier :"); 
    gets(a); 
    puts("Input multiplicand :"); 
    gets(b); 

    compute(a, b, c); 

    puts("Answer :"); 
    puts(c); 
    getchar(); 
} 

void compute(char *a, char *b, char *c) 
{ 
    int i, j, m, n; 
    long sum, carry; 

    m = strlen(a) - 1; 
    n = strlen(b) - 1; 
    for (i = m; i >= 0; i--) 
        a[i] -= '0'; 
    for (i = n; i >= 0; i--) 
        b[i] -= '0'; 
    c[m + n + 2] = '\0'; 

    carry = 0; 
    for (i = m + n; i >= 0; i--) /* i 为坐标和 */ 
    {
        sum = carry; 

        if ((j = i - m) < 0) 
            j = 0; 
        for ( ; j<=i && j<=n; j++) /* j 为纵坐标 */ 
            sum += a[i-j] * b[j]; /* 累计一组数的和 */ 

        c[i + 1] = sum % 10 + '0'; /* 算出保留的数字 */ 
        carry = sum / 10; /* 算出进位 */ 

    } 

    if ((c[0] = carry+'0') == '0') /* if no carry, */ 
        c[0] = '\040'; /* c[0] equals to space */ 
} 

效率分析：用以上算法计算 m位整数乘以n 位整数，需要先进行 m x n次乘法运算，再进行约 m + n次加法运算和 m + n次取模运算(实为整数除法)。把这个程序稍加修改，让它自己产生乘数与被乘数，然后计算随机的 7200位整数互乘，在Cyrix 6x86 pr166机器的纯DOS方式下耗时 7秒(用Borland C3.1编译)。

经过改进，此算法效率可以提高约9 倍。 
注意到以下事实：8216547 x 96785 将两数从个位起，每 3位分为节，列出乘法表，将斜线间的数字相加；

8   216   547 
     96     785

将表中最后一行进行如下处理：从个位数开始，每一个方格里只保留三位数字，超出 1000 的部
分进位到前一个方格里； 

所以8216547 x 96785 = 795238501395 

也就是说我们在计算生成这个二维表时，不必一位一位地乘，而可以三位三位地乘；在累加时也是满1000进位。这样，我们在计算 m位整数乘以 n位整数，只需要进行 m x n / 9次乘法运算，再进行约(m + n) / 3次加法运算和(m + n) /3 次取模运算。总体看来，效率约是前一种算法的 9倍。 
有人可能会想：既然能够三位三位地乘，为什么不4位 4位甚至5位5位地乘呢？那不是可以提高 16 乃至 25 倍效率吗？听我解来：本算法在累加表中斜线间的数字时，如果用无符号长整数(范围 0至~4294967295)作为累加变量，在最不利的情况下(两个乘数的所有数字均是 9)，能够累加约4294967295/(999*999)=4300 次，也就是能够准确计算任意两个均不超过 12900(每次累加的结果"值"三位，故 4300*3=12900)位的整数相乘。如果 4 位 4 位地乘，在最不利的情况下，能够累加约4294967295/(9999*9999)=43 次，仅能够确保任意两个不超过 172 位的整数相乘，没有什么实用价值，更不要说5位了。 

请看改进后的算法的实例程序： 
该程序随机产生两个72xx位的整数，把乘数与积保存在 result.txt中。在Borland C++ 3.1 中用 
BCC -3 -O2 -G -mh -Z -f287 -pr -T- dashu.cpp 编译生成的exe文件在Cyrix 6x86 pr166的机器上运行耗时0.82 秒。 
程序 2 清单： 
#include<conio.h> 
#include<string.h> 
#include<stdio.h> 
#include<stdlib.h> 
#include<time.h> 
#define N 7200 //作 72xx 位的整数乘法 
int max(int,int,int); 
int initarray(int a[]); 
void write(int a[],int l); 
FILE *fp; 

void main() 
{ 
    int a[5000]={0},b[5000]={0},k[10001]={0}; //声明存放乘数、被乘数与积的数组 
    clock_t start, end; //声明用于计时的变量 
    unsigned long c,d,e; //声明作累加用的无符号长整数变量 
    int i,j,la,lb,ma,mi,p,q,t; //声明其它变量 
    randomize(); //初始化随机数 

    la=initarray(a); //产生被乘数，并返回其长度 
    lb=initarray(b); //产生乘数，并返回其长度 

    if(la<lb) //如果被乘数长度小于乘数，则交换被乘数与乘数 
    {
        p=(lb>la)?lb:la; 
        for (q=0;q<p;q++) //交换被乘数与乘数 
        t=a[q],a[q]=b[q],b[q]=t; 
        t=la,la=lb,lb=t; //交换被乘数的长度与乘数的长度 
    } 

    start = clock();//开始计时
    c=d=0; //清空累加变量，其中 C 用于累加斜线间的数，d 用作进位标志
    for(i=la+lb-2;i>=0;i--) //累加斜线间的数，i 为横纵坐标之和
    {
        c=d; //将前一位的进位标志存入累加变量 c
        ma=max(0,i-la+1,i-lb+1); //求累加的下限
        mi=(i>la-1)?(la-1):i; //求累加的上限
        for(j=ma;j<=mi;j++) //计算出横纵坐标之和为 i 的单元内的数,并累加到 C 中
            c+=(long)a[j]*b[i-j];
        d=c/1000; //求进位标志
        if(c>999)
            c%=1000; //取 c 的末三位
        k[i]=c; //保存至表示乘积的数组 k[]
    }
    e=k[0]+1000*d; //求出乘积的最高位
    end = clock();//停止计时
    fp = fopen("result.txt", "w+"); //保存结果到 result.txt
    printf("\nThe elapsed time was: %3.4f\n", (end - start) / CLK_TCK);
    //打印消耗的时间
    fprintf(fp,"%d",a[0]); //打印被乘数最高位
    write(a,la); //打印被乘数其他位
    fprintf(fp,"%d",b[0]); //打印乘数最高位
    write(b,lb); //打印乘数其他位
    fprintf(fp,"%ld",e); //打印乘积最高位
    write(k,la+lb-1); //打印乘积其他位
    fclose(fp);
}

max(int a,int b,int c) 
{ 
    int d; 
    d=(a>b)?a:b; 
    return (d>c)?d:c; 
} 

int initarray(int a[]) 
{ 
    int q,p,i; 
    q=N+random(100); 
    if(q%3==0) 
        p=q/3; 
    else 
        p=q/3+1; 
    
    for(i=0;i<p;i++) 
        a[i]=random(1000); 
    
    if(q%3==0) 
        a[0]=100+random(900); 
    if(q%3==2) 
        a[0]=10+random(90); 
    if(q%3==1) 
        a[0]=1+random(9); 
    
    return p; 
} 

void write(int a[],int l) 
{ 
    int i; 
    char string[10]; 
    for(i=1;i<l;i++)
    { 
        itoa(a[i],string,10); 
        if (strlen(string)==1) 
            fprintf(fp,"00"); 
        if (strlen(string)==2) 
            fprintf(fp,"0"); 
        fprintf(fp,"%s",string); 
        if((i+1)%25==0) 
            fprintf(fp,"\n"); 
    } 
    fprintf(fp,"\n"); 
    fprintf(fp,"\n"); 
}