1 // SPDX-License-Identifier: GPL-2.0-or-later 1 // SPDX-License-Identifier: GPL-2.0-or-later
2 /* 2 /*
3 3
4 fp_arith.c: floating-point math routines fo 4 fp_arith.c: floating-point math routines for the Linux-m68k
5 floating point emulator. 5 floating point emulator.
6 6
7 Copyright (c) 1998-1999 David Huggins-Daine 7 Copyright (c) 1998-1999 David Huggins-Daines.
8 8
9 Somewhat based on the AlphaLinux floating p 9 Somewhat based on the AlphaLinux floating point emulator, by David
10 Mosberger-Tang. 10 Mosberger-Tang.
11 11
12 */ 12 */
13 13
14 #include "fp_emu.h" 14 #include "fp_emu.h"
15 #include "multi_arith.h" 15 #include "multi_arith.h"
16 #include "fp_arith.h" 16 #include "fp_arith.h"
17 17
18 const struct fp_ext fp_QNaN = 18 const struct fp_ext fp_QNaN =
19 { 19 {
20 .exp = 0x7fff, 20 .exp = 0x7fff,
21 .mant = { .m64 = ~0 } 21 .mant = { .m64 = ~0 }
22 }; 22 };
23 23
24 const struct fp_ext fp_Inf = 24 const struct fp_ext fp_Inf =
25 { 25 {
26 .exp = 0x7fff, 26 .exp = 0x7fff,
27 }; 27 };
28 28
29 /* let's start with the easy ones */ 29 /* let's start with the easy ones */
30 30
31 struct fp_ext *fp_fabs(struct fp_ext *dest, st 31 struct fp_ext *fp_fabs(struct fp_ext *dest, struct fp_ext *src)
32 { 32 {
33 dprint(PINSTR, "fabs\n"); 33 dprint(PINSTR, "fabs\n");
34 34
35 fp_monadic_check(dest, src); 35 fp_monadic_check(dest, src);
36 36
37 dest->sign = 0; 37 dest->sign = 0;
38 38
39 return dest; 39 return dest;
40 } 40 }
41 41
42 struct fp_ext *fp_fneg(struct fp_ext *dest, st 42 struct fp_ext *fp_fneg(struct fp_ext *dest, struct fp_ext *src)
43 { 43 {
44 dprint(PINSTR, "fneg\n"); 44 dprint(PINSTR, "fneg\n");
45 45
46 fp_monadic_check(dest, src); 46 fp_monadic_check(dest, src);
47 47
48 dest->sign = !dest->sign; 48 dest->sign = !dest->sign;
49 49
50 return dest; 50 return dest;
51 } 51 }
52 52
53 /* Now, the slightly harder ones */ 53 /* Now, the slightly harder ones */
54 54
55 /* fp_fadd: Implements the kernel of the FADD, 55 /* fp_fadd: Implements the kernel of the FADD, FSADD, FDADD, FSUB,
56 FDSUB, and FCMP instructions. */ 56 FDSUB, and FCMP instructions. */
57 57
58 struct fp_ext *fp_fadd(struct fp_ext *dest, st 58 struct fp_ext *fp_fadd(struct fp_ext *dest, struct fp_ext *src)
59 { 59 {
60 int diff; 60 int diff;
61 61
62 dprint(PINSTR, "fadd\n"); 62 dprint(PINSTR, "fadd\n");
63 63
64 fp_dyadic_check(dest, src); 64 fp_dyadic_check(dest, src);
65 65
66 if (IS_INF(dest)) { 66 if (IS_INF(dest)) {
67 /* infinity - infinity == NaN 67 /* infinity - infinity == NaN */
68 if (IS_INF(src) && (src->sign 68 if (IS_INF(src) && (src->sign != dest->sign))
69 fp_set_nan(dest); 69 fp_set_nan(dest);
70 return dest; 70 return dest;
71 } 71 }
72 if (IS_INF(src)) { 72 if (IS_INF(src)) {
73 fp_copy_ext(dest, src); 73 fp_copy_ext(dest, src);
74 return dest; 74 return dest;
75 } 75 }
76 76
77 if (IS_ZERO(dest)) { 77 if (IS_ZERO(dest)) {
78 if (IS_ZERO(src)) { 78 if (IS_ZERO(src)) {
79 if (src->sign != dest- 79 if (src->sign != dest->sign) {
80 if (FPDATA->rn 80 if (FPDATA->rnd == FPCR_ROUND_RM)
81 dest-> 81 dest->sign = 1;
82 else 82 else
83 dest-> 83 dest->sign = 0;
84 } 84 }
85 } else 85 } else
86 fp_copy_ext(dest, src) 86 fp_copy_ext(dest, src);
87 return dest; 87 return dest;
88 } 88 }
89 89
90 dest->lowmant = src->lowmant = 0; 90 dest->lowmant = src->lowmant = 0;
91 91
92 if ((diff = dest->exp - src->exp) > 0) 92 if ((diff = dest->exp - src->exp) > 0)
93 fp_denormalize(src, diff); 93 fp_denormalize(src, diff);
94 else if ((diff = -diff) > 0) 94 else if ((diff = -diff) > 0)
95 fp_denormalize(dest, diff); 95 fp_denormalize(dest, diff);
96 96
97 if (dest->sign == src->sign) { 97 if (dest->sign == src->sign) {
98 if (fp_addmant(dest, src)) 98 if (fp_addmant(dest, src))
99 if (!fp_addcarry(dest) 99 if (!fp_addcarry(dest))
100 return dest; 100 return dest;
101 } else { 101 } else {
102 if (dest->mant.m64 < src->mant 102 if (dest->mant.m64 < src->mant.m64) {
103 fp_submant(dest, src, 103 fp_submant(dest, src, dest);
104 dest->sign = !dest->si 104 dest->sign = !dest->sign;
105 } else 105 } else
106 fp_submant(dest, dest, 106 fp_submant(dest, dest, src);
107 } 107 }
108 108
109 return dest; 109 return dest;
110 } 110 }
111 111
112 /* fp_fsub: Implements the kernel of the FSUB, 112 /* fp_fsub: Implements the kernel of the FSUB, FSSUB, and FDSUB
113 instructions. 113 instructions.
114 114
115 Remember that the arguments are in assemble 115 Remember that the arguments are in assembler-syntax order! */
116 116
117 struct fp_ext *fp_fsub(struct fp_ext *dest, st 117 struct fp_ext *fp_fsub(struct fp_ext *dest, struct fp_ext *src)
118 { 118 {
119 dprint(PINSTR, "fsub "); 119 dprint(PINSTR, "fsub ");
120 120
121 src->sign = !src->sign; 121 src->sign = !src->sign;
122 return fp_fadd(dest, src); 122 return fp_fadd(dest, src);
123 } 123 }
124 124
125 125
126 struct fp_ext *fp_fcmp(struct fp_ext *dest, st 126 struct fp_ext *fp_fcmp(struct fp_ext *dest, struct fp_ext *src)
127 { 127 {
128 dprint(PINSTR, "fcmp "); 128 dprint(PINSTR, "fcmp ");
129 129
130 FPDATA->temp[1] = *dest; 130 FPDATA->temp[1] = *dest;
131 src->sign = !src->sign; 131 src->sign = !src->sign;
132 return fp_fadd(&FPDATA->temp[1], src); 132 return fp_fadd(&FPDATA->temp[1], src);
133 } 133 }
134 134
135 struct fp_ext *fp_ftst(struct fp_ext *dest, st 135 struct fp_ext *fp_ftst(struct fp_ext *dest, struct fp_ext *src)
136 { 136 {
137 dprint(PINSTR, "ftst\n"); 137 dprint(PINSTR, "ftst\n");
138 138
139 (void)dest; 139 (void)dest;
140 140
141 return src; 141 return src;
142 } 142 }
143 143
144 struct fp_ext *fp_fmul(struct fp_ext *dest, st 144 struct fp_ext *fp_fmul(struct fp_ext *dest, struct fp_ext *src)
145 { 145 {
146 union fp_mant128 temp; 146 union fp_mant128 temp;
147 int exp; 147 int exp;
148 148
149 dprint(PINSTR, "fmul\n"); 149 dprint(PINSTR, "fmul\n");
150 150
151 fp_dyadic_check(dest, src); 151 fp_dyadic_check(dest, src);
152 152
153 /* calculate the correct sign now, as 153 /* calculate the correct sign now, as it's necessary for infinities */
154 dest->sign = src->sign ^ dest->sign; 154 dest->sign = src->sign ^ dest->sign;
155 155
156 /* Handle infinities */ 156 /* Handle infinities */
157 if (IS_INF(dest)) { 157 if (IS_INF(dest)) {
158 if (IS_ZERO(src)) 158 if (IS_ZERO(src))
159 fp_set_nan(dest); 159 fp_set_nan(dest);
160 return dest; 160 return dest;
161 } 161 }
162 if (IS_INF(src)) { 162 if (IS_INF(src)) {
163 if (IS_ZERO(dest)) 163 if (IS_ZERO(dest))
164 fp_set_nan(dest); 164 fp_set_nan(dest);
165 else 165 else
166 fp_copy_ext(dest, src) 166 fp_copy_ext(dest, src);
167 return dest; 167 return dest;
168 } 168 }
169 169
170 /* Of course, as we all know, zero * a 170 /* Of course, as we all know, zero * anything = zero. You may
171 not have known that it might be a p 171 not have known that it might be a positive or negative
172 zero... */ 172 zero... */
173 if (IS_ZERO(dest) || IS_ZERO(src)) { 173 if (IS_ZERO(dest) || IS_ZERO(src)) {
174 dest->exp = 0; 174 dest->exp = 0;
175 dest->mant.m64 = 0; 175 dest->mant.m64 = 0;
176 dest->lowmant = 0; 176 dest->lowmant = 0;
177 177
178 return dest; 178 return dest;
179 } 179 }
180 180
181 exp = dest->exp + src->exp - 0x3ffe; 181 exp = dest->exp + src->exp - 0x3ffe;
182 182
183 /* shift up the mantissa for denormali 183 /* shift up the mantissa for denormalized numbers,
184 so that the highest bit is set, thi 184 so that the highest bit is set, this makes the
185 shift of the result below easier */ 185 shift of the result below easier */
186 if ((long)dest->mant.m32[0] >= 0) 186 if ((long)dest->mant.m32[0] >= 0)
187 exp -= fp_overnormalize(dest); 187 exp -= fp_overnormalize(dest);
188 if ((long)src->mant.m32[0] >= 0) 188 if ((long)src->mant.m32[0] >= 0)
189 exp -= fp_overnormalize(src); 189 exp -= fp_overnormalize(src);
190 190
191 /* now, do a 64-bit multiply with expa 191 /* now, do a 64-bit multiply with expansion */
192 fp_multiplymant(&temp, dest, src); 192 fp_multiplymant(&temp, dest, src);
193 193
194 /* normalize it back to 64 bits and st 194 /* normalize it back to 64 bits and stuff it back into the
195 destination struct */ 195 destination struct */
196 if ((long)temp.m32[0] > 0) { 196 if ((long)temp.m32[0] > 0) {
197 exp--; 197 exp--;
198 fp_putmant128(dest, &temp, 1); 198 fp_putmant128(dest, &temp, 1);
199 } else 199 } else
200 fp_putmant128(dest, &temp, 0); 200 fp_putmant128(dest, &temp, 0);
201 201
202 if (exp >= 0x7fff) { 202 if (exp >= 0x7fff) {
203 fp_set_ovrflw(dest); 203 fp_set_ovrflw(dest);
204 return dest; 204 return dest;
205 } 205 }
206 dest->exp = exp; 206 dest->exp = exp;
207 if (exp < 0) { 207 if (exp < 0) {
208 fp_set_sr(FPSR_EXC_UNFL); 208 fp_set_sr(FPSR_EXC_UNFL);
209 fp_denormalize(dest, -exp); 209 fp_denormalize(dest, -exp);
210 } 210 }
211 211
212 return dest; 212 return dest;
213 } 213 }
214 214
215 /* fp_fdiv: Implements the "kernel" of the FDI 215 /* fp_fdiv: Implements the "kernel" of the FDIV, FSDIV, FDDIV and
216 FSGLDIV instructions. 216 FSGLDIV instructions.
217 217
218 Note that the order of the operands is coun 218 Note that the order of the operands is counter-intuitive: instead
219 of src / dest, the result is actually dest 219 of src / dest, the result is actually dest / src. */
220 220
221 struct fp_ext *fp_fdiv(struct fp_ext *dest, st 221 struct fp_ext *fp_fdiv(struct fp_ext *dest, struct fp_ext *src)
222 { 222 {
223 union fp_mant128 temp; 223 union fp_mant128 temp;
224 int exp; 224 int exp;
225 225
226 dprint(PINSTR, "fdiv\n"); 226 dprint(PINSTR, "fdiv\n");
227 227
228 fp_dyadic_check(dest, src); 228 fp_dyadic_check(dest, src);
229 229
230 /* calculate the correct sign now, as 230 /* calculate the correct sign now, as it's necessary for infinities */
231 dest->sign = src->sign ^ dest->sign; 231 dest->sign = src->sign ^ dest->sign;
232 232
233 /* Handle infinities */ 233 /* Handle infinities */
234 if (IS_INF(dest)) { 234 if (IS_INF(dest)) {
235 /* infinity / infinity = NaN ( 235 /* infinity / infinity = NaN (quiet, as always) */
236 if (IS_INF(src)) 236 if (IS_INF(src))
237 fp_set_nan(dest); 237 fp_set_nan(dest);
238 /* infinity / anything else = 238 /* infinity / anything else = infinity (with appropriate sign) */
239 return dest; 239 return dest;
240 } 240 }
241 if (IS_INF(src)) { 241 if (IS_INF(src)) {
242 /* anything / infinity = zero 242 /* anything / infinity = zero (with appropriate sign) */
243 dest->exp = 0; 243 dest->exp = 0;
244 dest->mant.m64 = 0; 244 dest->mant.m64 = 0;
245 dest->lowmant = 0; 245 dest->lowmant = 0;
246 246
247 return dest; 247 return dest;
248 } 248 }
249 249
250 /* zeroes */ 250 /* zeroes */
251 if (IS_ZERO(dest)) { 251 if (IS_ZERO(dest)) {
252 /* zero / zero = NaN */ 252 /* zero / zero = NaN */
253 if (IS_ZERO(src)) 253 if (IS_ZERO(src))
254 fp_set_nan(dest); 254 fp_set_nan(dest);
255 /* zero / anything else = zero 255 /* zero / anything else = zero */
256 return dest; 256 return dest;
257 } 257 }
258 if (IS_ZERO(src)) { 258 if (IS_ZERO(src)) {
259 /* anything / zero = infinity 259 /* anything / zero = infinity (with appropriate sign) */
260 fp_set_sr(FPSR_EXC_DZ); 260 fp_set_sr(FPSR_EXC_DZ);
261 dest->exp = 0x7fff; 261 dest->exp = 0x7fff;
262 dest->mant.m64 = 0; 262 dest->mant.m64 = 0;
263 263
264 return dest; 264 return dest;
265 } 265 }
266 266
267 exp = dest->exp - src->exp + 0x3fff; 267 exp = dest->exp - src->exp + 0x3fff;
268 268
269 /* shift up the mantissa for denormali 269 /* shift up the mantissa for denormalized numbers,
270 so that the highest bit is set, thi 270 so that the highest bit is set, this makes lots
271 of things below easier */ 271 of things below easier */
272 if ((long)dest->mant.m32[0] >= 0) 272 if ((long)dest->mant.m32[0] >= 0)
273 exp -= fp_overnormalize(dest); 273 exp -= fp_overnormalize(dest);
274 if ((long)src->mant.m32[0] >= 0) 274 if ((long)src->mant.m32[0] >= 0)
275 exp -= fp_overnormalize(src); 275 exp -= fp_overnormalize(src);
276 276
277 /* now, do the 64-bit divide */ 277 /* now, do the 64-bit divide */
278 fp_dividemant(&temp, dest, src); 278 fp_dividemant(&temp, dest, src);
279 279
280 /* normalize it back to 64 bits and st 280 /* normalize it back to 64 bits and stuff it back into the
281 destination struct */ 281 destination struct */
282 if (!temp.m32[0]) { 282 if (!temp.m32[0]) {
283 exp--; 283 exp--;
284 fp_putmant128(dest, &temp, 32) 284 fp_putmant128(dest, &temp, 32);
285 } else 285 } else
286 fp_putmant128(dest, &temp, 31) 286 fp_putmant128(dest, &temp, 31);
287 287
288 if (exp >= 0x7fff) { 288 if (exp >= 0x7fff) {
289 fp_set_ovrflw(dest); 289 fp_set_ovrflw(dest);
290 return dest; 290 return dest;
291 } 291 }
292 dest->exp = exp; 292 dest->exp = exp;
293 if (exp < 0) { 293 if (exp < 0) {
294 fp_set_sr(FPSR_EXC_UNFL); 294 fp_set_sr(FPSR_EXC_UNFL);
295 fp_denormalize(dest, -exp); 295 fp_denormalize(dest, -exp);
296 } 296 }
297 297
298 return dest; 298 return dest;
299 } 299 }
300 300
301 struct fp_ext *fp_fsglmul(struct fp_ext *dest, 301 struct fp_ext *fp_fsglmul(struct fp_ext *dest, struct fp_ext *src)
302 { 302 {
303 int exp; 303 int exp;
304 304
305 dprint(PINSTR, "fsglmul\n"); 305 dprint(PINSTR, "fsglmul\n");
306 306
307 fp_dyadic_check(dest, src); 307 fp_dyadic_check(dest, src);
308 308
309 /* calculate the correct sign now, as 309 /* calculate the correct sign now, as it's necessary for infinities */
310 dest->sign = src->sign ^ dest->sign; 310 dest->sign = src->sign ^ dest->sign;
311 311
312 /* Handle infinities */ 312 /* Handle infinities */
313 if (IS_INF(dest)) { 313 if (IS_INF(dest)) {
314 if (IS_ZERO(src)) 314 if (IS_ZERO(src))
315 fp_set_nan(dest); 315 fp_set_nan(dest);
316 return dest; 316 return dest;
317 } 317 }
318 if (IS_INF(src)) { 318 if (IS_INF(src)) {
319 if (IS_ZERO(dest)) 319 if (IS_ZERO(dest))
320 fp_set_nan(dest); 320 fp_set_nan(dest);
321 else 321 else
322 fp_copy_ext(dest, src) 322 fp_copy_ext(dest, src);
323 return dest; 323 return dest;
324 } 324 }
325 325
326 /* Of course, as we all know, zero * a 326 /* Of course, as we all know, zero * anything = zero. You may
327 not have known that it might be a p 327 not have known that it might be a positive or negative
328 zero... */ 328 zero... */
329 if (IS_ZERO(dest) || IS_ZERO(src)) { 329 if (IS_ZERO(dest) || IS_ZERO(src)) {
330 dest->exp = 0; 330 dest->exp = 0;
331 dest->mant.m64 = 0; 331 dest->mant.m64 = 0;
332 dest->lowmant = 0; 332 dest->lowmant = 0;
333 333
334 return dest; 334 return dest;
335 } 335 }
336 336
337 exp = dest->exp + src->exp - 0x3ffe; 337 exp = dest->exp + src->exp - 0x3ffe;
338 338
339 /* do a 32-bit multiply */ 339 /* do a 32-bit multiply */
340 fp_mul64(dest->mant.m32[0], dest->mant 340 fp_mul64(dest->mant.m32[0], dest->mant.m32[1],
341 dest->mant.m32[0] & 0xffffff0 341 dest->mant.m32[0] & 0xffffff00,
342 src->mant.m32[0] & 0xffffff00 342 src->mant.m32[0] & 0xffffff00);
343 343
344 if (exp >= 0x7fff) { 344 if (exp >= 0x7fff) {
345 fp_set_ovrflw(dest); 345 fp_set_ovrflw(dest);
346 return dest; 346 return dest;
347 } 347 }
348 dest->exp = exp; 348 dest->exp = exp;
349 if (exp < 0) { 349 if (exp < 0) {
350 fp_set_sr(FPSR_EXC_UNFL); 350 fp_set_sr(FPSR_EXC_UNFL);
351 fp_denormalize(dest, -exp); 351 fp_denormalize(dest, -exp);
352 } 352 }
353 353
354 return dest; 354 return dest;
355 } 355 }
356 356
357 struct fp_ext *fp_fsgldiv(struct fp_ext *dest, 357 struct fp_ext *fp_fsgldiv(struct fp_ext *dest, struct fp_ext *src)
358 { 358 {
359 int exp; 359 int exp;
360 unsigned long quot, rem; 360 unsigned long quot, rem;
361 361
362 dprint(PINSTR, "fsgldiv\n"); 362 dprint(PINSTR, "fsgldiv\n");
363 363
364 fp_dyadic_check(dest, src); 364 fp_dyadic_check(dest, src);
365 365
366 /* calculate the correct sign now, as 366 /* calculate the correct sign now, as it's necessary for infinities */
367 dest->sign = src->sign ^ dest->sign; 367 dest->sign = src->sign ^ dest->sign;
368 368
369 /* Handle infinities */ 369 /* Handle infinities */
370 if (IS_INF(dest)) { 370 if (IS_INF(dest)) {
371 /* infinity / infinity = NaN ( 371 /* infinity / infinity = NaN (quiet, as always) */
372 if (IS_INF(src)) 372 if (IS_INF(src))
373 fp_set_nan(dest); 373 fp_set_nan(dest);
374 /* infinity / anything else = 374 /* infinity / anything else = infinity (with approprate sign) */
375 return dest; 375 return dest;
376 } 376 }
377 if (IS_INF(src)) { 377 if (IS_INF(src)) {
378 /* anything / infinity = zero 378 /* anything / infinity = zero (with appropriate sign) */
379 dest->exp = 0; 379 dest->exp = 0;
380 dest->mant.m64 = 0; 380 dest->mant.m64 = 0;
381 dest->lowmant = 0; 381 dest->lowmant = 0;
382 382
383 return dest; 383 return dest;
384 } 384 }
385 385
386 /* zeroes */ 386 /* zeroes */
387 if (IS_ZERO(dest)) { 387 if (IS_ZERO(dest)) {
388 /* zero / zero = NaN */ 388 /* zero / zero = NaN */
389 if (IS_ZERO(src)) 389 if (IS_ZERO(src))
390 fp_set_nan(dest); 390 fp_set_nan(dest);
391 /* zero / anything else = zero 391 /* zero / anything else = zero */
392 return dest; 392 return dest;
393 } 393 }
394 if (IS_ZERO(src)) { 394 if (IS_ZERO(src)) {
395 /* anything / zero = infinity 395 /* anything / zero = infinity (with appropriate sign) */
396 fp_set_sr(FPSR_EXC_DZ); 396 fp_set_sr(FPSR_EXC_DZ);
397 dest->exp = 0x7fff; 397 dest->exp = 0x7fff;
398 dest->mant.m64 = 0; 398 dest->mant.m64 = 0;
399 399
400 return dest; 400 return dest;
401 } 401 }
402 402
403 exp = dest->exp - src->exp + 0x3fff; 403 exp = dest->exp - src->exp + 0x3fff;
404 404
405 dest->mant.m32[0] &= 0xffffff00; 405 dest->mant.m32[0] &= 0xffffff00;
406 src->mant.m32[0] &= 0xffffff00; 406 src->mant.m32[0] &= 0xffffff00;
407 407
408 /* do the 32-bit divide */ 408 /* do the 32-bit divide */
409 if (dest->mant.m32[0] >= src->mant.m32 409 if (dest->mant.m32[0] >= src->mant.m32[0]) {
410 fp_sub64(dest->mant, src->mant 410 fp_sub64(dest->mant, src->mant);
411 fp_div64(quot, rem, dest->mant 411 fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
412 dest->mant.m32[0] = 0x80000000 412 dest->mant.m32[0] = 0x80000000 | (quot >> 1);
413 dest->mant.m32[1] = (quot & 1) 413 dest->mant.m32[1] = (quot & 1) | rem; /* only for rounding */
414 } else { 414 } else {
415 fp_div64(quot, rem, dest->mant 415 fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
416 dest->mant.m32[0] = quot; 416 dest->mant.m32[0] = quot;
417 dest->mant.m32[1] = rem; 417 dest->mant.m32[1] = rem; /* only for rounding */
418 exp--; 418 exp--;
419 } 419 }
420 420
421 if (exp >= 0x7fff) { 421 if (exp >= 0x7fff) {
422 fp_set_ovrflw(dest); 422 fp_set_ovrflw(dest);
423 return dest; 423 return dest;
424 } 424 }
425 dest->exp = exp; 425 dest->exp = exp;
426 if (exp < 0) { 426 if (exp < 0) {
427 fp_set_sr(FPSR_EXC_UNFL); 427 fp_set_sr(FPSR_EXC_UNFL);
428 fp_denormalize(dest, -exp); 428 fp_denormalize(dest, -exp);
429 } 429 }
430 430
431 return dest; 431 return dest;
432 } 432 }
433 433
434 /* fp_roundint: Internal rounding function for 434 /* fp_roundint: Internal rounding function for use by several of these
435 emulated instructions. 435 emulated instructions.
436 436
437 This one rounds off the fractional part usi 437 This one rounds off the fractional part using the rounding mode
438 specified. */ 438 specified. */
439 439
440 static void fp_roundint(struct fp_ext *dest, i 440 static void fp_roundint(struct fp_ext *dest, int mode)
441 { 441 {
442 union fp_mant64 oldmant; 442 union fp_mant64 oldmant;
443 unsigned long mask; 443 unsigned long mask;
444 444
445 if (!fp_normalize_ext(dest)) 445 if (!fp_normalize_ext(dest))
446 return; 446 return;
447 447
448 /* infinities and zeroes */ 448 /* infinities and zeroes */
449 if (IS_INF(dest) || IS_ZERO(dest)) 449 if (IS_INF(dest) || IS_ZERO(dest))
450 return; 450 return;
451 451
452 /* first truncate the lower bits */ 452 /* first truncate the lower bits */
453 oldmant = dest->mant; 453 oldmant = dest->mant;
454 switch (dest->exp) { 454 switch (dest->exp) {
455 case 0 ... 0x3ffe: 455 case 0 ... 0x3ffe:
456 dest->mant.m64 = 0; 456 dest->mant.m64 = 0;
457 break; 457 break;
458 case 0x3fff ... 0x401e: 458 case 0x3fff ... 0x401e:
459 dest->mant.m32[0] &= 0xfffffff 459 dest->mant.m32[0] &= 0xffffffffU << (0x401e - dest->exp);
460 dest->mant.m32[1] = 0; 460 dest->mant.m32[1] = 0;
461 if (oldmant.m64 == dest->mant. 461 if (oldmant.m64 == dest->mant.m64)
462 return; 462 return;
463 break; 463 break;
464 case 0x401f ... 0x403e: 464 case 0x401f ... 0x403e:
465 dest->mant.m32[1] &= 0xfffffff 465 dest->mant.m32[1] &= 0xffffffffU << (0x403e - dest->exp);
466 if (oldmant.m32[1] == dest->ma 466 if (oldmant.m32[1] == dest->mant.m32[1])
467 return; 467 return;
468 break; 468 break;
469 default: 469 default:
470 return; 470 return;
471 } 471 }
472 fp_set_sr(FPSR_EXC_INEX2); 472 fp_set_sr(FPSR_EXC_INEX2);
473 473
474 /* We might want to normalize upwards 474 /* We might want to normalize upwards here... however, since
475 we know that this is only called on 475 we know that this is only called on the output of fp_fdiv,
476 or with the input to fp_fint or fp_ 476 or with the input to fp_fint or fp_fintrz, and the inputs
477 to all these functions are either n 477 to all these functions are either normal or denormalized
478 (no subnormals allowed!), there's r 478 (no subnormals allowed!), there's really no need.
479 479
480 In the case of fp_fdiv, observe tha 480 In the case of fp_fdiv, observe that 0x80000000 / 0xffff =
481 0xffff8000, and the same holds for 481 0xffff8000, and the same holds for 128-bit / 64-bit. (i.e. the
482 smallest possible normal dividend a 482 smallest possible normal dividend and the largest possible normal
483 divisor will still produce a normal 483 divisor will still produce a normal quotient, therefore, (normal
484 << 64) / normal is normal in all ca 484 << 64) / normal is normal in all cases) */
485 485
486 switch (mode) { 486 switch (mode) {
487 case FPCR_ROUND_RN: 487 case FPCR_ROUND_RN:
488 switch (dest->exp) { 488 switch (dest->exp) {
489 case 0 ... 0x3ffd: 489 case 0 ... 0x3ffd:
490 return; 490 return;
491 case 0x3ffe: 491 case 0x3ffe:
492 /* As noted above, the 492 /* As noted above, the input is always normal, so the
493 guard bit (bit 63) 493 guard bit (bit 63) is always set. therefore, the
494 only case in which 494 only case in which we will NOT round to 1.0 is when
495 the input is exactl 495 the input is exactly 0.5. */
496 if (oldmant.m64 == (1U 496 if (oldmant.m64 == (1ULL << 63))
497 return; 497 return;
498 break; 498 break;
499 case 0x3fff ... 0x401d: 499 case 0x3fff ... 0x401d:
500 mask = 1 << (0x401d - 500 mask = 1 << (0x401d - dest->exp);
501 if (!(oldmant.m32[0] & 501 if (!(oldmant.m32[0] & mask))
502 return; 502 return;
503 if (oldmant.m32[0] & ( 503 if (oldmant.m32[0] & (mask << 1))
504 break; 504 break;
505 if (!(oldmant.m32[0] < 505 if (!(oldmant.m32[0] << (dest->exp - 0x3ffd)) &&
506 !oldma 506 !oldmant.m32[1])
507 return; 507 return;
508 break; 508 break;
509 case 0x401e: 509 case 0x401e:
510 if (oldmant.m32[1] & 0 510 if (oldmant.m32[1] & 0x80000000)
511 return; 511 return;
512 if (oldmant.m32[0] & 1 512 if (oldmant.m32[0] & 1)
513 break; 513 break;
514 if (!(oldmant.m32[1] < 514 if (!(oldmant.m32[1] << 1))
515 return; 515 return;
516 break; 516 break;
517 case 0x401f ... 0x403d: 517 case 0x401f ... 0x403d:
518 mask = 1 << (0x403d - 518 mask = 1 << (0x403d - dest->exp);
519 if (!(oldmant.m32[1] & 519 if (!(oldmant.m32[1] & mask))
520 return; 520 return;
521 if (oldmant.m32[1] & ( 521 if (oldmant.m32[1] & (mask << 1))
522 break; 522 break;
523 if (!(oldmant.m32[1] < 523 if (!(oldmant.m32[1] << (dest->exp - 0x401d)))
524 return; 524 return;
525 break; 525 break;
526 default: 526 default:
527 return; 527 return;
528 } 528 }
529 break; 529 break;
530 case FPCR_ROUND_RZ: 530 case FPCR_ROUND_RZ:
531 return; 531 return;
532 default: 532 default:
533 if (dest->sign ^ (mode - FPCR_ 533 if (dest->sign ^ (mode - FPCR_ROUND_RM))
534 break; 534 break;
535 return; 535 return;
536 } 536 }
537 537
538 switch (dest->exp) { 538 switch (dest->exp) {
539 case 0 ... 0x3ffe: 539 case 0 ... 0x3ffe:
540 dest->exp = 0x3fff; 540 dest->exp = 0x3fff;
541 dest->mant.m64 = 1ULL << 63; 541 dest->mant.m64 = 1ULL << 63;
542 break; 542 break;
543 case 0x3fff ... 0x401e: 543 case 0x3fff ... 0x401e:
544 mask = 1 << (0x401e - dest->ex 544 mask = 1 << (0x401e - dest->exp);
545 if (dest->mant.m32[0] += mask) 545 if (dest->mant.m32[0] += mask)
546 break; 546 break;
547 dest->mant.m32[0] = 0x80000000 547 dest->mant.m32[0] = 0x80000000;
548 dest->exp++; 548 dest->exp++;
549 break; 549 break;
550 case 0x401f ... 0x403e: 550 case 0x401f ... 0x403e:
551 mask = 1 << (0x403e - dest->ex 551 mask = 1 << (0x403e - dest->exp);
552 if (dest->mant.m32[1] += mask) 552 if (dest->mant.m32[1] += mask)
553 break; 553 break;
554 if (dest->mant.m32[0] += 1) 554 if (dest->mant.m32[0] += 1)
555 break; 555 break;
556 dest->mant.m32[0] = 0x80000000 556 dest->mant.m32[0] = 0x80000000;
557 dest->exp++; 557 dest->exp++;
558 break; 558 break;
559 } 559 }
560 } 560 }
561 561
562 /* modrem_kernel: Implementation of the FREM a 562 /* modrem_kernel: Implementation of the FREM and FMOD instructions
563 (which are exactly the same, except for the 563 (which are exactly the same, except for the rounding used on the
564 intermediate value) */ 564 intermediate value) */
565 565
566 static struct fp_ext *modrem_kernel(struct fp_ 566 static struct fp_ext *modrem_kernel(struct fp_ext *dest, struct fp_ext *src,
567 int mode) 567 int mode)
568 { 568 {
569 struct fp_ext tmp; 569 struct fp_ext tmp;
570 570
571 fp_dyadic_check(dest, src); 571 fp_dyadic_check(dest, src);
572 572
573 /* Infinities and zeros */ 573 /* Infinities and zeros */
574 if (IS_INF(dest) || IS_ZERO(src)) { 574 if (IS_INF(dest) || IS_ZERO(src)) {
575 fp_set_nan(dest); 575 fp_set_nan(dest);
576 return dest; 576 return dest;
577 } 577 }
578 if (IS_ZERO(dest) || IS_INF(src)) 578 if (IS_ZERO(dest) || IS_INF(src))
579 return dest; 579 return dest;
580 580
581 /* FIXME: there is almost certainly a 581 /* FIXME: there is almost certainly a smarter way to do this */
582 fp_copy_ext(&tmp, dest); 582 fp_copy_ext(&tmp, dest);
583 fp_fdiv(&tmp, src); /* NOT 583 fp_fdiv(&tmp, src); /* NOTE: src might be modified */
584 fp_roundint(&tmp, mode); 584 fp_roundint(&tmp, mode);
585 fp_fmul(&tmp, src); 585 fp_fmul(&tmp, src);
586 fp_fsub(dest, &tmp); 586 fp_fsub(dest, &tmp);
587 587
588 /* set the quotient byte */ 588 /* set the quotient byte */
589 fp_set_quotient((dest->mant.m64 & 0x7f 589 fp_set_quotient((dest->mant.m64 & 0x7f) | (dest->sign << 7));
590 return dest; 590 return dest;
591 } 591 }
592 592
593 /* fp_fmod: Implements the kernel of the FMOD 593 /* fp_fmod: Implements the kernel of the FMOD instruction.
594 594
595 Again, the argument order is backwards. Th 595 Again, the argument order is backwards. The result, as defined in
596 the Motorola manuals, is: 596 the Motorola manuals, is:
597 597
598 fmod(src,dest) = (dest - (src * floor(dest 598 fmod(src,dest) = (dest - (src * floor(dest / src))) */
599 599
600 struct fp_ext *fp_fmod(struct fp_ext *dest, st 600 struct fp_ext *fp_fmod(struct fp_ext *dest, struct fp_ext *src)
601 { 601 {
602 dprint(PINSTR, "fmod\n"); 602 dprint(PINSTR, "fmod\n");
603 return modrem_kernel(dest, src, FPCR_R 603 return modrem_kernel(dest, src, FPCR_ROUND_RZ);
604 } 604 }
605 605
606 /* fp_frem: Implements the kernel of the FREM 606 /* fp_frem: Implements the kernel of the FREM instruction.
607 607
608 frem(src,dest) = (dest - (src * round(dest 608 frem(src,dest) = (dest - (src * round(dest / src)))
609 */ 609 */
610 610
611 struct fp_ext *fp_frem(struct fp_ext *dest, st 611 struct fp_ext *fp_frem(struct fp_ext *dest, struct fp_ext *src)
612 { 612 {
613 dprint(PINSTR, "frem\n"); 613 dprint(PINSTR, "frem\n");
614 return modrem_kernel(dest, src, FPCR_R 614 return modrem_kernel(dest, src, FPCR_ROUND_RN);
615 } 615 }
616 616
617 struct fp_ext *fp_fint(struct fp_ext *dest, st 617 struct fp_ext *fp_fint(struct fp_ext *dest, struct fp_ext *src)
618 { 618 {
619 dprint(PINSTR, "fint\n"); 619 dprint(PINSTR, "fint\n");
620 620
621 fp_copy_ext(dest, src); 621 fp_copy_ext(dest, src);
622 622
623 fp_roundint(dest, FPDATA->rnd); 623 fp_roundint(dest, FPDATA->rnd);
624 624
625 return dest; 625 return dest;
626 } 626 }
627 627
628 struct fp_ext *fp_fintrz(struct fp_ext *dest, 628 struct fp_ext *fp_fintrz(struct fp_ext *dest, struct fp_ext *src)
629 { 629 {
630 dprint(PINSTR, "fintrz\n"); 630 dprint(PINSTR, "fintrz\n");
631 631
632 fp_copy_ext(dest, src); 632 fp_copy_ext(dest, src);
633 633
634 fp_roundint(dest, FPCR_ROUND_RZ); 634 fp_roundint(dest, FPCR_ROUND_RZ);
635 635
636 return dest; 636 return dest;
637 } 637 }
638 638
639 struct fp_ext *fp_fscale(struct fp_ext *dest, 639 struct fp_ext *fp_fscale(struct fp_ext *dest, struct fp_ext *src)
640 { 640 {
641 int scale, oldround; 641 int scale, oldround;
642 642
643 dprint(PINSTR, "fscale\n"); 643 dprint(PINSTR, "fscale\n");
644 644
645 fp_dyadic_check(dest, src); 645 fp_dyadic_check(dest, src);
646 646
647 /* Infinities */ 647 /* Infinities */
648 if (IS_INF(src)) { 648 if (IS_INF(src)) {
649 fp_set_nan(dest); 649 fp_set_nan(dest);
650 return dest; 650 return dest;
651 } 651 }
652 if (IS_INF(dest)) 652 if (IS_INF(dest))
653 return dest; 653 return dest;
654 654
655 /* zeroes */ 655 /* zeroes */
656 if (IS_ZERO(src) || IS_ZERO(dest)) 656 if (IS_ZERO(src) || IS_ZERO(dest))
657 return dest; 657 return dest;
658 658
659 /* Source exponent out of range */ 659 /* Source exponent out of range */
660 if (src->exp >= 0x400c) { 660 if (src->exp >= 0x400c) {
661 fp_set_ovrflw(dest); 661 fp_set_ovrflw(dest);
662 return dest; 662 return dest;
663 } 663 }
664 664
665 /* src must be rounded with round to z 665 /* src must be rounded with round to zero. */
666 oldround = FPDATA->rnd; 666 oldround = FPDATA->rnd;
667 FPDATA->rnd = FPCR_ROUND_RZ; 667 FPDATA->rnd = FPCR_ROUND_RZ;
668 scale = fp_conv_ext2long(src); 668 scale = fp_conv_ext2long(src);
669 FPDATA->rnd = oldround; 669 FPDATA->rnd = oldround;
670 670
671 /* new exponent */ 671 /* new exponent */
672 scale += dest->exp; 672 scale += dest->exp;
673 673
674 if (scale >= 0x7fff) { 674 if (scale >= 0x7fff) {
675 fp_set_ovrflw(dest); 675 fp_set_ovrflw(dest);
676 } else if (scale <= 0) { 676 } else if (scale <= 0) {
677 fp_set_sr(FPSR_EXC_UNFL); 677 fp_set_sr(FPSR_EXC_UNFL);
678 fp_denormalize(dest, -scale); 678 fp_denormalize(dest, -scale);
679 } else 679 } else
680 dest->exp = scale; 680 dest->exp = scale;
681 681
682 return dest; 682 return dest;
683 } 683 }
684 684
685 685
Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.