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"
  4 b3f61791 2004-03-21 devnull 
  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))
  7 b3f61791 2004-03-21 devnull 
  8 b3f61791 2004-03-21 devnull static void
  9 b3f61791 2004-03-21 devnull mpdigmul(mpdigit a, mpdigit b, mpdigit *p)
 10 b3f61791 2004-03-21 devnull {
 11 b3f61791 2004-03-21 devnull 	mpdigit x, ah, al, bh, bl, p1, p2, p3, p4;
 12 b3f61791 2004-03-21 devnull 	int carry;
 13 b3f61791 2004-03-21 devnull 
 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);
 19 b3f61791 2004-03-21 devnull 
 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;
 25 b3f61791 2004-03-21 devnull 
 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);
 29 b3f61791 2004-03-21 devnull 	p3 += x;
 30 b3f61791 2004-03-21 devnull 	if(p3 < x)
 31 b3f61791 2004-03-21 devnull 		carry++;
 32 b3f61791 2004-03-21 devnull 	x = p2<<(Dbits/2);
 33 b3f61791 2004-03-21 devnull 	p3 += x;
 34 b3f61791 2004-03-21 devnull 	if(p3 < x)
 35 b3f61791 2004-03-21 devnull 		carry++;
 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;
 39 b3f61791 2004-03-21 devnull }
 40 b3f61791 2004-03-21 devnull 
 41 cbeb0b26 2006-04-01 devnull /* prereq: p must have room for n+1 digits */
 42 b3f61791 2004-03-21 devnull void
 43 b3f61791 2004-03-21 devnull mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
 44 b3f61791 2004-03-21 devnull {
 45 b3f61791 2004-03-21 devnull 	int i;
 46 b3f61791 2004-03-21 devnull 	mpdigit carry, x, y, part[2];
 47 b3f61791 2004-03-21 devnull 
 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;
 54 b3f61791 2004-03-21 devnull 		else
 55 b3f61791 2004-03-21 devnull 			carry = 0;
 56 b3f61791 2004-03-21 devnull 		y = *p;
 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])
 60 b3f61791 2004-03-21 devnull 			carry++;
 61 b3f61791 2004-03-21 devnull 		x += y;
 62 b3f61791 2004-03-21 devnull 		if(x < y)
 63 b3f61791 2004-03-21 devnull 			carry++;
 64 b3f61791 2004-03-21 devnull 		*p++ = x;
 65 b3f61791 2004-03-21 devnull 	}
 66 b3f61791 2004-03-21 devnull 	*p = part[1] + carry;
 67 b3f61791 2004-03-21 devnull }
 68 b3f61791 2004-03-21 devnull 
 69 cbeb0b26 2006-04-01 devnull /* prereq: p must have room for n+1 digits */
 70 b3f61791 2004-03-21 devnull int
 71 b3f61791 2004-03-21 devnull mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
 72 b3f61791 2004-03-21 devnull {
 73 b3f61791 2004-03-21 devnull 	int i;
 74 b3f61791 2004-03-21 devnull 	mpdigit x, y, part[2], borrow;
 75 b3f61791 2004-03-21 devnull 
 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++){
 79 b3f61791 2004-03-21 devnull 		x = *p;
 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;
 83 b3f61791 2004-03-21 devnull 		else
 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;
 94 b3f61791 2004-03-21 devnull 	}
 95 b3f61791 2004-03-21 devnull 
 96 b3f61791 2004-03-21 devnull 	x = *p;
 97 b3f61791 2004-03-21 devnull 	y = x - borrow - part[1];
 98 b3f61791 2004-03-21 devnull 	*p = y;
 99 b3f61791 2004-03-21 devnull 	if(y > x)
100 b3f61791 2004-03-21 devnull 		return -1;
101 b3f61791 2004-03-21 devnull 	else
102 b3f61791 2004-03-21 devnull 		return 1;
103 b3f61791 2004-03-21 devnull }