1*8ccd4a63SDavid du Colombier #include "os.h"
2*8ccd4a63SDavid du Colombier #include <mp.h>
3*8ccd4a63SDavid du Colombier #include "dat.h"
4*8ccd4a63SDavid du Colombier
5*8ccd4a63SDavid du Colombier #define LO(x) ((x) & ((1<<(Dbits/2))-1))
6*8ccd4a63SDavid du Colombier #define HI(x) ((x) >> (Dbits/2))
7*8ccd4a63SDavid du Colombier
8*8ccd4a63SDavid du Colombier static void
mpdigmul(mpdigit a,mpdigit b,mpdigit * p)9*8ccd4a63SDavid du Colombier mpdigmul(mpdigit a, mpdigit b, mpdigit *p)
10*8ccd4a63SDavid du Colombier {
11*8ccd4a63SDavid du Colombier mpdigit x, ah, al, bh, bl, p1, p2, p3, p4;
12*8ccd4a63SDavid du Colombier int carry;
13*8ccd4a63SDavid du Colombier
14*8ccd4a63SDavid du Colombier // half digits
15*8ccd4a63SDavid du Colombier ah = HI(a);
16*8ccd4a63SDavid du Colombier al = LO(a);
17*8ccd4a63SDavid du Colombier bh = HI(b);
18*8ccd4a63SDavid du Colombier bl = LO(b);
19*8ccd4a63SDavid du Colombier
20*8ccd4a63SDavid du Colombier // partial products
21*8ccd4a63SDavid du Colombier p1 = ah*bl;
22*8ccd4a63SDavid du Colombier p2 = bh*al;
23*8ccd4a63SDavid du Colombier p3 = bl*al;
24*8ccd4a63SDavid du Colombier p4 = ah*bh;
25*8ccd4a63SDavid du Colombier
26*8ccd4a63SDavid du Colombier // p = ((p1+p2)<<(Dbits/2)) + (p4<<Dbits) + p3
27*8ccd4a63SDavid du Colombier carry = 0;
28*8ccd4a63SDavid du Colombier x = p1<<(Dbits/2);
29*8ccd4a63SDavid du Colombier p3 += x;
30*8ccd4a63SDavid du Colombier if(p3 < x)
31*8ccd4a63SDavid du Colombier carry++;
32*8ccd4a63SDavid du Colombier x = p2<<(Dbits/2);
33*8ccd4a63SDavid du Colombier p3 += x;
34*8ccd4a63SDavid du Colombier if(p3 < x)
35*8ccd4a63SDavid du Colombier carry++;
36*8ccd4a63SDavid du Colombier p4 += carry + HI(p1) + HI(p2); // can't carry out of the high digit
37*8ccd4a63SDavid du Colombier p[0] = p3;
38*8ccd4a63SDavid du Colombier p[1] = p4;
39*8ccd4a63SDavid du Colombier }
40*8ccd4a63SDavid du Colombier
41*8ccd4a63SDavid du Colombier // prereq: p must have room for n+1 digits
42*8ccd4a63SDavid du Colombier void
mpvecdigmuladd(mpdigit * b,int n,mpdigit m,mpdigit * p)43*8ccd4a63SDavid du Colombier mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
44*8ccd4a63SDavid du Colombier {
45*8ccd4a63SDavid du Colombier int i;
46*8ccd4a63SDavid du Colombier mpdigit carry, x, y, part[2];
47*8ccd4a63SDavid du Colombier
48*8ccd4a63SDavid du Colombier carry = 0;
49*8ccd4a63SDavid du Colombier part[1] = 0;
50*8ccd4a63SDavid du Colombier for(i = 0; i < n; i++){
51*8ccd4a63SDavid du Colombier x = part[1] + carry;
52*8ccd4a63SDavid du Colombier if(x < carry)
53*8ccd4a63SDavid du Colombier carry = 1;
54*8ccd4a63SDavid du Colombier else
55*8ccd4a63SDavid du Colombier carry = 0;
56*8ccd4a63SDavid du Colombier y = *p;
57*8ccd4a63SDavid du Colombier mpdigmul(*b++, m, part);
58*8ccd4a63SDavid du Colombier x += part[0];
59*8ccd4a63SDavid du Colombier if(x < part[0])
60*8ccd4a63SDavid du Colombier carry++;
61*8ccd4a63SDavid du Colombier x += y;
62*8ccd4a63SDavid du Colombier if(x < y)
63*8ccd4a63SDavid du Colombier carry++;
64*8ccd4a63SDavid du Colombier *p++ = x;
65*8ccd4a63SDavid du Colombier }
66*8ccd4a63SDavid du Colombier *p = part[1] + carry;
67*8ccd4a63SDavid du Colombier }
68*8ccd4a63SDavid du Colombier
69*8ccd4a63SDavid du Colombier // prereq: p must have room for n+1 digits
70*8ccd4a63SDavid du Colombier int
mpvecdigmulsub(mpdigit * b,int n,mpdigit m,mpdigit * p)71*8ccd4a63SDavid du Colombier mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
72*8ccd4a63SDavid du Colombier {
73*8ccd4a63SDavid du Colombier int i;
74*8ccd4a63SDavid du Colombier mpdigit x, y, part[2], borrow;
75*8ccd4a63SDavid du Colombier
76*8ccd4a63SDavid du Colombier borrow = 0;
77*8ccd4a63SDavid du Colombier part[1] = 0;
78*8ccd4a63SDavid du Colombier for(i = 0; i < n; i++){
79*8ccd4a63SDavid du Colombier x = *p;
80*8ccd4a63SDavid du Colombier y = x - borrow;
81*8ccd4a63SDavid du Colombier if(y > x)
82*8ccd4a63SDavid du Colombier borrow = 1;
83*8ccd4a63SDavid du Colombier else
84*8ccd4a63SDavid du Colombier borrow = 0;
85*8ccd4a63SDavid du Colombier x = part[1];
86*8ccd4a63SDavid du Colombier mpdigmul(*b++, m, part);
87*8ccd4a63SDavid du Colombier x += part[0];
88*8ccd4a63SDavid du Colombier if(x < part[0])
89*8ccd4a63SDavid du Colombier borrow++;
90*8ccd4a63SDavid du Colombier x = y - x;
91*8ccd4a63SDavid du Colombier if(x > y)
92*8ccd4a63SDavid du Colombier borrow++;
93*8ccd4a63SDavid du Colombier *p++ = x;
94*8ccd4a63SDavid du Colombier }
95*8ccd4a63SDavid du Colombier
96*8ccd4a63SDavid du Colombier x = *p;
97*8ccd4a63SDavid du Colombier y = x - borrow - part[1];
98*8ccd4a63SDavid du Colombier *p = y;
99*8ccd4a63SDavid du Colombier if(y > x)
100*8ccd4a63SDavid du Colombier return -1;
101*8ccd4a63SDavid du Colombier else
102*8ccd4a63SDavid du Colombier return 1;
103*8ccd4a63SDavid du Colombier }
104