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