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