1 /* $OpenBSD: fpu_implode.c,v 1.9 2022/10/16 01:22:39 jsg Exp $ */ 2 /* $NetBSD: fpu_implode.c,v 1.7 2000/08/03 18:32:08 eeh Exp $ */ 3 4 /* 5 * Copyright (c) 1992, 1993 6 * The Regents of the University of California. All rights reserved. 7 * 8 * This software was developed by the Computer Systems Engineering group 9 * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and 10 * contributed to Berkeley. 11 * 12 * All advertising materials mentioning features or use of this software 13 * must display the following acknowledgement: 14 * This product includes software developed by the University of 15 * California, Lawrence Berkeley Laboratory. 16 * 17 * Redistribution and use in source and binary forms, with or without 18 * modification, are permitted provided that the following conditions 19 * are met: 20 * 1. Redistributions of source code must retain the above copyright 21 * notice, this list of conditions and the following disclaimer. 22 * 2. Redistributions in binary form must reproduce the above copyright 23 * notice, this list of conditions and the following disclaimer in the 24 * documentation and/or other materials provided with the distribution. 25 * 3. Neither the name of the University nor the names of its contributors 26 * may be used to endorse or promote products derived from this software 27 * without specific prior written permission. 28 * 29 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 30 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 31 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 32 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 33 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 34 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 35 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 36 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 37 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 38 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 39 * SUCH DAMAGE. 40 * 41 * @(#)fpu_implode.c 8.1 (Berkeley) 6/11/93 42 */ 43 44 /* 45 * FPU subroutines: `implode' internal format numbers into the machine's 46 * `packed binary' format. 47 */ 48 49 #include <sys/types.h> 50 #include <sys/systm.h> 51 52 #include <machine/ieee.h> 53 #include <machine/instr.h> 54 #include <machine/reg.h> 55 56 #include <sparc64/fpu/fpu_arith.h> 57 #include <sparc64/fpu/fpu_emu.h> 58 #include <sparc64/fpu/fpu_extern.h> 59 60 static int fpu_round(register struct fpemu *, register struct fpn *); 61 static int toinf(struct fpemu *, int); 62 63 /* 64 * Round a number (algorithm from Motorola MC68882 manual, modified for 65 * our internal format). Set inexact exception if rounding is required. 66 * Return true iff we rounded up. 67 * 68 * After rounding, we discard the guard and round bits by shifting right 69 * 2 bits (a la fpu_shr(), but we do not bother with fp->fp_sticky). 70 * This saves effort later. 71 * 72 * Note that we may leave the value 2.0 in fp->fp_mant; it is the caller's 73 * responsibility to fix this if necessary. 74 */ 75 static int 76 fpu_round(register struct fpemu *fe, register struct fpn *fp) 77 { 78 register u_int m0, m1, m2, m3; 79 register int gr, s; 80 81 m0 = fp->fp_mant[0]; 82 m1 = fp->fp_mant[1]; 83 m2 = fp->fp_mant[2]; 84 m3 = fp->fp_mant[3]; 85 gr = m3 & 3; 86 s = fp->fp_sticky; 87 88 /* mant >>= FP_NG */ 89 m3 = (m3 >> FP_NG) | (m2 << (32 - FP_NG)); 90 m2 = (m2 >> FP_NG) | (m1 << (32 - FP_NG)); 91 m1 = (m1 >> FP_NG) | (m0 << (32 - FP_NG)); 92 m0 >>= FP_NG; 93 94 if ((gr | s) == 0) /* result is exact: no rounding needed */ 95 goto rounddown; 96 97 fe->fe_cx |= FSR_NX; /* inexact */ 98 99 /* Go to rounddown to round down; break to round up. */ 100 switch ((fe->fe_fsr >> FSR_RD_SHIFT) & FSR_RD_MASK) { 101 102 case FSR_RD_RN: 103 default: 104 /* 105 * Round only if guard is set (gr & 2). If guard is set, 106 * but round & sticky both clear, then we want to round 107 * but have a tie, so round to even, i.e., add 1 iff odd. 108 */ 109 if ((gr & 2) == 0) 110 goto rounddown; 111 if ((gr & 1) || fp->fp_sticky || (m3 & 1)) 112 break; 113 goto rounddown; 114 115 case FSR_RD_RZ: 116 /* Round towards zero, i.e., down. */ 117 goto rounddown; 118 119 case FSR_RD_RM: 120 /* Round towards -Inf: up if negative, down if positive. */ 121 if (fp->fp_sign) 122 break; 123 goto rounddown; 124 125 case FSR_RD_RP: 126 /* Round towards +Inf: up if positive, down otherwise. */ 127 if (!fp->fp_sign) 128 break; 129 goto rounddown; 130 } 131 132 /* Bump low bit of mantissa, with carry. */ 133 FPU_ADDS(m3, m3, 1); 134 FPU_ADDCS(m2, m2, 0); 135 FPU_ADDCS(m1, m1, 0); 136 FPU_ADDC(m0, m0, 0); 137 fp->fp_mant[0] = m0; 138 fp->fp_mant[1] = m1; 139 fp->fp_mant[2] = m2; 140 fp->fp_mant[3] = m3; 141 return (1); 142 143 rounddown: 144 fp->fp_mant[0] = m0; 145 fp->fp_mant[1] = m1; 146 fp->fp_mant[2] = m2; 147 fp->fp_mant[3] = m3; 148 return (0); 149 } 150 151 /* 152 * For overflow: return true if overflow is to go to +/-Inf, according 153 * to the sign of the overflowing result. If false, overflow is to go 154 * to the largest magnitude value instead. 155 */ 156 static int 157 toinf(struct fpemu *fe, int sign) 158 { 159 int inf; 160 161 /* look at rounding direction */ 162 switch ((fe->fe_fsr >> FSR_RD_SHIFT) & FSR_RD_MASK) { 163 164 default: 165 case FSR_RD_RN: /* the nearest value is always Inf */ 166 inf = 1; 167 break; 168 169 case FSR_RD_RZ: /* toward 0 => never towards Inf */ 170 inf = 0; 171 break; 172 173 case FSR_RD_RP: /* toward +Inf iff positive */ 174 inf = sign == 0; 175 break; 176 177 case FSR_RD_RM: /* toward -Inf iff negative */ 178 inf = sign; 179 break; 180 } 181 return (inf); 182 } 183 184 /* 185 * fpn -> int (int value returned as return value). 186 * 187 * N.B.: this conversion always rounds towards zero (this is a peculiarity 188 * of the SPARC instruction set). 189 */ 190 u_int 191 fpu_ftoi(struct fpemu *fe, register struct fpn *fp) 192 { 193 register u_int i; 194 register int sign, exp; 195 196 sign = fp->fp_sign; 197 switch (fp->fp_class) { 198 199 case FPC_ZERO: 200 return (0); 201 202 case FPC_NUM: 203 /* 204 * If exp >= 2^32, overflow. Otherwise shift value right 205 * into last mantissa word (this will not exceed 0xffffffff), 206 * shifting any guard and round bits out into the sticky 207 * bit. Then ``round'' towards zero, i.e., just set an 208 * inexact exception if sticky is set (see fpu_round()). 209 * If the result is > 0x80000000, or is positive and equals 210 * 0x80000000, overflow; otherwise the last fraction word 211 * is the result. 212 */ 213 if ((exp = fp->fp_exp) >= 32) 214 break; 215 /* NB: the following includes exp < 0 cases */ 216 if (fpu_shr(fp, FP_NMANT - 1 - exp) != 0) 217 fe->fe_cx |= FSR_NX; 218 i = fp->fp_mant[3]; 219 if (i >= ((u_int)0x80000000 + sign)) 220 break; 221 return (sign ? -i : i); 222 223 default: /* Inf, qNaN, sNaN */ 224 break; 225 } 226 /* overflow: replace any inexact exception with invalid */ 227 fe->fe_cx = (fe->fe_cx & ~FSR_NX) | FSR_NV; 228 return (0x7fffffff + sign); 229 } 230 231 /* 232 * fpn -> extended int (high bits of int value returned as return value). 233 * 234 * N.B.: this conversion always rounds towards zero (this is a peculiarity 235 * of the SPARC instruction set). 236 */ 237 u_int 238 fpu_ftox(struct fpemu *fe, register struct fpn *fp, u_int *res) 239 { 240 register u_int64_t i; 241 register int sign, exp; 242 243 sign = fp->fp_sign; 244 switch (fp->fp_class) { 245 246 case FPC_ZERO: 247 i = 0; 248 goto out; 249 250 case FPC_NUM: 251 /* 252 * If exp >= 2^64, overflow. Otherwise shift value right 253 * into last mantissa word (this will not exceed 0xffffffffffffffff), 254 * shifting any guard and round bits out into the sticky 255 * bit. Then ``round'' towards zero, i.e., just set an 256 * inexact exception if sticky is set (see fpu_round()). 257 * If the result is > 0x8000000000000000, or is positive and equals 258 * 0x8000000000000000, overflow; otherwise the last fraction word 259 * is the result. 260 */ 261 if ((exp = fp->fp_exp) >= 64) 262 break; 263 /* NB: the following includes exp < 0 cases */ 264 if (fpu_shr(fp, FP_NMANT - 1 - exp) != 0) 265 fe->fe_cx |= FSR_NX; 266 i = ((u_int64_t)fp->fp_mant[2]<<32)|fp->fp_mant[3]; 267 if (i >= ((u_int64_t)0x8000000000000000LL + sign)) 268 break; 269 if (sign) 270 i = -i; 271 goto out; 272 273 default: /* Inf, qNaN, sNaN */ 274 break; 275 } 276 /* overflow: replace any inexact exception with invalid */ 277 fe->fe_cx = (fe->fe_cx & ~FSR_NX) | FSR_NV; 278 i = 0x7fffffffffffffffLL + sign; 279 out: 280 res[1] = i & 0xffffffff; 281 return (i >> 32); 282 } 283 284 /* 285 * fpn -> single (32 bit single returned as return value). 286 * We assume <= 29 bits in a single-precision fraction (1.f part). 287 */ 288 u_int 289 fpu_ftos(struct fpemu *fe, register struct fpn *fp) 290 { 291 register u_int sign = fp->fp_sign << 31; 292 register int exp; 293 294 #define SNG_EXP(e) ((e) << SNG_FRACBITS) /* makes e an exponent */ 295 #define SNG_MASK (SNG_EXP(1) - 1) /* mask for fraction */ 296 297 /* Take care of non-numbers first. */ 298 if (ISNAN(fp)) { 299 /* 300 * Preserve upper bits of NaN, per SPARC V8 appendix N. 301 * Note that fp->fp_mant[0] has the quiet bit set, 302 * even if it is classified as a signalling NaN. 303 */ 304 (void) fpu_shr(fp, FP_NMANT - 1 - SNG_FRACBITS); 305 exp = SNG_EXP_INFNAN; 306 goto done; 307 } 308 if (ISINF(fp)) 309 return (sign | SNG_EXP(SNG_EXP_INFNAN)); 310 if (ISZERO(fp)) 311 return (sign); 312 313 /* 314 * Normals (including subnormals). Drop all the fraction bits 315 * (including the explicit ``implied'' 1 bit) down into the 316 * single-precision range. If the number is subnormal, move 317 * the ``implied'' 1 into the explicit range as well, and shift 318 * right to introduce leading zeroes. Rounding then acts 319 * differently for normals and subnormals: the largest subnormal 320 * may round to the smallest normal (1.0 x 2^minexp), or may 321 * remain subnormal. In the latter case, signal an underflow 322 * if the result was inexact or if underflow traps are enabled. 323 * 324 * Rounding a normal, on the other hand, always produces another 325 * normal (although either way the result might be too big for 326 * single precision, and cause an overflow). If rounding a 327 * normal produces 2.0 in the fraction, we need not adjust that 328 * fraction at all, since both 1.0 and 2.0 are zero under the 329 * fraction mask. 330 * 331 * Note that the guard and round bits vanish from the number after 332 * rounding. 333 */ 334 if ((exp = fp->fp_exp + SNG_EXP_BIAS) <= 0) { /* subnormal */ 335 /* -NG for g,r; -SNG_FRACBITS-exp for fraction */ 336 (void) fpu_shr(fp, FP_NMANT - FP_NG - SNG_FRACBITS - exp); 337 if (fpu_round(fe, fp) && fp->fp_mant[3] == SNG_EXP(1)) 338 return (sign | SNG_EXP(1) | 0); 339 if ((fe->fe_cx & FSR_NX) || 340 (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT))) 341 fe->fe_cx |= FSR_UF; 342 return (sign | SNG_EXP(0) | fp->fp_mant[3]); 343 } 344 /* -FP_NG for g,r; -1 for implied 1; -SNG_FRACBITS for fraction */ 345 (void) fpu_shr(fp, FP_NMANT - FP_NG - 1 - SNG_FRACBITS); 346 #ifdef DIAGNOSTIC 347 if ((fp->fp_mant[3] & SNG_EXP(1 << FP_NG)) == 0) 348 panic("fpu_ftos"); 349 #endif 350 if (fpu_round(fe, fp) && fp->fp_mant[3] == SNG_EXP(2)) 351 exp++; 352 if (exp >= SNG_EXP_INFNAN) { 353 /* overflow to inf or to max single */ 354 fe->fe_cx |= FSR_OF | FSR_NX; 355 if (toinf(fe, sign)) 356 return (sign | SNG_EXP(SNG_EXP_INFNAN)); 357 return (sign | SNG_EXP(SNG_EXP_INFNAN - 1) | SNG_MASK); 358 } 359 done: 360 /* phew, made it */ 361 return (sign | SNG_EXP(exp) | (fp->fp_mant[3] & SNG_MASK)); 362 } 363 364 /* 365 * fpn -> double (32 bit high-order result returned; 32-bit low order result 366 * left in res[1]). Assumes <= 61 bits in double precision fraction. 367 * 368 * This code mimics fpu_ftos; see it for comments. 369 */ 370 u_int 371 fpu_ftod(struct fpemu *fe, register struct fpn *fp, u_int *res) 372 { 373 register u_int sign = fp->fp_sign << 31; 374 register int exp; 375 376 #define DBL_EXP(e) ((e) << (DBL_FRACBITS & 31)) 377 #define DBL_MASK (DBL_EXP(1) - 1) 378 379 if (ISNAN(fp)) { 380 (void) fpu_shr(fp, FP_NMANT - 1 - DBL_FRACBITS); 381 exp = DBL_EXP_INFNAN; 382 goto done; 383 } 384 if (ISINF(fp)) { 385 sign |= DBL_EXP(DBL_EXP_INFNAN); 386 goto zero; 387 } 388 if (ISZERO(fp)) { 389 zero: res[1] = 0; 390 return (sign); 391 } 392 393 if ((exp = fp->fp_exp + DBL_EXP_BIAS) <= 0) { 394 (void) fpu_shr(fp, FP_NMANT - FP_NG - DBL_FRACBITS - exp); 395 if (fpu_round(fe, fp) && fp->fp_mant[2] == DBL_EXP(1)) { 396 res[1] = 0; 397 return (sign | DBL_EXP(1) | 0); 398 } 399 if ((fe->fe_cx & FSR_NX) || 400 (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT))) 401 fe->fe_cx |= FSR_UF; 402 exp = 0; 403 goto done; 404 } 405 (void) fpu_shr(fp, FP_NMANT - FP_NG - 1 - DBL_FRACBITS); 406 if (fpu_round(fe, fp) && fp->fp_mant[2] == DBL_EXP(2)) 407 exp++; 408 if (exp >= DBL_EXP_INFNAN) { 409 fe->fe_cx |= FSR_OF | FSR_NX; 410 if (toinf(fe, sign)) { 411 res[1] = 0; 412 return (sign | DBL_EXP(DBL_EXP_INFNAN) | 0); 413 } 414 res[1] = ~0; 415 return (sign | DBL_EXP(DBL_EXP_INFNAN) | DBL_MASK); 416 } 417 done: 418 res[1] = fp->fp_mant[3]; 419 return (sign | DBL_EXP(exp) | (fp->fp_mant[2] & DBL_MASK)); 420 } 421 422 /* 423 * fpn -> extended (32 bit high-order result returned; low-order fraction 424 * words left in res[1]..res[3]). Like ftod, which is like ftos ... but 425 * our internal format *is* extended precision, plus 2 bits for guard/round, 426 * so we can avoid a small bit of work. 427 */ 428 u_int 429 fpu_ftoq(struct fpemu *fe, register struct fpn *fp, u_int *res) 430 { 431 register u_int sign = fp->fp_sign << 31; 432 register int exp; 433 434 #define EXT_EXP(e) ((e) << (EXT_FRACBITS & 31)) 435 #define EXT_MASK (EXT_EXP(1) - 1) 436 437 if (ISNAN(fp)) { 438 (void) fpu_shr(fp, 2); /* since we are not rounding */ 439 exp = EXT_EXP_INFNAN; 440 goto done; 441 } 442 if (ISINF(fp)) { 443 sign |= EXT_EXP(EXT_EXP_INFNAN); 444 goto zero; 445 } 446 if (ISZERO(fp)) { 447 zero: res[1] = res[2] = res[3] = 0; 448 return (sign); 449 } 450 451 if ((exp = fp->fp_exp + EXT_EXP_BIAS) <= 0) { 452 (void) fpu_shr(fp, FP_NMANT - FP_NG - EXT_FRACBITS - exp); 453 if (fpu_round(fe, fp) && fp->fp_mant[0] == EXT_EXP(1)) { 454 res[1] = res[2] = res[3] = 0; 455 return (sign | EXT_EXP(1) | 0); 456 } 457 if ((fe->fe_cx & FSR_NX) || 458 (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT))) 459 fe->fe_cx |= FSR_UF; 460 exp = 0; 461 goto done; 462 } 463 /* Since internal == extended, no need to shift here. */ 464 if (fpu_round(fe, fp) && fp->fp_mant[0] == EXT_EXP(2)) 465 exp++; 466 if (exp >= EXT_EXP_INFNAN) { 467 fe->fe_cx |= FSR_OF | FSR_NX; 468 if (toinf(fe, sign)) { 469 res[1] = res[2] = res[3] = 0; 470 return (sign | EXT_EXP(EXT_EXP_INFNAN) | 0); 471 } 472 res[1] = res[2] = res[3] = ~0; 473 return (sign | EXT_EXP(EXT_EXP_INFNAN) | EXT_MASK); 474 } 475 done: 476 res[1] = fp->fp_mant[1]; 477 res[2] = fp->fp_mant[2]; 478 res[3] = fp->fp_mant[3]; 479 return (sign | EXT_EXP(exp) | (fp->fp_mant[0] & EXT_MASK)); 480 } 481 482 /* 483 * Implode an fpn, writing the result into the given space. 484 */ 485 void 486 fpu_implode(struct fpemu *fe, register struct fpn *fp, int type, 487 register u_int *space) 488 { 489 DPRINTF(FPE_INSN, ("fpu_implode: ")); 490 switch (type) { 491 492 case FTYPE_LNG: 493 space[0] = fpu_ftox(fe, fp, space); 494 DPRINTF(FPE_INSN, ("LNG %x %x\n", space[0], space[1])); 495 break; 496 497 case FTYPE_INT: 498 space[0] = fpu_ftoi(fe, fp); 499 DPRINTF(FPE_INSN, ("INT %x\n", space[0])); 500 break; 501 502 case FTYPE_SNG: 503 space[0] = fpu_ftos(fe, fp); 504 DPRINTF(FPE_INSN, ("SNG %x\n", space[0])); 505 break; 506 507 case FTYPE_DBL: 508 space[0] = fpu_ftod(fe, fp, space); 509 DPRINTF(FPE_INSN, ("DBL %x %x\n", space[0], space[1])); 510 break; 511 512 case FTYPE_EXT: 513 /* funky rounding precision options ?? */ 514 space[0] = fpu_ftoq(fe, fp, space); 515 DPRINTF(FPE_INSN, ("EXT %x %x %x %x\n", space[0], space[1], 516 space[2], space[3])); 517 break; 518 519 default: 520 panic("fpu_implode"); 521 } 522 } 523