1 /* Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
2 /* All Rights Reserved */
3
4
5 /*
6 * Copyright (c) 1980 Regents of the University of California.
7 * All rights reserved. The Berkeley software License Agreement
8 * specifies the terms and conditions for redistribution.
9 */
10 /* Portions Copyright(c) 1988, Sun Microsystems Inc. */
11 /* All Rights Reserved */
12
13 /*
14 * Copyright (c) 1997, by Sun Microsystems, Inc.
15 * All rights reserved.
16 */
17
18 /*
19 * Copyright (c) 2018, Joyent, Inc.
20 */
21
22 /* LINTLIBRARY */
23
24 #include <mp.h>
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <sys/types.h>
28 #include "libmp.h"
29
30 static void m_div(MINT *, MINT *, MINT *, MINT *);
31
32 void
33 mp_mdiv(MINT *a, MINT *b, MINT *q, MINT *r)
34 {
35 MINT x, y;
36 int sign;
37
38 sign = 1;
39 x.len = y.len = 0;
40 _mp_move(a, &x);
41 _mp_move(b, &y);
42 if (x.len < 0) {
43 sign = -1;
44 x.len = -x.len;
45 }
46 if (y.len < 0) {
47 sign = -sign;
48 y.len = -y.len;
49 }
50 _mp_xfree(q);
51 _mp_xfree(r);
52 m_div(&x, &y, q, r);
53 if (sign == -1) {
54 q->len = -q->len;
55 r->len = -r->len;
56 }
57 _mp_xfree(&x);
58 _mp_xfree(&y);
59 }
60
61 static int
62 m_dsb(int qx, int n, short *a, short *b)
63 {
64 int borrow;
65 int s3b2shit;
66 int j;
67 short fifteen = 15;
68 short *aptr, *bptr;
69 #ifdef DEBUGDSB
70 (void) printf("m_dsb %d %d %d %d\n", qx, n, *a, *b);
71 #endif
72
73 borrow = 0;
74 aptr = a;
75 bptr = b;
76 for (j = n; j > 0; j--) {
77 #ifdef DEBUGDSB
78 (void) printf("1 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
79 *bptr, *aptr);
80 #endif
81 borrow -= (*aptr++) * qx - *bptr;
82 #ifdef DEBUGDSB
83 (void) printf("2 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
84 *bptr, *aptr);
85 #endif
86 *bptr++ = (short)(borrow & 077777);
87 #ifdef DEBUGDSB
88 (void) printf("3 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
89 *bptr, *aptr);
90 #endif
91 if (borrow >= 0) borrow >>= fifteen; /* 3b2 */
92 else borrow = 0xfffe0000 | (borrow >> fifteen);
93 #ifdef DEBUGDSB
94 (void) printf("4 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
95 *bptr, *aptr);
96 #endif
97 }
98 borrow += *bptr;
99 *bptr = (short)(borrow & 077777);
100 if (borrow >= 0) s3b2shit = borrow >> fifteen; /* 3b2 */
101 else s3b2shit = 0xfffe0000 | (borrow >> fifteen);
102 if (s3b2shit == 0) {
103 #ifdef DEBUGDSB
104 (void) printf("mdsb 0\n");
105 #endif
106 return (0);
107 }
108 borrow = 0;
109 aptr = a;
110 bptr = b;
111 for (j = n; j > 0; j--) {
112 borrow += *aptr++ + *bptr;
113 *bptr++ = (short)(borrow & 077777);
114 if (borrow >= 0) borrow >>= fifteen; /* 3b2 */
115 else borrow = 0xfffe0000 | (borrow >>fifteen);
116 }
117 #ifdef DEBUGDSB
118 (void) printf("mdsb 1\n");
119 #endif
120 return (1);
121 }
122
123 static int
124 m_trq(short v1, short v2, short u1, short u2, short u3)
125 {
126 short d;
127 int x1;
128 int c1;
129
130 c1 = u1 * 0100000 + u2;
131 if (u1 == v1) {
132 d = 077777;
133 } else {
134 d = (short)(c1 / v1);
135 }
136 do {
137 x1 = c1 - v1 * d;
138 x1 = x1 * 0100000 + u3 - v2 * d;
139 --d;
140 } while (x1 < 0);
141 #ifdef DEBUGMTRQ
142 (void) printf("mtrq %d %d %d %d %d %d\n", v1, v2, u1, u2, u3, (d+1));
143 #endif
144 return ((int)d + 1);
145 }
146
147 static void
148 m_div(MINT *a, MINT *b, MINT *q, MINT *r)
149 {
150 MINT u, v, x, w;
151 short d;
152 short *qval;
153 short *uval;
154 int j;
155 int qq;
156 int n;
157 short v1;
158 short v2;
159
160 u.len = v.len = x.len = w.len = 0;
161 if (b->len == 0) {
162 _mp_fatal("mdiv divide by zero");
163 return;
164 }
165 if (b->len == 1) {
166 r->val = _mp_xalloc(1, "m_div1");
167 mp_sdiv(a, b->val[0], q, r->val);
168 if (r->val[0] == 0) {
169 free(r->val);
170 r->len = 0;
171 } else {
172 r->len = 1;
173 }
174 return;
175 }
176 if (a -> len < b -> len) {
177 q->len = 0;
178 r->len = a->len;
179 r->val = _mp_xalloc(r->len, "m_div2");
180 for (qq = 0; qq < r->len; qq++) {
181 r->val[qq] = a->val[qq];
182 }
183 return;
184 }
185 x.len = 1;
186 x.val = &d;
187 n = b->len;
188 d = 0100000 / (b->val[n - 1] + 1);
189 mp_mult(a, &x, &u); /* subtle: relies on mult allocing extra space */
190 mp_mult(b, &x, &v);
191 #ifdef DEBUG_MDIV
192 (void) printf(" u=%s\n", mtox(&u));
193 (void) printf(" v=%s\n", mtox(&v));
194 #endif
195 v1 = v.val[n - 1];
196 v2 = v.val[n - 2];
197 qval = _mp_xalloc(a -> len - n + 1, "m_div3");
198 uval = u.val;
199 for (j = a->len - n; j >= 0; j--) {
200 qq = m_trq(v1, v2, uval[j + n], uval[j + n - 1],
201 uval[j + n - 2]);
202 if (m_dsb(qq, n, v.val, uval + j))
203 qq -= 1;
204 qval[j] = (short)qq;
205 }
206 x.len = n;
207 x.val = u.val;
208 _mp_mcan(&x);
209 #ifdef DEBUG_MDIV
210 (void) printf(" x=%s\n", mtox(&x));
211 (void) printf(" d(in)=%d\n", (d));
212 #endif
213 mp_sdiv(&x, d, &w, &d);
214 #ifdef DEBUG_MDIV
215 (void) printf(" w=%s\n", mtox(&w));
216 (void) printf(" d(out)=%d\n", (d));
217 #endif
218 r->len = w.len;
219 r->val = w.val;
220 q->val = qval;
221 qq = a->len - n + 1;
222 if (qq > 0 && qval[qq - 1] == 0)
223 qq -= 1;
224 q->len = qq;
225 if (qq == 0)
226 free(qval);
227 if (x.len != 0)
228 _mp_xfree(&u);
229 _mp_xfree(&v);
230 }