大数乘法之multiplication using a rectangle

Declaration

　　该文为本人原创，如需转载，请注名原作者和本站地址，但未经本人同意，严禁用于商业途径。
　　整理自我的另外一个博客，有删减。
　　原文地址：http://oss.lzu.edu.cn/blog/blog.php?do_showone/tid_355.html

Text

　　以前做的一道ACM题，这里再整理一下。那道题目的原文是[北大ACM(acm.pku.edu.cn)的第1001号题]：

Description
Problems involving the computation of exact values of very large magnitude and precision are common. For example, the computation of the national debt is a taxing experience for many computer systems.

This problem requires that you write a program to compute the exact value of Rn where R is a real number ( 0.0 < R < 99.999 ) and n is an integer such that 0 < n <= 25.

Input
The input will consist of a set of pairs of values for R and n. The R value will occupy columns 1 through 6, and the n value will be in columns 8 and 9.

Output
The output will consist of one line for each line of input giving the exact value of R^n. Leading zeros should be suppressed in the output. Insignificant trailing zeros must not be printed. Don't print the decimal point if the result is an integer.

Sample Input

95.123 12
0.4321 20
5.1234 15
6.7592 9
98.999 10
1.0100 12

Sample Output

548815620517731830194541.899025343415715973535967221869852721
.00000005148554641076956121994511276767154838481760200726351203835429763013462401
43992025569.928573701266488041146654993318703707511666295476720493953024
29448126.764121021618164430206909037173276672
90429072743629540498.107596019456651774561044010001
1.126825030131969720661201

　　这里要求我们求出一个浮点数的大整数次方，用硬件是无法直接实现的，我们必须采用软件方面的技术来解决。首先，为了保持精度，我们需要把浮点数转化为整数，之后就是大整数乘法的问题了。因此，我们这里就需要实现大整数的乘法。本来我们可以直接使用分段法来做，但是我在一本算法书里头找到了一个更妙的算法：multiplication using a rectangle。不知道翻译成什么比较好，先看看这个算法的原理再说。

　　下面这个表是用这个算法计算981　X　1234　=　1210554的实例（蓝色字体标记的是乘数和被乘数，红色背景部分是中间过程，红色字体是结果）。

（这个图用OpenOffice和dia做出来，有点麻烦哦，不知道有没有什么更方便的工具？）

　　很容易看出规律出来：

　　１、粉红色格式里头的数为对应列和行位置上数的乘积。如第一行，9=9 X 1，这样就可以求出所有粉红色格子里头的数。
　　２、把粉红色格子里头的数的不同位分别放在红色斜线的两侧，如果十位没有，那么用零补充。
　　３、从右下角开始，我们把同处于斜线一侧的数相加，结果存放到对应的左下部位的格子里头。如4=4, 5=2+0+3, 5=6+3+4+0+2 mod 10[进位为一], 0=1[进位]+3+7+2+6+0+1 mod 10,以此类推,直到左上角。

　　实际上，仔细分析以下，我们就会这种乘法和我们小学就学习过的“竖式乘法”原理是一样的，所以肯定是正确的。但是，这种描述方式给我们带来了算法设计和实现上的突破。

　　下面我们就需要把它实现了。
　　首先，我们定义数据的存储方式，乘数和被乘数直接存放到两个一維数组里头，中间过程（即粉红色背景部分）用一个三维数据来存放，存放形式如下：

　　

　　假设该三维数组的名称为Rect，则有如下关系：

Rect[3][2][1] = 4
Rect[3][2][0] = 0 Rect[3][1][1] = 2 Rect[2][2][1] = 3
...

　　很容易发现这样的规律:
　　１、计算结果的最后一位的下标最大，其下标之和为6；倒数第二位的下标之和为5，以此类推，第一位的下标之和是0。因此，下标之和的顺序刚好是结果位的位序，我们可以据此求出结果的各个位。
　　２、最大的下标之和刚好是乘数位数和被乘数位数之和减一，结果位数是乘数位数和被乘数位数之和。据此，我们可以定义结果数组和一维数组，大小是乘数和被乘数位数之和的一维数组来存放结果。

　　下面我们就根据这种思路写出我们的大数乘法的function，然后再通过调用该function，解决上面的题目——即求出大整数的幂。

Code

[1] lnmp.c -- count the product of two large numbers' multiplication

/* set ts=4 */
/*
 * FILE: lnmp.c
 * Funciton: Count the product of two large numbers' multiplication
 */

/*
 * Author: Wu Zhangjin,
 * (c) 2006-12-31 dslab.lzu.edu.cn Lanzhou University
 * License GPL Version 2
 *
 */

#include <stdio.h>        /* standard input/output functions */
#include <string.h>        /* needed by strlen function */
#define MAX_DIGIT 2000        /* define the max digit of number here*/

