2007-10-29

大数乘法巧解

关键字: 转贴
下面这个表是用这个算法计算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, <zhangjinw@gmail.com>
 * (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, <zhangjinw@gmail.com>
 * (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
评论
发表评论

您还没有登录，请登录后发表评论

younker
搜索本博客
最近加入圈子
存档
最新评论