Blob
1 #include "os.h"2 #include <mp.h>3 #include "dat.h"5 #define LO(x) ((x) & ((1<<(Dbits/2))-1))6 #define HI(x) ((x) >> (Dbits/2))8 static void9 mpdigmul(mpdigit a, mpdigit b, mpdigit *p)10 {11 mpdigit x, ah, al, bh, bl, p1, p2, p3, p4;12 int carry;14 // half digits15 ah = HI(a);16 al = LO(a);17 bh = HI(b);18 bl = LO(b);20 // partial products21 p1 = ah*bl;22 p2 = bh*al;23 p3 = bl*al;24 p4 = ah*bh;26 // p = ((p1+p2)<<(Dbits/2)) + (p4<<Dbits) + p327 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 digit37 p[0] = p3;38 p[1] = p4;39 }41 // prereq: p must have room for n+1 digits42 void43 mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)44 {45 int i;46 mpdigit carry, x, y, part[2];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 else55 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 }69 // prereq: p must have room for n+1 digits70 int71 mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)72 {73 int i;74 mpdigit x, y, part[2], borrow;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 else84 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 }96 x = *p;97 y = x - borrow - part[1];98 *p = y;99 if(y > x)100 return -1;101 else102 return 1;103 }