/**
 * mulof_ln - count the product of two large numbers' multiplication
 * @ln_arr: the address of two large numbers
 * @pofm: the address of the result array
 *
 * Description:
 * the function use the algorithm: multiplication using a rectangle
 * like this: 981 X 1234 = 1210554
 * 9 8 1
 * 1 0/9 0/8 0/1 1
 * 2 1/8 1/6 0/2 2
 * 1 2/7 2/4 0/3 3
 * 0 3/6 3/2 0/4 4
 * 5 5 4
 * you can find the the rule of it easily.
 *
 * Returns:
 * the address of the product of multiplication.
 *
 */

char* mulof_ln(char ln_arr[2][MAX_DIGIT], char pofm[2*MAX_DIGIT-1])
{
    int digit_arr[2], digit_pofm, i, j, k, m, tmp, carry;
    char td_arr[MAX_DIGIT][MAX_DIGIT][2];
    
    carry = 0;
    digit_arr[0] = strlen(ln_arr[0]);
    digit_arr[1] = strlen(ln_arr[1]);
    digit_pofm = digit_arr[0] + digit_arr[1] - 1;

    /*count the intermediate result, save them to a three-dimensional array*/
    for (i = digit_arr[0] - 1; i >= 0; i --)
    for (j = digit_arr[1] - 1; j >= 0; j --)
    {
        tmp = (ln_arr[0][i] - '0') * (ln_arr[1][j] - '0');
        td_arr[j][i][0] = tmp / 10;
        td_arr[j][i][1] = tmp % 10;
    }

    /*count the result, save it to the result array*/
    for (m = digit_pofm; m >= 0; m--)
    {
        tmp = carry;
        for (i = digit_arr[0] - 1; i >= 0; i --)
        for (j = digit_arr[1] - 1; j >= 0; j --)    
        for (k = 1; k >= 0; k --)
        if (m == i+j+k) tmp = tmp + td_arr[j][i][k];
        carry = tmp / 10;
        pofm[m] = tmp%10 + '0';
    }

    /*take out the ZERO before the result*/
    tmp = 0;
    for (i = 0; i < digit_pofm; i ++) {
        if (pofm[i] != '0') break;
        else tmp ++;
    }
    for (j = tmp, i = 0; j <= digit_pofm + 1; j ++, i ++)
        pofm[i] = pofm[j];

    return pofm;        
}

/**
 * main - the main function for testing the function: mulof_ln
 */

int main()
{
    char ln_arr[2][MAX_DIGIT];
    char pofm[2*MAX_DIGIT-1];    /* the product of the multiplication */

    printf("Please input two large numbers: \n");
    scanf("%s %s", ln_arr[0], ln_arr[1]);
        
    printf("The product of the muliplication is:\n%s\n", mulof_ln(ln_arr, pofm));
    
    return 0;
}


Demo

shell> make lnmp
cc     lnmp.c   -o lnmp
shell> ./lnmp
Please input two large numbers:
2
5
The product of the muliplication is:
10
shell> ./lnmp
Please input two large numbers:
22222222222222222222222222222222222222222222222222222222222222222222222222
5
The product of the muliplication is:
111111111111111111111111111111111111111111111111111111111111111111111111110
shell> ./lnmp
Please input two large numbers:
0000000000000000000000000000000000000000000000002
000000000000000000000000000000000000000000000055
The product of the muliplication is:
110

[2] rton.c -- count the result of r to the nth(r ^ n)
　　说明：这里并不是直接调用上面的大数乘积的function，为了更方便地操作和节省内存，进行了一定的调整。

/* set ts=4 */
/*
 * FILE: rton.c
 * Funciton: count the r to the nth
 *
 */

/*
 * Author: Wu Zhangjin,
 * (c) 2006-12-31 dslab.lzu.edu.cn Lanzhou University
 * License GPL Version 2
 *
 */

#include <stdio.h>        /* standard input/output functions */
#include <string.h>        /* needed by strlen function */
#define MAX_DIGIT 10        /* define the max digit of number here*/
#define MAX_N 25
/**
 * mulof_ln - count the product of two large numbers' multiplication
 * @m0: the address of two large numbers
 * @m1: the address of two large numbers
 * @pofm: the address of the result array
 *
 * Description:
 * the function use the algorithm: multiplication using a rectangle
 * like this: 981 X 1234 = 1210554
 * 9 8 1
 * 1 0/9 0/8 0/1 1
 * 2 1/8 1/6 0/2 2
 * 1 2/7 2/4 0/3 3
 * 0 3/6 3/2 0/4 4
 * 5 5 4
 * you can find the the rule of it easily.
 *
 * Returns:
 * the address of the product of multiplication.
 *
 * Notes:
 * the index of td_arr must be engough large, when you use it to count the power like R ^ N
 */

