mpvecdigmuladd.c (1609B)
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 }