大数乘法的C代码实现

Posted 2020-09-26 veli

在C语言中，宽度最大的无符号整数类型是unsigned long long, 占8个字节。那么，如果整数超过8个字节，如何进行大数乘法呢？ 例如：

$ python
Python 2.7.6 (default, Oct 26 2016, 20:32:47) 
...<snip>....
>>> a = 0x123456781234567812345678
>>> b = 0x876543211234567887654321
>>> print "a * b = 0x%x" % (a * b)
a * b = 0x9a0cd057ba4c159a33a669f0a522711984e32bd70b88d78

用C语言实现大数乘法，基本思路是采用分而治之的策略，难点就是进位处理相对复杂一些。本文尝试给出C代码实现，并使用Python脚本验证计算结果。

1. foo.c

  1 #include <stdio.h>
  2 #include <stdlib.h>
  3 #include <string.h>
  4 
  5 typedef unsigned char           byte;   /* 1 byte */
  6 typedef unsigned short          word;   /* 2 bytes */
  7 typedef unsigned int            dword;  /* 4 bytes */
  8 typedef unsigned long long      qword;  /* 8 bytes */
  9 
 10 typedef struct big_number_s {
 11         dword   *data;
 12         dword   size;
 13 } big_number_t;
 14 
 15 static void
 16 dump(char *tag, big_number_t *p)
 17 {
 18         if (p == NULL)
 19                 return;
 20 
 21         printf("%s : data=%p : size=%d:\t", tag, p, p->size);
 22         for (dword i = 0; i < p->size; i++)
 23                 printf("0x%08x ", (p->data)[i]);
 24         printf("\n");
 25 }
 26 
 27 /*
 28  * Add 64-bit number (8 bytes) to a[] whose element is 32-bit int (4 bytes)
 29  *
 30  * e.g.
 31  *      a[] = {0x12345678,0x87654321,0x0}; n = 3;
 32  *      n64 =  0xffffffff12345678
 33  *
 34  *      The whole process of add64() looks like:
 35  *
 36  *             0x12345678 0x87654321 0x00000000
 37  *          +  0x12345678 0xffffffff
 38  *          -----------------------------------
 39  *          =  0x2468acf0 0x87654321 0x00000000
 40  *          +             0xffffffff
 41  *          -----------------------------------
 42  *          =  0x2468acf0 0x87654320 0x00000001
 43  *
 44  *      Finally,
 45  *      a[] = {0x2468acf0,0x87654320,0x00000001}
 46  */
 47 static void
 48 add64(dword a[], dword n, qword n64)
 49 {
 50         dword carry = 0;
 51 
 52         carry = n64 & 0xFFFFFFFF; /* low 32 bits of n64 */
 53         for (dword i = 0; i < n; i++) {
 54                 if (carry == 0x0)
 55                         break;
 56 
 57                 qword t = (qword)a[i] + (qword)carry;
 58                 a[i] = t & 0xFFFFFFFF;
 59                 carry = (word)(t >> 32); /* next carry */
 60         }
 61 
 62         carry = (dword)(n64 >> 32); /* high 32 bits of n64 */
 63         for (dword i = 1; i < n; i++) {
 64                 if (carry == 0x0)
 65                         break;
 66 
 67                 qword t = (qword)a[i] + (qword)carry;
 68                 a[i] = t & 0xFFFFFFFF;
 69                 carry = (word)(t >> 32); /* next carry */
 70         }
 71 }
 72 
 73 static big_number_t *
 74 big_number_mul(big_number_t *a, big_number_t *b)
 75 {
 76         big_number_t *c = (big_number_t *)malloc(sizeof(big_number_t));
 77         if (c == NULL) /* malloc error */
 78                 return NULL;
 79 
 80         c->size = a->size + b->size;
 81         c->data = (dword *)malloc(sizeof(dword) * c->size);
 82         if (c->data == NULL) /* malloc error */
 83                 return NULL;
 84 
 85         memset(c->data, 0, sizeof(dword) * c->size);
 86 
 87         dword *adp = a->data;
 88         dword *bdp = b->data;
 89         dword *cdp = c->data;
 90         for (dword i = 0; i < a->size; i++) {
 91                 if (adp[i] == 0x0)
 92                         continue;
 93 
 94                 for (dword j = 0; j < b->size; j++) {
 95                         if (bdp[j] == 0x0)
 96                                 continue;
 97 
 98                         qword n64 = (qword)adp[i] * (qword)bdp[j];
 99                         dword *dst = cdp + i + j;
100                         add64(dst, c->size - (i + j), n64);
101                 }
102         }
103 
104         return c;
105 }
106 
107 static void
108 free_big_number(big_number_t *p)
109 {
110         if (p == NULL)
111                 return;
112 
113         if (p->data != NULL)
114                 free(p->data);
115 
116         free(p);
117 }
118 
119 int
120 main(int argc, char *argv[])
121 {
122         dword a_data[] = {0x12345678, 0x9abcdef0, 0xffffffff, 0x9abcdefa, 0x0};
123         dword b_data[] = {0xfedcba98, 0x76543210, 0x76543210, 0xfedcba98, 0x0};
124 
125         big_number_t a;
126         a.data = (dword *)a_data;
127         a.size = sizeof(a_data) / sizeof(dword);
128 
129         big_number_t b;
130         b.data = (dword *)b_data;
131         b.size = sizeof(b_data) / sizeof(dword);
132 
133         dump("BigNumber A", &a);
134         dump("BigNumber B", &b);
135         big_number_t *c = big_number_mul(&a, &b);
136         dump("  C = A * B", c);
137         free_big_number(c);
138 
139         return 0;
140 }

