1 #include "os.h"
2 #include <mp.h>
3 #include "dat.h"
4 
5 #define LO(x) ((x) & (((mpdigit)1<<(Dbits/2))-1))
6 #define HI(x) ((x) >> (Dbits/2))
7 
8 static void
9 mpdigmul(mpdigit a, mpdigit b, mpdigit *p)
10 {
11 	mpdigit x, ah, al, bh, bl, p1, p2, p3, p4;
12 	int carry;
13 
14 	/* half digits */
15 	ah = HI(a);
16 	al = LO(a);
17 	bh = HI(b);
18 	bl = LO(b);
19 
20 	/* partial products */
21 	p1 = ah*bl;
22 	p2 = bh*al;
23 	p3 = bl*al;
24 	p4 = ah*bh;
25 
26 	/* p = ((p1+p2)<<(Dbits/2)) + (p4<<Dbits) + p3 */
27 	carry = 0;
28 	x = p1<<(Dbits/2);
29 	p3 += x;
30 	if(p3 < x)
31 		carry++;
32 	x = p2<<(Dbits/2);
33 	p3 += x;
34 	if(p3 < x)
35 		carry++;
36 	p4 += carry + HI(p1) + HI(p2);	/* can't carry out of the high digit */
37 	p[0] = p3;
38 	p[1] = p4;
39 }
40 
41 /* prereq: p must have room for n+1 digits */
42 void
43 mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
44 {
45 	int i;
46 	mpdigit carry, x, y, part[2];
47 
48 	carry = 0;
49 	part[1] = 0;
50 	for(i = 0; i < n; i++){
51 		x = part[1] + carry;
52 		if(x < carry)
53 			carry = 1;
54 		else
55 			carry = 0;
56 		y = *p;
57 		mpdigmul(*b++, m, part);
58 		x += part[0];
59 		if(x < part[0])
60 			carry++;
61 		x += y;
62 		if(x < y)
63 			carry++;
64 		*p++ = x;
65 	}
66 	*p = part[1] + carry;
67 }
68 
69 /* prereq: p must have room for n+1 digits */
70 int
71 mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
72 {
73 	int i;
74 	mpdigit x, y, part[2], borrow;
75 
76 	borrow = 0;
77 	part[1] = 0;
78 	for(i = 0; i < n; i++){
79 		x = *p;
80 		y = x - borrow;
81 		if(y > x)
82 			borrow = 1;
83 		else
84 			borrow = 0;
85 		x = part[1];
86 		mpdigmul(*b++, m, part);
87 		x += part[0];
88 		if(x < part[0])
89 			borrow++;
90 		x = y - x;
91 		if(x > y)
92 			borrow++;
93 		*p++ = x;
94 	}
95 
96 	x = *p;
97 	y = x - borrow - part[1];
98 	*p = y;
99 	if(y > x)
100 		return -1;
101 	else
102 		return 1;
103 }