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