|
|
1.1 ! root 1: /* $NetBSD: fpu_implode.c,v 1.15 2013/03/26 11:30:21 isaki Exp $ */ ! 2: ! 3: /* ! 4: * Copyright (c) 1992, 1993 ! 5: * The Regents of the University of California. All rights reserved. ! 6: * ! 7: * This software was developed by the Computer Systems Engineering group ! 8: * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and ! 9: * contributed to Berkeley. ! 10: * ! 11: * All advertising materials mentioning features or use of this software ! 12: * must display the following acknowledgement: ! 13: * This product includes software developed by the University of ! 14: * California, Lawrence Berkeley Laboratory. ! 15: * ! 16: * Redistribution and use in source and binary forms, with or without ! 17: * modification, are permitted provided that the following conditions ! 18: * are met: ! 19: * 1. Redistributions of source code must retain the above copyright ! 20: * notice, this list of conditions and the following disclaimer. ! 21: * 2. Redistributions in binary form must reproduce the above copyright ! 22: * notice, this list of conditions and the following disclaimer in the ! 23: * documentation and/or other materials provided with the distribution. ! 24: * 3. Neither the name of the University nor the names of its contributors ! 25: * may be used to endorse or promote products derived from this software ! 26: * without specific prior written permission. ! 27: * ! 28: * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND ! 29: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE ! 30: * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ! 31: * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE ! 32: * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL ! 33: * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS ! 34: * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) ! 35: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT ! 36: * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY ! 37: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF ! 38: * SUCH DAMAGE. ! 39: * ! 40: * @(#)fpu_implode.c 8.1 (Berkeley) 6/11/93 ! 41: */ ! 42: ! 43: /* ! 44: * FPU subroutines: `implode' internal format numbers into the machine's ! 45: * `packed binary' format. ! 46: */ ! 47: ! 48: #include "fpu_emulate.h" ! 49: ! 50: /* Conversion from internal format -- note asymmetry. */ ! 51: static uint32_t fpu_ftoi(struct fpemu *fe, struct fpn *fp); ! 52: static uint32_t fpu_ftos(struct fpemu *fe, struct fpn *fp); ! 53: static uint32_t fpu_ftod(struct fpemu *fe, struct fpn *fp, uint32_t *); ! 54: static uint32_t fpu_ftox(struct fpemu *fe, struct fpn *fp, uint32_t *); ! 55: static void fpu_round_chkmin(struct fpemu *fe, struct fpn *fp, int, int); ! 56: static void fpu_round_chkinf(struct fpemu *fe, struct fpn *fp, int); ! 57: ! 58: /* ! 59: * Round a number (algorithm from Motorola MC68882 manual, modified for ! 60: * our internal format). Set inexact exception if rounding is required. ! 61: * Return true iff we rounded up. ! 62: * ! 63: * After rounding, we discard the guard and round bits by shifting right ! 64: * 2 bits (a la fpu_shr(), but we do not bother with fp->fp_sticky). ! 65: * This saves effort later. ! 66: * ! 67: * Note that we may leave the value 2.0 in fp->fp_mant; it is the caller's ! 68: * responsibility to fix this if necessary. ! 69: */ ! 70: // FPCR_ROUND (NR, RZ, RP, RM) によって fp (拡張精度) をラウンディングする。 ! 71: // 戻り値はラウンドアップしたら 1(true)。 ! 72: // ラウンディングが発生すると内部 FPSR:EXCP に INEX2 を立てる。 ! 73: // 任意精度のラウンディングは fpu_round_prec() 参照のこと。XXX どうしたもんか ! 74: // ! 75: // 実行後の fp_mant[2] は GR ビットを捨ててあり、つまり fp_mant 全体が当初 ! 76: // に比べて2ビット右シフトされた状態となる。 ! 77: // ラウンドアップによって生じた整数部 2 は何もせずそのまま返されるので、 ! 78: // これを必要に応じてどうにかするのは呼び出し側の責任。 ! 79: // (ただし整数部(小数点位置)も2ビットずれてるはず、でいいのかな?) ! 80: // fp_sticky は常にゼロクリアされる。 ! 81: // fp_exp は変更しない。 ! 82: int ! 83: fpu_round(struct fpemu *fe, struct fpn *fp) ! 84: { ! 85: uint32_t m0, m1, m2; ! 86: int gr, s; ! 87: ! 88: m0 = fp->fp_mant[0]; ! 89: m1 = fp->fp_mant[1]; ! 90: m2 = fp->fp_mant[2]; ! 91: gr = m2 & 3; ! 92: s = fp->fp_sticky; ! 93: ! 94: /* mant >>= FP_NG */ ! 95: m2 = (m2 >> FP_NG) | (m1 << (32 - FP_NG)); ! 96: m1 = (m1 >> FP_NG) | (m0 << (32 - FP_NG)); ! 97: m0 >>= FP_NG; ! 98: ! 99: if ((gr | s) == 0) /* result is exact: no rounding needed */ ! 100: goto rounddown; ! 101: ! 102: fe->fe_fpsr |= FPSR_INEX2; /* inexact */ ! 103: ! 104: /* Go to rounddown to round down; break to round up. */ ! 105: switch (fe->fe_fpcr & FPCR_ROUND) { ! 106: ! 107: case FPCR_NEAR: ! 108: default: ! 109: /* ! 110: * Round only if guard is set (gr & 2). If guard is set, ! 111: * but round & sticky both clear, then we want to round ! 112: * but have a tie, so round to even, i.e., add 1 iff odd. ! 113: */ ! 114: if ((gr & 2) == 0) ! 115: goto rounddown; ! 116: if ((gr & 1) || fp->fp_sticky || (m2 & 1)) ! 117: break; ! 118: goto rounddown; ! 119: ! 120: case FPCR_ZERO: ! 121: /* Round towards zero, i.e., down. */ ! 122: goto rounddown; ! 123: ! 124: case FPCR_MINF: ! 125: /* Round towards -Inf: up if negative, down if positive. */ ! 126: if (fp->fp_sign) ! 127: break; ! 128: goto rounddown; ! 129: ! 130: case FPCR_PINF: ! 131: /* Round towards +Inf: up if positive, down otherwise. */ ! 132: if (!fp->fp_sign) ! 133: break; ! 134: goto rounddown; ! 135: } ! 136: ! 137: /* Bump low bit of mantissa, with carry. */ ! 138: if (++m2 == 0 && ++m1 == 0) ! 139: m0++; ! 140: fp->fp_sticky = 0; ! 141: fp->fp_mant[0] = m0; ! 142: fp->fp_mant[1] = m1; ! 143: fp->fp_mant[2] = m2; ! 144: return (1); ! 145: ! 146: rounddown: ! 147: fp->fp_sticky = 0; ! 148: fp->fp_mant[0] = m0; ! 149: fp->fp_mant[1] = m1; ! 150: fp->fp_mant[2] = m2; ! 151: return (0); ! 152: } ! 153: ! 154: /* ! 155: * For overflow: return true if overflow is to go to +/-Inf, according ! 156: * to the sign of the overflowing result. If false, overflow is to go ! 157: * to the largest magnitude value instead. ! 158: */ ! 159: static int ! 160: toinf(struct fpemu *fe, int sign) ! 161: { ! 162: int inf; ! 163: ! 164: /* look at rounding direction */ ! 165: switch (fe->fe_fpcr & FPCR_ROUND) { ! 166: ! 167: default: ! 168: case FPCR_NEAR: /* the nearest value is always Inf */ ! 169: inf = 1; ! 170: break; ! 171: ! 172: case FPCR_ZERO: /* toward 0 => never towards Inf */ ! 173: inf = 0; ! 174: break; ! 175: ! 176: case FPCR_PINF: /* toward +Inf iff positive */ ! 177: inf = (sign == 0); ! 178: break; ! 179: ! 180: case FPCR_MINF: /* toward -Inf iff negative */ ! 181: inf = sign; ! 182: break; ! 183: } ! 184: return (inf); ! 185: } ! 186: ! 187: /* ! 188: * fpn -> int (int value returned as return value). ! 189: * ! 190: * N.B.: this conversion always rounds towards zero (this is a peculiarity ! 191: * of the SPARC instruction set). ! 192: */ ! 193: static uint32_t ! 194: fpu_ftoi(struct fpemu *fe, struct fpn *fp) ! 195: { ! 196: uint32_t i; ! 197: int sign, exp; ! 198: ! 199: sign = fp->fp_sign; ! 200: switch (fp->fp_class) { ! 201: case FPC_ZERO: ! 202: return (0); ! 203: ! 204: case FPC_NUM: ! 205: /* ! 206: * If exp >= 2^32, overflow. Otherwise shift value right ! 207: * into last mantissa word (this will not exceed 0xffffffff), ! 208: * shifting any guard and round bits out into the sticky ! 209: * bit. Then ``round'' towards zero, i.e., just set an ! 210: * inexact exception if sticky is set (see fpu_round()). ! 211: * If the result is > 0x80000000, or is positive and equals ! 212: * 0x80000000, overflow; otherwise the last fraction word ! 213: * is the result. ! 214: */ ! 215: if ((exp = fp->fp_exp) >= 32) ! 216: break; ! 217: /* NB: the following includes exp < 0 cases */ ! 218: if (fpu_shr(fp, FP_NMANT - 1 - FP_NG - exp) != 0) { ! 219: /* ! 220: * m68881/2 do not underflow when ! 221: * converting to integer ! 222: */ ! 223: ; ! 224: } ! 225: fpu_round(fe, fp); ! 226: i = fp->fp_mant[2]; ! 227: if (i >= ((uint32_t)0x80000000 + sign)) ! 228: break; ! 229: return (sign ? -i : i); ! 230: ! 231: default: /* Inf, qNaN, sNaN */ ! 232: break; ! 233: } ! 234: /* overflow: replace any inexact exception with invalid */ ! 235: fe->fe_fpsr = (fe->fe_fpsr & ~FPSR_INEX2) | FPSR_OPERR; ! 236: return (0x7fffffff + sign); ! 237: } ! 238: ! 239: /* ! 240: * fpn -> single (32 bit single returned as return value). ! 241: * We assume <= 29 bits in a single-precision fraction (1.f part). ! 242: */ ! 243: static uint32_t ! 244: fpu_ftos(struct fpemu *fe, struct fpn *fp) ! 245: { ! 246: uint32_t sign = fp->fp_sign << 31; ! 247: int exp; ! 248: ! 249: #define SNG_EXP(e) ((e) << SNG_FRACBITS) /* makes e an exponent */ ! 250: #define SNG_MASK (SNG_EXP(1) - 1) /* mask for fraction */ ! 251: ! 252: /* Take care of non-numbers first. */ ! 253: if (ISNAN(fp)) { ! 254: /* ! 255: * Preserve upper bits of NaN, per SPARC V8 appendix N. ! 256: * Note that fp->fp_mant[0] has the quiet bit set, ! 257: * even if it is classified as a signalling NaN. ! 258: */ ! 259: (void) fpu_shr(fp, FP_NMANT - 1 - SNG_FRACBITS); ! 260: exp = SNG_EXP_INFNAN; ! 261: goto done; ! 262: } ! 263: if (ISINF(fp)) ! 264: return (sign | SNG_EXP(SNG_EXP_INFNAN)); ! 265: if (ISZERO(fp)) ! 266: return (sign); ! 267: ! 268: /* ! 269: * Normals (including subnormals). Drop all the fraction bits ! 270: * (including the explicit ``implied'' 1 bit) down into the ! 271: * single-precision range. If the number is subnormal, move ! 272: * the ``implied'' 1 into the explicit range as well, and shift ! 273: * right to introduce leading zeroes. Rounding then acts ! 274: * differently for normals and subnormals: the largest subnormal ! 275: * may round to the smallest normal (1.0 x 2^minexp), or may ! 276: * remain subnormal. In the latter case, signal an underflow ! 277: * if the result was inexact or if underflow traps are enabled. ! 278: * ! 279: * Rounding a normal, on the other hand, always produces another ! 280: * normal (although either way the result might be too big for ! 281: * single precision, and cause an overflow). If rounding a ! 282: * normal produces 2.0 in the fraction, we need not adjust that ! 283: * fraction at all, since both 1.0 and 2.0 are zero under the ! 284: * fraction mask. ! 285: * ! 286: * Note that the guard and round bits vanish from the number after ! 287: * rounding. ! 288: */ ! 289: if ((exp = fp->fp_exp + SNG_EXP_BIAS) <= 0) { /* subnormal */ ! 290: fe->fe_fpsr |= FPSR_UNFL; ! 291: /* -NG for g,r; -SNG_FRACBITS-exp for fraction */ ! 292: (void) fpu_shr(fp, FP_NMANT - FP_NG - SNG_FRACBITS - exp); ! 293: if (fpu_round(fe, fp) && fp->fp_mant[2] == SNG_EXP(1)) ! 294: return (sign | SNG_EXP(1) | 0); ! 295: if (fe->fe_fpsr & FPSR_INEX2) { ! 296: /* mc68881/2 don't underflow when converting */ ! 297: fe->fe_fpsr |= FPSR_UNFL; ! 298: } ! 299: return (sign | SNG_EXP(0) | fp->fp_mant[2]); ! 300: } ! 301: /* -FP_NG for g,r; -1 for implied 1; -SNG_FRACBITS for fraction */ ! 302: (void) fpu_shr(fp, FP_NMANT - FP_NG - 1 - SNG_FRACBITS); ! 303: #ifdef DIAGNOSTIC ! 304: if ((fp->fp_mant[2] & SNG_EXP(1 << FP_NG)) == 0) ! 305: panic("fpu_ftos"); ! 306: #endif ! 307: if (fpu_round(fe, fp) && fp->fp_mant[2] == SNG_EXP(2)) ! 308: exp++; ! 309: if (exp >= SNG_EXP_INFNAN) { ! 310: /* overflow to inf or to max single */ ! 311: fe->fe_fpsr |= FPSR_OVFL; ! 312: if (toinf(fe, sign)) { ! 313: fp->fp_class = FPC_INF; ! 314: return (sign | SNG_EXP(SNG_EXP_INFNAN)); ! 315: } ! 316: fe->fe_fpsr |= FPSR_OPERR; ! 317: return (sign | SNG_EXP(SNG_EXP_INFNAN - 1) | SNG_MASK); ! 318: } ! 319: done: ! 320: if ((fp->fp_mant[2] & SNG_MASK) == 0) ! 321: fp->fp_class = FPC_ZERO; ! 322: /* phew, made it */ ! 323: return (sign | SNG_EXP(exp) | (fp->fp_mant[2] & SNG_MASK)); ! 324: } ! 325: ! 326: /* ! 327: * fpn -> double (32 bit high-order result returned; 32-bit low order result ! 328: * left in res[1]). Assumes <= 61 bits in double precision fraction. ! 329: * ! 330: * This code mimics fpu_ftos; see it for comments. ! 331: */ ! 332: static uint32_t ! 333: fpu_ftod(struct fpemu *fe, struct fpn *fp, uint32_t *res) ! 334: { ! 335: uint32_t sign = fp->fp_sign << 31; ! 336: int exp; ! 337: ! 338: #define DBL_EXP(e) ((e) << (DBL_FRACBITS & 31)) ! 339: #define DBL_MASK (DBL_EXP(1) - 1) ! 340: ! 341: if (ISNAN(fp)) { ! 342: (void) fpu_shr(fp, FP_NMANT - 1 - DBL_FRACBITS); ! 343: exp = DBL_EXP_INFNAN; ! 344: goto done; ! 345: } ! 346: if (ISINF(fp)) { ! 347: sign |= DBL_EXP(DBL_EXP_INFNAN); ! 348: res[1] = 0; ! 349: return (sign); ! 350: } ! 351: if (ISZERO(fp)) { ! 352: res[1] = 0; ! 353: return (sign); ! 354: } ! 355: ! 356: if ((exp = fp->fp_exp + DBL_EXP_BIAS) <= 0) { ! 357: fe->fe_fpsr |= FPSR_UNFL; ! 358: (void) fpu_shr(fp, FP_NMANT - FP_NG - DBL_FRACBITS - exp); ! 359: if (fpu_round(fe, fp) && fp->fp_mant[1] == DBL_EXP(1)) { ! 360: res[1] = 0; ! 361: return (sign | DBL_EXP(1) | 0); ! 362: } ! 363: if (fe->fe_fpsr & FPSR_INEX2) { ! 364: /* mc68881/2 don't underflow when converting */ ! 365: fe->fe_fpsr |= FPSR_UNFL; ! 366: } ! 367: exp = 0; ! 368: goto done; ! 369: } ! 370: (void) fpu_shr(fp, FP_NMANT - FP_NG - 1 - DBL_FRACBITS); ! 371: if (fpu_round(fe, fp) && fp->fp_mant[1] == DBL_EXP(2)) ! 372: exp++; ! 373: if (exp >= DBL_EXP_INFNAN) { ! 374: fe->fe_fpsr |= FPSR_OVFL; ! 375: if (toinf(fe, sign)) { ! 376: fp->fp_class = FPC_INF; ! 377: res[1] = 0; ! 378: return (sign | DBL_EXP(DBL_EXP_INFNAN) | 0); ! 379: } ! 380: fe->fe_fpsr |= FPSR_OPERR; ! 381: res[1] = ~0; ! 382: return (sign | DBL_EXP(DBL_EXP_INFNAN) | DBL_MASK); ! 383: } ! 384: done: ! 385: res[1] = fp->fp_mant[2]; ! 386: if (((fp->fp_mant[1] & DBL_MASK) | res[1]) == 0) ! 387: fp->fp_class = FPC_ZERO; ! 388: return (sign | DBL_EXP(exp) | (fp->fp_mant[1] & DBL_MASK)); ! 389: } ! 390: ! 391: /* ! 392: * fpn -> 68k extended (32 bit high-order result returned; two 32-bit low ! 393: * order result left in res[1] & res[2]). Assumes == 64 bits in extended ! 394: * precision fraction. ! 395: * ! 396: * This code mimics fpu_ftos; see it for comments. ! 397: */ ! 398: static uint32_t ! 399: fpu_ftox(struct fpemu *fe, struct fpn *fp, uint32_t *res) ! 400: { ! 401: uint32_t sign = fp->fp_sign << 31; ! 402: int exp; ! 403: ! 404: #define EXT_EXP(e) ((e) << 16) ! 405: /* ! 406: * on m68k extended prec, significand does not share the same long ! 407: * word with exponent ! 408: */ ! 409: #define EXT_MASK 0 ! 410: #define EXT_EXPLICIT1 (1UL << (63 & 31)) ! 411: #define EXT_EXPLICIT2 (1UL << (64 & 31)) ! 412: ! 413: if (ISNAN(fp)) { ! 414: (void) fpu_shr(fp, FP_NMANT - EXT_FRACBITS); ! 415: exp = EXT_EXP_INFNAN; ! 416: goto done; ! 417: } ! 418: if (ISINF(fp)) { ! 419: sign |= EXT_EXP(EXT_EXP_INFNAN); ! 420: res[1] = res[2] = 0; ! 421: return (sign); ! 422: } ! 423: if (ISZERO(fp)) { ! 424: res[1] = res[2] = 0; ! 425: return (sign); ! 426: } ! 427: ! 428: if ((exp = fp->fp_exp + EXT_EXP_BIAS) < 0) { ! 429: fe->fe_fpsr |= FPSR_UNFL; ! 430: /* ! 431: * I'm not sure about this <=... exp==0 doesn't mean ! 432: * it's a denormal in extended format ! 433: */ ! 434: (void) fpu_shr(fp, FP_NMANT - FP_NG - EXT_FRACBITS - exp); ! 435: if (fpu_round(fe, fp) && fp->fp_mant[1] == EXT_EXPLICIT1) { ! 436: res[1] = res[2] = 0; ! 437: return (sign | EXT_EXP(1) | 0); ! 438: } ! 439: if (fe->fe_fpsr & FPSR_INEX2) { ! 440: /* mc68881/2 don't underflow */ ! 441: fe->fe_fpsr |= FPSR_UNFL; ! 442: } ! 443: exp = 0; ! 444: goto done; ! 445: } ! 446: #if (FP_NMANT - FP_NG - EXT_FRACBITS) > 0 ! 447: (void) fpu_shr(fp, FP_NMANT - FP_NG - EXT_FRACBITS); ! 448: #endif ! 449: if (fpu_round(fe, fp) && fp->fp_mant[0] == EXT_EXPLICIT2) { ! 450: exp++; ! 451: fpu_shr(fp, 1); ! 452: } ! 453: if (exp >= EXT_EXP_INFNAN) { ! 454: fe->fe_fpsr |= FPSR_OVFL; ! 455: if (toinf(fe, sign)) { ! 456: fp->fp_class = FPC_INF; ! 457: res[1] = res[2] = 0; ! 458: return (sign | EXT_EXP(EXT_EXP_INFNAN) | 0); ! 459: } ! 460: fe->fe_fpsr |= FPSR_OPERR; ! 461: res[1] = res[2] = ~0; ! 462: return (sign | EXT_EXP(EXT_EXP_INFNAN) | EXT_MASK); ! 463: } ! 464: done: ! 465: res[1] = fp->fp_mant[1]; ! 466: res[2] = fp->fp_mant[2]; ! 467: if ((res[1] | res[2]) == 0) ! 468: fp->fp_class = FPC_ZERO; ! 469: return (sign | EXT_EXP(exp)); ! 470: } ! 471: ! 472: /* ! 473: * Implode an fpn, writing the result into the given space. ! 474: */ ! 475: void ! 476: fpu_implode(struct fpemu *fe, struct fpn *fp, int type, uint32_t *space) ! 477: { ! 478: /* XXX Dont delete exceptions set here: fe->fe_fpsr &= ~FPSR_EXCP; */ ! 479: ! 480: switch (type) { ! 481: case FTYPE_LNG: ! 482: space[0] = fpu_ftoi(fe, fp); ! 483: break; ! 484: ! 485: case FTYPE_SNG: ! 486: space[0] = fpu_ftos(fe, fp); ! 487: break; ! 488: ! 489: case FTYPE_DBL: ! 490: space[0] = fpu_ftod(fe, fp, space); ! 491: break; ! 492: ! 493: case FTYPE_EXT: ! 494: /* funky rounding precision options ?? */ ! 495: space[0] = fpu_ftox(fe, fp, space); ! 496: break; ! 497: ! 498: default: ! 499: /* 何も出来ることがない */ ! 500: break; ! 501: } ! 502: } ! 503: ! 504: #if defined(XM6i_FPE) ! 505: /* ! 506: * Shift the given number left lsh bits. ! 507: * Note that the sticky filed is cleared. ! 508: */ ! 509: static void ! 510: fpu_shl(struct fpn *fp, int lsh) ! 511: { ! 512: uint32_t m0, m1, m2; ! 513: int rsh; ! 514: ! 515: m0 = fp->fp_mant[0]; ! 516: m1 = fp->fp_mant[1]; ! 517: m2 = fp->fp_mant[2]; ! 518: ! 519: while (lsh >= 32) { ! 520: m0 = m1; ! 521: m1 = m2; ! 522: m2 = (fp->fp_sticky != 0) ? 0x80000000 : 0; ! 523: fp->fp_sticky = 0; ! 524: lsh -= 32; ! 525: } ! 526: if (lsh != 0) { ! 527: rsh = 32 - lsh; ! 528: m0 = (m0 << lsh) | (m1 >> rsh); ! 529: m1 = (m1 << lsh) | (m2 >> rsh); ! 530: m2 = (m2 << lsh); ! 531: m2 |= (fp->fp_sticky != 0) ? (1 << (lsh - 1)) : 0; ! 532: fp->fp_sticky = 0; ! 533: } ! 534: fp->fp_mant[0] = m0 & (FP_2 - 1); ! 535: fp->fp_mant[1] = m1; ! 536: fp->fp_mant[2] = m2; ! 537: } ! 538: ! 539: /* ! 540: * Round a number according to FPCR_MODE. ! 541: * (fpu_round() rounds according to FPCR_ROUND as FPCR_MODE = FPCR_EXTD) ! 542: */ ! 543: // fp を FPCR_PREC(精度)/ FPCR_MODE(RN,RZ,RP,RM) によってラウンディングする。 ! 544: // FPSR:EXCP に UNFL、INEX2 を立てる場合がある。 ! 545: // 元からある fpu_round() は拡張精度限定のラウンディング。どうしたもんか。 ! 546: void ! 547: fpu_round_prec(struct fpemu *fe, struct fpn *fp) ! 548: { ! 549: if (fp->fp_class != FPC_NUM) { ! 550: PRINTF("fpu_round_prec class != NUM\n"); ! 551: return; ! 552: } ! 553: ! 554: switch ((fe->fe_fpcr & FPCR_PREC)) { ! 555: case FPCR_SNGL: ! 556: fpu_round_chkmin(fe, fp, SNG_EXP_BIAS, SNG_FRACBITS); ! 557: fpu_round_chkinf(fe, fp, SNG_EXP_BIAS); ! 558: break; ! 559: ! 560: case FPCR_DBL: ! 561: fpu_round_chkmin(fe, fp, DBL_EXP_BIAS, DBL_FRACBITS); ! 562: fpu_round_chkinf(fe, fp, DBL_EXP_BIAS); ! 563: break; ! 564: ! 565: case FPCR_EXTD: ! 566: default: ! 567: // ! 568: // 拡張精度だけ fpu_round_chk{min,inf}() と微妙に処理が異なる…。 ! 569: // ! 570: if (fp->fp_exp < -EXT_EXP_BIAS - EXT_FRACBITS) { ! 571: // 拡張精度にすると指数が小さすぎて表現できない場合 ! 572: fe->fe_fpsr |= FPSR_UNFL; ! 573: DUMPFP("e1:start", fp); ! 574: ! 575: // sticky だけの状態から FPCR_MODE によるラウンディングで ! 576: // fp_mant だけ作る。ここで仮数部ゼロならゼロ。 ! 577: fp->fp_mant[0] = 0; ! 578: fp->fp_mant[1] = 0; ! 579: fp->fp_mant[2] = 0; ! 580: fp->fp_sticky = 1; ! 581: fpu_round(fe, fp); ! 582: if (fp->fp_mant[2] == 0) ! 583: fp->fp_class = FPC_ZERO; ! 584: // fp_mant の LSB に立ってるはずのビットを整数部になるよう正規化。 ! 585: fp->fp_exp = FP_NMANT - EXT_EXP_BIAS - EXT_FRACBITS; ! 586: fpu_norm(fp); ! 587: } else if (fp->fp_exp < -EXT_EXP_BIAS) { ! 588: // 拡張精度にすると非正規化数になる場合 ! 589: fe->fe_fpsr |= FPSR_UNFL; ! 590: DUMPFP("e2:start", fp); ! 591: ! 592: int shift = FP_NMANT - EXT_FRACBITS - 0/*integer bit included*/; ! 593: int effbits = EXT_EXP_BIAS + EXT_FRACBITS + fp->fp_exp; ! 594: PRINTF("effbits=%d\n", effbits); ! 595: shift += EXT_FRACBITS - effbits; ! 596: PRINTF("shift=%d\n", shift); ! 597: fpu_shr(fp, shift - FP_NG); ! 598: DUMPFP("e2:shr ", fp); ! 599: ! 600: int rup = fpu_round(fe, fp); ! 601: uint64_t m = (((uint64_t)fp->fp_mant[1]) << 32) ! 602: | (uint64_t)fp->fp_mant[2]; ! 603: if (rup && m == (1ULL << effbits)) { ! 604: fp->fp_exp++; ! 605: shift--; ! 606: } ! 607: DUMPFP("e2:round", fp); ! 608: ! 609: fpu_shl(fp, shift); ! 610: DUMPFP("e2:shl ", fp); ! 611: ! 612: if (fp->fp_exp <= -EXT_EXP_BIAS - EXT_FRACBITS) { ! 613: PRINTF("e2:zero\n"); ! 614: fp->fp_class = FPC_ZERO; ! 615: } ! 616: } else { ! 617: int shift = FP_NMANT - EXT_FRACBITS - 0/*Integer bit included*/; ! 618: PRINTF("shift=%d\n", shift); ! 619: fpu_shr(fp, shift - FP_NG); ! 620: DUMPFP("e3:shr ", fp); ! 621: ! 622: int rup = fpu_round(fe, fp); ! 623: uint32_t m = fp->fp_mant[1] | fp->fp_mant[2]; ! 624: if (rup && fp->fp_mant[0] == 1 && m == 0) { /* mant == 2.0 */ ! 625: fp->fp_exp++; ! 626: shift--; ! 627: DUMPFP("e3:rup ", fp); ! 628: } else ! 629: DUMPFP("e3:round", fp); ! 630: ! 631: fpu_shl(fp, shift); ! 632: DUMPFP("e3:shl ", fp); ! 633: } ! 634: ! 635: if (fp->fp_exp > EXT_EXP_BIAS) { ! 636: PRINTF("fpu_round_prec EXTD\n"); ! 637: fe->fe_fpsr |= FPSR_OVFL; ! 638: ! 639: // 拡張精度で表現できない大きい値の場合は ! 640: // 正で RN/RP か、負で RN/RM なら Inf、 ! 641: // そうでなければ拡張精度で表現できる最大値。 ! 642: if (fp->fp_sign == 0) { ! 643: switch ((fe->fe_fpcr & FPCR_ROUND)) { ! 644: case FPCR_NEAR: ! 645: case FPCR_PINF: ! 646: fp->fp_class = FPC_INF; ! 647: break; ! 648: } ! 649: } else { ! 650: switch ((fe->fe_fpcr & FPCR_ROUND)) { ! 651: case FPCR_NEAR: ! 652: case FPCR_MINF: ! 653: fp->fp_class = FPC_INF; ! 654: break; ! 655: } ! 656: } ! 657: if (fp->fp_class != FPC_INF) { ! 658: PRINTF("fpu_round_prec make INF\n"); ! 659: // 最大値 ! 660: fp->fp_exp = EXT_EXP_BIAS; ! 661: fp->fp_mant[0] = FP_2 -1; ! 662: fp->fp_mant[1] = 0xffffffff; ! 663: fp->fp_mant[2] = 0xfff80000; ! 664: } ! 665: } ! 666: break; ! 667: } ! 668: } ! 669: ! 670: void ! 671: fpu_round_chkmin(struct fpemu *fe, struct fpn *fp, int EXP_BIAS, int FRACBITS) ! 672: { ! 673: int shift; ! 674: int effbits; ! 675: int rup; ! 676: ! 677: if (fp->fp_exp < -EXP_BIAS - FRACBITS) { ! 678: // 指定精度にすると指数が小さすぎて表現できない場合 ! 679: fe->fe_fpsr |= FPSR_UNFL; ! 680: DUMPFP("s1:start", fp); ! 681: ! 682: // sticky だけの状態から FPCR_MODE によるラウンディングで ! 683: // fp_mant だけ作る。ここで仮数部ゼロならゼロ。 ! 684: fp->fp_mant[0] = 0; ! 685: fp->fp_mant[1] = 0; ! 686: fp->fp_mant[2] = 0; ! 687: fp->fp_sticky = 1; ! 688: fpu_round(fe, fp); ! 689: DUMPFP("s1:round", fp); ! 690: if (fp->fp_mant[2] == 0) ! 691: fp->fp_class = FPC_ZERO; ! 692: // fp_mant の LSB に立ってるはずのビットを整数部になるよう正規化。 ! 693: fp->fp_exp = FP_NMANT - EXP_BIAS - FRACBITS; ! 694: fpu_norm(fp); ! 695: DUMPFP("s1:norm ", fp); ! 696: return; ! 697: } ! 698: if (fp->fp_exp <= -EXP_BIAS) { ! 699: // 指定精度にすると非正規化数になる場合 ! 700: fe->fe_fpsr |= FPSR_UNFL; ! 701: DUMPFP("s2:start", fp); ! 702: ! 703: shift = FP_NMANT - FRACBITS - 0/*No integer bit*/; ! 704: /* ! 705: * 1.XX~XX * 2^-127 => effbits = 22 ! 706: * 1.X * 2^-148 => effbits = 2 ! 707: * 1. * 2^-149 => effbits = 1 ! 708: */ ! 709: effbits = EXP_BIAS + FRACBITS + fp->fp_exp; ! 710: PRINTF("effbits=%d\n", effbits); ! 711: shift += FRACBITS - effbits; ! 712: PRINTF("shift=%d\n", shift); ! 713: fpu_shr(fp, shift - FP_NG); ! 714: DUMPFP("s2:shr ", fp); ! 715: ! 716: rup = fpu_round(fe, fp); ! 717: uint64_t m = (((uint64_t)fp->fp_mant[1]) << 32) ! 718: | (uint64_t)fp->fp_mant[2]; ! 719: if (rup && m == (1ULL << effbits)) { ! 720: fp->fp_exp++; ! 721: shift--; ! 722: } ! 723: DUMPFP("s2:round", fp); ! 724: ! 725: fpu_shl(fp, shift); ! 726: DUMPFP("s2:shl ", fp); ! 727: ! 728: if (fp->fp_exp <= -EXP_BIAS - FRACBITS) { ! 729: PRINTF("s2:zero\n"); ! 730: fp->fp_class = FPC_ZERO; ! 731: } ! 732: return; ! 733: } ! 734: ! 735: // 指定精度の正規化数で表現できる場合 ! 736: DUMPFP("s3:start", fp); ! 737: ! 738: // [0] [1] [2] ! 739: // 8765432109876543210 0 0 0 0 ! 740: // IMMMMMMMMMMMMMMMMMM m..m MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMGR : origin ! 741: // IMMMMMMMMMMMMMMMMMMmmmmmmm : shifted ! 742: // 21098765432109876543210GR ! 743: shift = FP_NMANT - FRACBITS - 1/*Integer bit*/; ! 744: fpu_shr(fp, shift - FP_NG); ! 745: DUMPFP("s3:shr ", fp); ! 746: ! 747: rup = fpu_round(fe, fp); ! 748: uint64_t m = (((uint64_t)fp->fp_mant[1]) << 32) ! 749: | (uint64_t)fp->fp_mant[2]; ! 750: if (rup && (m == (2ULL << FRACBITS))) { ! 751: fp->fp_exp++; ! 752: shift--; ! 753: } ! 754: DUMPFP("s3:round", fp); ! 755: ! 756: fpu_shl(fp, shift); ! 757: DUMPFP("s3:shl ", fp); ! 758: } ! 759: ! 760: void ! 761: fpu_round_chkinf(struct fpemu *fe, struct fpn *fp, int EXP_BIAS) ! 762: { ! 763: if (fp->fp_exp > EXP_BIAS) { ! 764: fe->fe_fpsr |= FPSR_OVFL | FPSR_INEX2; ! 765: ! 766: // 指定精度で表現できない大きい値の場合は ! 767: // 正で RN/RP か、負で RN/RM なら Inf、 ! 768: // そうでなければ指定精度で表現できる最大値。 ! 769: if (fp->fp_sign == 0) { ! 770: switch ((fe->fe_fpcr & FPCR_ROUND)) { ! 771: case FPCR_NEAR: ! 772: case FPCR_PINF: ! 773: fp->fp_class = FPC_INF; ! 774: break; ! 775: } ! 776: } else { ! 777: switch ((fe->fe_fpcr & FPCR_ROUND)) { ! 778: case FPCR_NEAR: ! 779: case FPCR_MINF: ! 780: fp->fp_class = FPC_INF; ! 781: break; ! 782: } ! 783: } ! 784: if (fp->fp_class != FPC_INF) { ! 785: // 最大値 ! 786: fp->fp_exp = EXP_BIAS; ! 787: fp->fp_mant[0] = FP_2 -1; ! 788: switch ((fe->fe_fpcr & FPCR_PREC)) { ! 789: case FPCR_SNGL: ! 790: fp->fp_mant[1] = 0xf8000000; ! 791: fp->fp_mant[2] = 0; ! 792: break; ! 793: case FPCR_DBL: ! 794: fp->fp_mant[1] = 0xffffffff; ! 795: fp->fp_mant[2] = 0xc0000000; ! 796: break; ! 797: case FPCR_EXTD: ! 798: fp->fp_mant[1] = 0xffffffff; ! 799: fp->fp_mant[2] = 0xfff80000; ! 800: break; ! 801: } ! 802: } ! 803: } ! 804: } ! 805: #endif /* XM6i_FPE */
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.