char* mulof_ln(char* m0, char* m1, char pofm[2*MAX_DIGIT-1])
{
    int digit_arr[2], digit_pofm, i, j, k, m, tmp, carry;
    char td_arr[MAX_DIGIT*MAX_N][MAX_DIGIT*MAX_N][2];
    
    carry = 0;
    digit_arr[0] = strlen(m0);
    digit_arr[1] = strlen(m1);
    digit_pofm = digit_arr[0] + digit_arr[1] - 1;

    /*count the intermediate result, save them to a three-dimensional array*/
    for (i = digit_arr[0] - 1; i >= 0; i --)
    for (j = digit_arr[1] - 1; j >= 0; j --)
    {    
        tmp = (m0[i] - '0') * (m1[j] - '0');
        td_arr[j][i][0] = tmp / 10;
        td_arr[j][i][1] = tmp % 10;
    }

    /*count the result, save it to the result array*/
    for (m = digit_pofm; m >= 0; m--)
    {
        tmp = carry;
        for (i = digit_arr[0] - 1; i >= 0; i --)
        for (j = digit_arr[1] - 1; j >= 0; j --)    
        for (k = 1; k >= 0; k --)
        if (m == i+j+k) tmp = tmp + td_arr[j][i][k];
        carry = tmp / 10;
        pofm[m] = tmp%10 + '0';
    }

    /*take out the ZERO before the result*/
    tmp = 0;
    for (i = 0; i < digit_pofm; i ++) {
        if (pofm[i] != '0') break;
        else tmp ++;
    }
    if (tmp > 0)
    for (j = tmp, i = 0; j <= digit_pofm + 1; j ++, i ++)
        pofm[i] = pofm[j];

    return pofm;        
}

/**
 * mypower - count the R to the nth
 * @r: the base of the power, it is a float number.
 * @n: the exponent of the power
 * @p: the address of the result
 *
 * Descript:
 * call the mulof_ln function to count the large float number R to the nth
 */

char* mypower(char* r, int n, char p[MAX_DIGIT*MAX_N])
{
    int digit_r, digit_p, i, j, pos_point, pos_point_result, tmp;

    memset(p, 0, MAX_DIGIT*MAX_N);
    digit_r = strlen(r);
    for (i = 0; i < digit_r; i ++)
        if (r[i] == '.') break;
    pos_point = i;

    /*if R is a decimal, change it to an integer firstly*/
    if (pos_point != digit_r) {
            pos_point_result = (digit_r - i - 1) * n;
            for (i = pos_point; i <= digit_r; i++)
                r[i] = r[i+1];
    }
    /*take out the ZERO before*/
    tmp = 0;
    for (i = 0; i < digit_r -1; i ++) {
        if (r[i] != '0') break;
        else tmp ++;
    }
    if (tmp > 0)
    for (j = tmp, i = 0; j <= digit_r + 1; j ++, i ++)
        r[i] = r[j];

    strcpy(p, r);
    for (i = 1; i < n; i++)
        mulof_ln(r, p, p);
    
    //normalize the output

    digit_p = strlen(p);
    if (pos_point != digit_r)
    {
        tmp = digit_p - pos_point_result;
        if (tmp >= 0) {
            j = 1;
            for (i = digit_p; i > tmp; i --) {
                if (j && p[i-1] == '0') p[i-1] = '\0';
                else {
                    p[i] = p[i-1];
                    j = 0;
                }
            }
            if (j==0) p[tmp] = '.';
        } else {
            tmp = -tmp;
            j = 1;
            for (i = digit_p; i >= 0; i --) {
                if (j && p[i-1] == '0') p[i-1] = '\0';
                else {
                    p[i+tmp] = p[i-1];
                    j = 0;
                }
            }
            for (i = 1; i <= tmp; i ++)
                p[i] = '0';
            p[0] = '.';
        }
    }

    return p;
}

/**
 * main - the main function for testing the function: mulof_ln
 */

int main()
{
    char p[MAX_DIGIT*MAX_N];
    char r[MAX_DIGIT];
    int n;

    while(scanf("%s %d", r, &n)!=EOF){
        printf("%s\n", mypower(r, n, p));
        memset(p, 0, MAX_DIGIT*MAX_N);
        memset(r, 0, MAX_DIGIT);
    }
    
    return 0;
}


Demo
　　测试数据：文件名test.txt（说明：下面的测试结果我增加了两条测试数据，非常重要的两条）

95.123 12
0.4321 20
5.1234 15
6.7592 9
98.999 10
1.0100 12
0.0001 20
3.0000 25

　　
　　测试结果：

shell>make rton
cc     rton.c   -o rton
shell> time ./rton 548815620517731830194541.899025343415715973535967221869852721
.00000005148554641076956121994511276767154838481760200726351203835429763013462401
43992025569.928573701266488041146654993318703707511666295476720493953024
29448126.764121021618164430206909037173276672
90429072743629540498.107596019456651774561044010001
1.126825030131969720661201
.00000000000000000000000000000000000000000000000000000000000000000000000000000001
847288609443

real    0m0.023s
user    0m0.020s
sys     0m0.000s

Extention

　　感觉上面的代码写得比较乱，尤其是变量声明，内存使用都没有得到很好的控制。有时间得再修改一下。也欢迎你提出批评和意见。

Download
　　

文件:	lnmp-0.1.tar.gz [大数乘法]
大小:	1KB
下载:	下载

　　

文件:	rton-0.1.tar.gz [求大数的幂]
大小:	1KB
下载:	下载