harvey/sys/src/libmp/port/mpvecdigmuladd.c

#include "os.h"
#include <mp.h>
#include "dat.h"

#define LO(x) ((x) & ((1<<(Dbits/2))-1))
#define HI(x) ((x) >> (Dbits/2))

static void
mpdigmul(mpdigit a, mpdigit b, mpdigit *p)
{
	mpdigit x, ah, al, bh, bl, p1, p2, p3, p4;
	int carry;

	// half digits
	ah = HI(a);
	al = LO(a);
	bh = HI(b);
	bl = LO(b);

	// partial products
	p1 = ah*bl;
	p2 = bh*al;
	p3 = bl*al;
	p4 = ah*bh;

	// p = ((p1+p2)<<(Dbits/2)) + (p4<<Dbits) + p3
	carry = 0;
	x = p1<<(Dbits/2);
	p3 += x;
	if(p3 < x)
		carry++;
	x = p2<<(Dbits/2);
	p3 += x;
	if(p3 < x)
		carry++;
	p4 += carry + HI(p1) + HI(p2);	// can't carry out of the high digit
	p[0] = p3;
	p[1] = p4;
}

// prereq: p must have room for n+1 digits
void
mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
{
	int i;
	mpdigit carry, x, y, part[2];

	carry = 0;
	part[1] = 0;
	for(i = 0; i < n; i++){
		x = part[1] + carry;
		if(x < carry)
			carry = 1;
		else
			carry = 0;
		y = *p;
		mpdigmul(*b++, m, part);
		x += part[0];
		if(x < part[0])
			carry++;
		x += y;
		if(x < y)
			carry++;
		*p++ = x;
	}
	*p = part[1] + carry;
}

// prereq: p must have room for n+1 digits
int
mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
{
	int i;
	mpdigit x, y, part[2], borrow;

	borrow = 0;
	part[1] = 0;
	for(i = 0; i < n; i++){
		x = *p;
		y = x - borrow;
		if(y > x)
			borrow = 1;
		else
			borrow = 0;
		x = part[1];
		mpdigmul(*b++, m, part);
		x += part[0];
		if(x < part[0])
			borrow++;
		x = y - x;
		if(x > y)
			borrow++;
		*p++ = x;
	}

	x = *p;
	y = x - borrow - part[1];
	*p = y;
	if(y > x)
		return -1;
	else
		return 1;
}