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