| rsc | b3f6179 | 2004-03-21 14:06:38 +0000 | [diff] [blame] | 1 | #include "os.h" |
| 2 | #include <mp.h> |
| 3 | #include "dat.h" |
| 4 | |
| 5 | #define LO(x) ((x) & ((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 | } |
