Annotation of nono/fpe/fpu_implode.c, revision 1.1

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 */

unix.superglobalmegacorp.com

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.