1/*	$NetBSD: fpu_rem.c,v 1.18 2023/11/19 03:58:15 isaki Exp $	*/
2
3/*
4 * Copyright (c) 1995  Ken Nakata
5 *	All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code must retain the above copyright
11 *    notice, this list of conditions and the following disclaimer.
12 * 2. Redistributions in binary form must reproduce the above copyright
13 *    notice, this list of conditions and the following disclaimer in the
14 *    documentation and/or other materials provided with the distribution.
15 * 3. Neither the name of the author nor the names of its contributors
16 *    may be used to endorse or promote products derived from this software
17 *    without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
29 * SUCH DAMAGE.
30 *
31 *	@(#)fpu_rem.c	10/24/95
32 */
33
34#include <sys/cdefs.h>
35__KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v 1.18 2023/11/19 03:58:15 isaki Exp $");
36
37#include <sys/types.h>
38#include <sys/signal.h>
39#include <machine/frame.h>
40
41#include "fpu_emulate.h"
42
43/*
44 *       ALGORITHM
45 *
46 *       Step 1.  Save and strip signs of X and Y: signX := sign(X),
47 *                signY := sign(Y), X := *X*, Y := *Y*,
48 *                signQ := signX EOR signY. Record whether MOD or REM
49 *                is requested.
50 *
51 *       Step 2.  Set L := expo(X)-expo(Y), Q := 0.
52 *                If (L < 0) then
53 *                   R := X, go to Step 4.
54 *                else
55 *                   R := 2^(-L)X, j := L.
56 *                endif
57 *
58 *       Step 3.  Perform MOD(X,Y)
59 *            3.1 If R = Y, then { Q := Q + 1, R := 0, go to Step 7. }
60 *            3.2 If R > Y, then { R := R - Y, Q := Q + 1}
61 *            3.3 If j = 0, go to Step 4.
62 *            3.4 j := j - 1, Q := 2Q, R := 2R. Go to Step 3.1.
63 *
64 *       Step 4.  R := signX*R.
65 *
66 *       Step 5.  If MOD is requested, go to Step 7.
67 *
68 *       Step 6.  Now, R = MOD(X,Y), convert to REM(X,Y) is requested.
69 *                Do banker's rounding.
70 *                If   abs(R) > Y/2
71 *                 || (abs(R) == Y/2 && Q % 2 == 1) then
72 *                 { Q := Q + 1, R := R - signX * Y }.
73 *
74 *       Step 7.  Return signQ, last 7 bits of Q, and R as required.
75 */
76
77static struct fpn * __fpu_modrem(struct fpemu *fe, int is_mod);
78static int abscmp3(const struct fpn *a, const struct fpn *b);
79
80/* Absolute FORTRAN Compare */
81static int
82abscmp3(const struct fpn *a, const struct fpn *b)
83{
84	int i;
85
86	if (a->fp_exp < b->fp_exp) {
87		return -1;
88	} else if (a->fp_exp > b->fp_exp) {
89		return 1;
90	} else {
91		for (i = 0; i < 3; i++) {
92			if (a->fp_mant[i] < b->fp_mant[i])
93				return -1;
94			else if (a->fp_mant[i] > b->fp_mant[i])
95				return 1;
96		}
97	}
98	return 0;
99}
100
101static struct fpn *
102__fpu_modrem(struct fpemu *fe, int is_mod)
103{
104	static struct fpn X, Y;
105	struct fpn *x, *y, *r;
106	uint32_t signX, signY, signQ;
107	int j, l, q;
108	int cmp;
109
110	if (ISNAN(&fe->fe_f1) || ISNAN(&fe->fe_f2))
111		return fpu_newnan(fe);
112	if (ISINF(&fe->fe_f1) || ISZERO(&fe->fe_f2))
113		return fpu_newnan(fe);
114
115	CPYFPN(&X, &fe->fe_f1);
116	CPYFPN(&Y, &fe->fe_f2);
117	x = &X;
118	y = &Y;
119	q = 0;
120	r = &fe->fe_f2;
121
122	/*
123	 * Step 1
124	 */
125	signX = x->fp_sign;
126	signY = y->fp_sign;
127	signQ = (signX ^ signY);
128	x->fp_sign = y->fp_sign = 0;
129
130	/* Special treatment that just return input value but Q is necessary */
131	if (ISZERO(x) || ISINF(y)) {
132		r = &fe->fe_f1;
133		goto Step7;
134	}
135
136	/*
137	 * Step 2
138	 */
139	l = x->fp_exp - y->fp_exp;
140	CPYFPN(r, x);
141	if (l >= 0) {
142		r->fp_exp -= l;
143		j = l;
144
145		/*
146		 * Step 3
147		 */
148		for (;;) {
149			cmp = abscmp3(r, y);
150
151			/* Step 3.1 */
152			if (cmp == 0)
153				break;
154
155			/* Step 3.2 */
156			if (cmp > 0) {
157				CPYFPN(&fe->fe_f1, r);
158				CPYFPN(&fe->fe_f2, y);
159				fe->fe_f2.fp_sign = 1;
160				r = fpu_add(fe);
161				q++;
162			}
163
164			/* Step 3.3 */
165			if (j == 0)
166				goto Step4;
167
168			/* Step 3.4 */
169			j--;
170			q += q;
171			r->fp_exp++;
172		}
173		/* R == Y */
174		q++;
175		r->fp_class = FPC_ZERO;
176		goto Step7;
177	}
178 Step4:
179	r->fp_sign = signX;
180
181	/*
182	 * Step 5
183	 */
184	if (is_mod)
185		goto Step7;
186
187	/*
188	 * Step 6
189	 */
190	/* y = y / 2 */
191	y->fp_exp--;
192	/* abscmp3 ignore sign */
193	cmp = abscmp3(r, y);
194	/* revert y */
195	y->fp_exp++;
196
197	if (cmp > 0 || (cmp == 0 && q % 2)) {
198		q++;
199		CPYFPN(&fe->fe_f1, r);
200		CPYFPN(&fe->fe_f2, y);
201		fe->fe_f2.fp_sign = !signX;
202		r = fpu_add(fe);
203	}
204
205	/*
206	 * Step 7
207	 */
208 Step7:
209	q &= 0x7f;
210	q |= (signQ << 7);
211	fe->fe_fpframe->fpf_fpsr =
212	fe->fe_fpsr =
213	    (fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
214	return r;
215}
216
217struct fpn *
218fpu_rem(struct fpemu *fe)
219{
220	return __fpu_modrem(fe, 0);
221}
222
223struct fpn *
224fpu_mod(struct fpemu *fe)
225{
226	return __fpu_modrem(fe, 1);
227}
228