1 b3f61791 2004-03-21 devnull #include "os.h"
2 b3f61791 2004-03-21 devnull #include <mp.h>
3 b3f61791 2004-03-21 devnull #include "dat.h"
5 9fe7e1a1 2005-12-31 devnull #define LO(x) ((x) & (((mpdigit)1<<(Dbits/2))-1))
6 b3f61791 2004-03-21 devnull #define HI(x) ((x) >> (Dbits/2))
8 b3f61791 2004-03-21 devnull static void
9 b3f61791 2004-03-21 devnull mpdigmul(mpdigit a, mpdigit b, mpdigit *p)
11 b3f61791 2004-03-21 devnull mpdigit x, ah, al, bh, bl, p1, p2, p3, p4;
12 b3f61791 2004-03-21 devnull int carry;
14 cbeb0b26 2006-04-01 devnull /* half digits */
15 b3f61791 2004-03-21 devnull ah = HI(a);
16 b3f61791 2004-03-21 devnull al = LO(a);
17 b3f61791 2004-03-21 devnull bh = HI(b);
18 b3f61791 2004-03-21 devnull bl = LO(b);
20 cbeb0b26 2006-04-01 devnull /* partial products */
21 b3f61791 2004-03-21 devnull p1 = ah*bl;
22 b3f61791 2004-03-21 devnull p2 = bh*al;
23 b3f61791 2004-03-21 devnull p3 = bl*al;
24 b3f61791 2004-03-21 devnull p4 = ah*bh;
26 cbeb0b26 2006-04-01 devnull /* p = ((p1+p2)<<(Dbits/2)) + (p4<<Dbits) + p3 */
27 b3f61791 2004-03-21 devnull carry = 0;
28 b3f61791 2004-03-21 devnull x = p1<<(Dbits/2);
30 b3f61791 2004-03-21 devnull if(p3 < x)
32 b3f61791 2004-03-21 devnull x = p2<<(Dbits/2);
34 b3f61791 2004-03-21 devnull if(p3 < x)
36 cbeb0b26 2006-04-01 devnull p4 += carry + HI(p1) + HI(p2); /* can't carry out of the high digit */
37 b3f61791 2004-03-21 devnull p[0] = p3;
38 b3f61791 2004-03-21 devnull p[1] = p4;
41 cbeb0b26 2006-04-01 devnull /* prereq: p must have room for n+1 digits */
43 b3f61791 2004-03-21 devnull mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
46 b3f61791 2004-03-21 devnull mpdigit carry, x, y, part[2];
48 b3f61791 2004-03-21 devnull carry = 0;
49 b3f61791 2004-03-21 devnull part[1] = 0;
50 b3f61791 2004-03-21 devnull for(i = 0; i < n; i++){
51 b3f61791 2004-03-21 devnull x = part[1] + carry;
52 b3f61791 2004-03-21 devnull if(x < carry)
53 b3f61791 2004-03-21 devnull carry = 1;
55 b3f61791 2004-03-21 devnull carry = 0;
57 b3f61791 2004-03-21 devnull mpdigmul(*b++, m, part);
58 b3f61791 2004-03-21 devnull x += part[0];
59 b3f61791 2004-03-21 devnull if(x < part[0])
62 b3f61791 2004-03-21 devnull if(x < y)
64 b3f61791 2004-03-21 devnull *p++ = x;
66 b3f61791 2004-03-21 devnull *p = part[1] + carry;
69 cbeb0b26 2006-04-01 devnull /* prereq: p must have room for n+1 digits */
71 b3f61791 2004-03-21 devnull mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
74 b3f61791 2004-03-21 devnull mpdigit x, y, part[2], borrow;
76 b3f61791 2004-03-21 devnull borrow = 0;
77 b3f61791 2004-03-21 devnull part[1] = 0;
78 b3f61791 2004-03-21 devnull for(i = 0; i < n; i++){
80 b3f61791 2004-03-21 devnull y = x - borrow;
81 b3f61791 2004-03-21 devnull if(y > x)
82 b3f61791 2004-03-21 devnull borrow = 1;
84 b3f61791 2004-03-21 devnull borrow = 0;
85 b3f61791 2004-03-21 devnull x = part[1];
86 b3f61791 2004-03-21 devnull mpdigmul(*b++, m, part);
87 b3f61791 2004-03-21 devnull x += part[0];
88 b3f61791 2004-03-21 devnull if(x < part[0])
89 b3f61791 2004-03-21 devnull borrow++;
90 b3f61791 2004-03-21 devnull x = y - x;
91 b3f61791 2004-03-21 devnull if(x > y)
92 b3f61791 2004-03-21 devnull borrow++;
93 b3f61791 2004-03-21 devnull *p++ = x;
97 b3f61791 2004-03-21 devnull y = x - borrow - part[1];
99 b3f61791 2004-03-21 devnull if(y > x)
100 b3f61791 2004-03-21 devnull return -1;
102 b3f61791 2004-03-21 devnull return 1;