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
77 static struct fpn * __fpu_modrem(struct fpemu *fe, int is_mod);
78 static int abscmp3(const struct fpn *a, const struct fpn *b);
79
80 /* Absolute FORTRAN Compare */
81 static int
abscmp3(const struct fpn * a,const struct fpn * b)82 abscmp3(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
101 static struct fpn *
__fpu_modrem(struct fpemu * fe,int is_mod)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
217 struct fpn *
fpu_rem(struct fpemu * fe)218 fpu_rem(struct fpemu *fe)
219 {
220 return __fpu_modrem(fe, 0);
221 }
222
223 struct fpn *
fpu_mod(struct fpemu * fe)224 fpu_mod(struct fpemu *fe)
225 {
226 return __fpu_modrem(fe, 1);
227 }
228