104 lines
1.5 KiB
C
104 lines
1.5 KiB
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;
|
||
|
|
}
|