2. bar.py

 1 #!/usr/bin/python
 2 
 3 import sys
 4 
 5 def str2hex(s):
 6     l = s.split(‘ ‘)
 7 
 8     i = len(l)
 9     out = ""
10     while i > 0:
11         i -= 1
12         e = l[i]
13         if e.startswith("0x"):
14             e = e[2:]
15         out += e
16 
17     out = "0x%s" % out
18     n = eval("%s * %d" % (out, 0x1))
19     return n
20 
21 def main(argc, argv):
22     if argc != 4:
23         sys.stderr.write("Usage: %s <a> <b> <c>\n" % argv[0])
24         return -1
25 
26     a = argv[1]
27     b = argv[2]
28     c = argv[3]
29     ax = str2hex(a)
30     bx = str2hex(b)
31 
32     axbx = ax * bx
33 
34     s_hex = "%x" % axbx
35     if s_hex.startswith("0x"):
36         s_hex = s_hex[2:]
37 
38     n = len(s_hex)
39     m = n % 8
40     if m != 0:
41         s_hex = ‘0‘ * (8 - m) + s_hex
42         n += (8 - m)
43     i = n
44     l = []
45     while i >= 8:
46         l.append(‘0x‘ + s_hex[i-8:i])
47         i -= 8
48     got = "%s" % ‘ ‘.join(l)
49 
50     print "0x%x * 0x%x = " % (ax, bx)
51     print "0x%x = " % axbx
52     print "%s <-- got" % got
53     print "%s <-- exp" % c
54     if got != c:
55         print "\nFalse"
56         return 1
57     else:
58         print "\nTrue"
59 
60     return 0
61 
62 if __name__ == ‘__main__‘:
63     argv = sys.argv
64     argc = len(argv)
65     sys.exit(main(argc, argv))

3. Makefile

CC        = gcc
CFLAGS        = -g -Wall -m32 -std=c99

TARGETS        = foo bar

all: $(TARGETS)

foo: foo.c
    $(CC) $(CFLAGS) -o [email protected] $<

bar: bar.py
    cp $< [email protected] && chmod +x [email protected]

clean:
    rm -f *.o
clobber: clean
    rm -f $(TARGETS)
cl: clobber

4. 编译并测试

$ make
gcc -g -Wall -m32 -std=c99 -o foo foo.c
cp bar.py bar && chmod +x bar
$ ./foo
BigNumber A : data=0xbfffc058 : size=5: 0x12345678 0x9abcdef0 0xffffffff 0x9abcdefa 0x00000000
BigNumber B : data=0xbfffc060 : size=5: 0xfedcba98 0x76543210 0x76543210 0xfedcba98 0x00000000
  C = A * B : data=0x8234008 : size=10: 0x35068740 0xee07360a 0x053bd8c9 0x2895f6cd 0xb973e57e 0x4e6cfe66 0x0b60b60b 0x9a0cd056 0x00000000 0x00000000
$ A="0x12345678 0x9abcdef0 0xffffffff 0x9abcdefa 0x00000000"
$ B="0xfedcba98 0x76543210 0x76543210 0xfedcba98 0x00000000"
$ C="0x35068740 0xee07360a 0x053bd8c9 0x2895f6cd 0xb973e57e 0x4e6cfe66 0x0b60b60b 0x9a0cd056"
$ ./bar "$A" "$B" "$C"
0x9abcdefaffffffff9abcdef012345678 * 0xfedcba987654321076543210fedcba98 =
0x9a0cd0560b60b60b4e6cfe66b973e57e2895f6cd053bd8c9ee07360a35068740 =
0x35068740 0xee07360a 0x053bd8c9 0x2895f6cd 0xb973e57e 0x4e6cfe66 0x0b60b60b 0x9a0cd056 <-- got
0x35068740 0xee07360a 0x053bd8c9 0x2895f6cd 0xb973e57e 0x4e6cfe66 0x0b60b60b 0x9a0cd056 <-- exp

True
$

结束语： TBD