xref: /titanic_41/usr/src/lib/libmp/common/mdiv.c (revision fa9e4066f08beec538e775443c5be79dd423fcab)
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 }
229