xref: /illumos-gate/usr/src/lib/libmp/common/mdiv.c (revision 981fe1b1a6ca86e9c41a706462f976f947e6f465)
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
mp_mdiv(MINT * a,MINT * b,MINT * q,MINT * r)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
m_dsb(int qx,int n,short * a,short * b)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
m_trq(short v1,short v2,short u1,short u2,short u3)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
m_div(MINT * a,MINT * b,MINT * q,MINT * r)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 }
231