Return to fpu_sqrt.c CVS log

Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe

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

1.1     ! root        1: /*     $NetBSD: fpu_sqrt.c,v 1.8 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_sqrt.c  8.1 (Berkeley) 6/11/93
        !            41:  */
        !            42: 
        !            43: /*
        !            44:  * Perform an FPU square root (return sqrt(x)).
        !            45:  */
        !            46: 
        !            47: #include "fpu_arith.h"
        !            48: #include "fpu_emulate.h"
        !            49: 
        !            50: /*
        !            51:  * Our task is to calculate the square root of a floating point number x0.
        !            52:  * This number x normally has the form:
        !            53:  *
        !            54:  *                 exp
        !            55:  *     x = mant * 2            (where 1 <= mant < 2 and exp is an integer)
        !            56:  *
        !            57:  * This can be left as it stands, or the mantissa can be doubled and the
        !            58:  * exponent decremented:
        !            59:  *
        !            60:  *                       exp-1
        !            61:  *     x = (2 * mant) * 2      (where 2 <= 2 * mant < 4)
        !            62:  *
        !            63:  * If the exponent `exp' is even, the square root of the number is best
        !            64:  * handled using the first form, and is by definition equal to:
        !            65:  *
        !            66:  *                             exp/2
        !            67:  *     sqrt(x) = sqrt(mant) * 2
        !            68:  *
        !            69:  * If exp is odd, on the other hand, it is convenient to use the second
        !            70:  * form, giving:
        !            71:  *
        !            72:  *                                 (exp-1)/2
        !            73:  *     sqrt(x) = sqrt(2 * mant) * 2
        !            74:  *
        !            75:  * In the first case, we have
        !            76:  *
        !            77:  *     1 <= mant < 2
        !            78:  *
        !            79:  * and therefore
        !            80:  *
        !            81:  *     sqrt(1) <= sqrt(mant) < sqrt(2)
        !            82:  *
        !            83:  * while in the second case we have
        !            84:  *
        !            85:  *     2 <= 2*mant < 4
        !            86:  *
        !            87:  * and therefore
        !            88:  *
        !            89:  *     sqrt(2) <= sqrt(2*mant) < sqrt(4)
        !            90:  *
        !            91:  * so that in any case, we are sure that
        !            92:  *
        !            93:  *     sqrt(1) <= sqrt(n * mant) < sqrt(4),    n = 1 or 2
        !            94:  *
        !            95:  * or
        !            96:  *
        !            97:  *     1 <= sqrt(n * mant) < 2,                n = 1 or 2.
        !            98:  *
        !            99:  * This root is therefore a properly formed mantissa for a floating
        !           100:  * point number.  The exponent of sqrt(x) is either exp/2 or (exp-1)/2
        !           101:  * as above.  This leaves us with the problem of finding the square root
        !           102:  * of a fixed-point number in the range [1..4).
        !           103:  *
        !           104:  * Though it may not be instantly obvious, the following square root
        !           105:  * algorithm works for any integer x of an even number of bits, provided
        !           106:  * that no overflows occur:
        !           107:  *
        !           108:  *     let q = 0
        !           109:  *     for k = NBITS-1 to 0 step -1 do -- for each digit in the answer...
        !           110:  *             x *= 2                  -- multiply by radix, for next digit
        !           111:  *             if x >= 2q + 2^k then   -- if adding 2^k does not
        !           112:  *                     x -= 2q + 2^k   -- exceed the correct root,
        !           113:  *                     q += 2^k        -- add 2^k and adjust x
        !           114:  *             fi
        !           115:  *     done
        !           116:  *     sqrt = q / 2^(NBITS/2)          -- (and any remainder is in x)
        !           117:  *
        !           118:  * If NBITS is odd (so that k is initially even), we can just add another
        !           119:  * zero bit at the top of x.  Doing so means that q is not going to acquire
        !           120:  * a 1 bit in the first trip around the loop (since x0 < 2^NBITS).  If the
        !           121:  * final value in x is not needed, or can be off by a factor of 2, this is
        !           122:  * equivalant to moving the `x *= 2' step to the bottom of the loop:
        !           123:  *
        !           124:  *     for k = NBITS-1 to 0 step -1 do if ... fi; x *= 2; done
        !           125:  *
        !           126:  * and the result q will then be sqrt(x0) * 2^floor(NBITS / 2).
        !           127:  * (Since the algorithm is destructive on x, we will call x's initial
        !           128:  * value, for which q is some power of two times its square root, x0.)
        !           129:  *
        !           130:  * If we insert a loop invariant y = 2q, we can then rewrite this using
        !           131:  * C notation as:
        !           132:  *
        !           133:  *     q = y = 0; x = x0;
        !           134:  *     for (k = NBITS; --k >= 0;) {
        !           135:  * #if (NBITS is even)
        !           136:  *             x *= 2;
        !           137:  * #endif
        !           138:  *             t = y + (1 << k);
        !           139:  *             if (x >= t) {
        !           140:  *                     x -= t;
        !           141:  *                     q += 1 << k;
        !           142:  *                     y += 1 << (k + 1);
        !           143:  *             }
        !           144:  * #if (NBITS is odd)
        !           145:  *             x *= 2;
        !           146:  * #endif
        !           147:  *     }
        !           148:  *
        !           149:  * If x0 is fixed point, rather than an integer, we can simply alter the
        !           150:  * scale factor between q and sqrt(x0).  As it happens, we can easily arrange
        !           151:  * for the scale factor to be 2**0 or 1, so that sqrt(x0) == q.
        !           152:  *
        !           153:  * In our case, however, x0 (and therefore x, y, q, and t) are multiword
        !           154:  * integers, which adds some complication.  But note that q is built one
        !           155:  * bit at a time, from the top down, and is not used itself in the loop
        !           156:  * (we use 2q as held in y instead).  This means we can build our answer
        !           157:  * in an integer, one word at a time, which saves a bit of work.  Also,
        !           158:  * since 1 << k is always a `new' bit in q, 1 << k and 1 << (k+1) are
        !           159:  * `new' bits in y and we can set them with an `or' operation rather than
        !           160:  * a full-blown multiword add.
        !           161:  *
        !           162:  * We are almost done, except for one snag.  We must prove that none of our
        !           163:  * intermediate calculations can overflow.  We know that x0 is in [1..4)
        !           164:  * and therefore the square root in q will be in [1..2), but what about x,
        !           165:  * y, and t?
        !           166:  *
        !           167:  * We know that y = 2q at the beginning of each loop.  (The relation only
        !           168:  * fails temporarily while y and q are being updated.)  Since q < 2, y < 4.
        !           169:  * The sum in t can, in our case, be as much as y+(1<<1) = y+2 < 6, and.
        !           170:  * Furthermore, we can prove with a bit of work that x never exceeds y by
        !           171:  * more than 2, so that even after doubling, 0 <= x < 8.  (This is left as
        !           172:  * an exercise to the reader, mostly because I have become tired of working
        !           173:  * on this comment.)
        !           174:  *
        !           175:  * If our floating point mantissas (which are of the form 1.frac) occupy
        !           176:  * B+1 bits, our largest intermediary needs at most B+3 bits, or two extra.
        !           177:  * In fact, we want even one more bit (for a carry, to avoid compares), or
        !           178:  * three extra.  There is a comment in fpu_emu.h reminding maintainers of
        !           179:  * this, so we have some justification in assuming it.
        !           180:  */
        !           181: struct fpn *
        !           182: fpu_sqrt(struct fpemu *fe)
        !           183: {
        !           184:        struct fpn *x = &fe->fe_f2;
        !           185:        uint32_t bit, q, tt;
        !           186:        uint32_t x0, x1, x2;
        !           187:        uint32_t y0, y1, y2;
        !           188:        uint32_t d0, d1, d2;
        !           189:        int e;
        !           190:        FPU_DECL_CARRY
        !           191: 
        !           192:        /*
        !           193:         * Take care of special cases first.  In order:
        !           194:         *
        !           195:         *      sqrt(NaN) = NaN
        !           196:         *      sqrt(+0) = +0
        !           197:         *      sqrt(-0) = -0
        !           198:         *      sqrt(x < 0) = NaN       (including sqrt(-Inf))
        !           199:         *      sqrt(+Inf) = +Inf
        !           200:         *
        !           201:         * Then all that remains are numbers with mantissas in [1..2).
        !           202:         */
        !           203:        if (ISNAN(x) || ISZERO(x))
        !           204:                return (x);
        !           205:        if (x->fp_sign)
        !           206:                return (fpu_newnan(fe));
        !           207:        if (ISINF(x))
        !           208:                return (x);
        !           209: 
        !           210:        /*
        !           211:         * Calculate result exponent.  As noted above, this may involve
        !           212:         * doubling the mantissa.  We will also need to double x each
        !           213:         * time around the loop, so we define a macro for this here, and
        !           214:         * we break out the multiword mantissa.
        !           215:         */
        !           216: #ifdef FPU_SHL1_BY_ADD
        !           217: #define        DOUBLE_X { \
        !           218:        FPU_ADDS(x2, x2, x2); \
        !           219:        FPU_ADDCS(x1, x1, x1); FPU_ADDC(x0, x0, x0); \
        !           220: }
        !           221: #else
        !           222: #define        DOUBLE_X { \
        !           223:        x0 = (x0 << 1) | (x1 >> 31); x1 = (x1 << 1) | (x2 >> 31); \
        !           224:        x2 <<= 1; \
        !           225: }
        !           226: #endif
        !           227: #if (FP_NMANT & 1) != 0
        !           228: # define ODD_DOUBLE    DOUBLE_X
        !           229: # define EVEN_DOUBLE   /* nothing */
        !           230: #else
        !           231: # define ODD_DOUBLE    /* nothing */
        !           232: # define EVEN_DOUBLE   DOUBLE_X
        !           233: #endif
        !           234:        x0 = x->fp_mant[0];
        !           235:        x1 = x->fp_mant[1];
        !           236:        x2 = x->fp_mant[2];
        !           237:        e = x->fp_exp;
        !           238:        if (e & 1)              /* exponent is odd; use sqrt(2mant) */
        !           239:                DOUBLE_X;
        !           240:        /* THE FOLLOWING ASSUMES THAT RIGHT SHIFT DOES SIGN EXTENSION */
        !           241:        x->fp_exp = e >> 1;     /* calculates (e&1 ? (e-1)/2 : e/2 */
        !           242: 
        !           243:        /*
        !           244:         * Now calculate the mantissa root.  Since x is now in [1..4),
        !           245:         * we know that the first trip around the loop will definitely
        !           246:         * set the top bit in q, so we can do that manually and start
        !           247:         * the loop at the next bit down instead.  We must be sure to
        !           248:         * double x correctly while doing the `known q=1.0'.
        !           249:         *
        !           250:         * We do this one mantissa-word at a time, as noted above, to
        !           251:         * save work.  To avoid `(1 << 31) << 1', we also do the top bit
        !           252:         * outside of each per-word loop.
        !           253:         *
        !           254:         * The calculation `t = y + bit' breaks down into `t0 = y0, ...,
        !           255:         * t2 = y2, t? |= bit' for the appropriate word.  Since the bit
        !           256:         * is always a `new' one, this means that three of the `t?'s are
        !           257:         * just the corresponding `y?'; we use `#define's here for this.
        !           258:         * The variable `tt' holds the actual `t?' variable.
        !           259:         */
        !           260: 
        !           261:        /* calculate q0 */
        !           262: #define        t0 tt
        !           263:        bit = FP_1;
        !           264:        EVEN_DOUBLE;
        !           265:        /* if (x >= (t0 = y0 | bit)) { */       /* always true */
        !           266:                q = bit;
        !           267:                x0 -= bit;
        !           268:                y0 = bit << 1;
        !           269:        /* } */
        !           270:        ODD_DOUBLE;
        !           271:        while ((bit >>= 1) != 0) {      /* for remaining bits in q0 */
        !           272:                EVEN_DOUBLE;
        !           273:                t0 = y0 | bit;          /* t = y + bit */
        !           274:                if (x0 >= t0) {         /* if x >= t then */
        !           275:                        x0 -= t0;       /*      x -= t */
        !           276:                        q |= bit;       /*      q += bit */
        !           277:                        y0 |= bit << 1; /*      y += bit << 1 */
        !           278:                }
        !           279:                ODD_DOUBLE;
        !           280:        }
        !           281:        x->fp_mant[0] = q;
        !           282: #undef t0
        !           283: 
        !           284:        /* calculate q1.  note (y0&1)==0. */
        !           285: #define t0 y0
        !           286: #define t1 tt
        !           287:        q = 0;
        !           288:        y1 = 0;
        !           289:        bit = 1 << 31;
        !           290:        EVEN_DOUBLE;
        !           291:        t1 = bit;
        !           292:        FPU_SUBS(d1, x1, t1);
        !           293:        FPU_SUBC(d0, x0, t0);           /* d = x - t */
        !           294:        if ((int)d0 >= 0) {             /* if d >= 0 (i.e., x >= t) then */
        !           295:                x0 = d0, x1 = d1;       /*      x -= t */
        !           296:                q = bit;                /*      q += bit */
        !           297:                y0 |= 1;                /*      y += bit << 1 */
        !           298:        }
        !           299:        ODD_DOUBLE;
        !           300:        while ((bit >>= 1) != 0) {      /* for remaining bits in q1 */
        !           301:                EVEN_DOUBLE;            /* as before */
        !           302:                t1 = y1 | bit;
        !           303:                FPU_SUBS(d1, x1, t1);
        !           304:                FPU_SUBC(d0, x0, t0);
        !           305:                if ((int)d0 >= 0) {
        !           306:                        x0 = d0, x1 = d1;
        !           307:                        q |= bit;
        !           308:                        y1 |= bit << 1;
        !           309:                }
        !           310:                ODD_DOUBLE;
        !           311:        }
        !           312:        x->fp_mant[1] = q;
        !           313: #undef t1
        !           314: 
        !           315:        /* calculate q2.  note (y1&1)==0; y0 (aka t0) is fixed. */
        !           316: #define t1 y1
        !           317: #define t2 tt
        !           318:        q = 0;
        !           319:        y2 = 0;
        !           320:        bit = 1 << 31;
        !           321:        EVEN_DOUBLE;
        !           322:        t2 = bit;
        !           323:        FPU_SUBS(d2, x2, t2);
        !           324:        FPU_SUBCS(d1, x1, t1);
        !           325:        FPU_SUBC(d0, x0, t0);
        !           326:        if ((int)d0 >= 0) {
        !           327:                x0 = d0, x1 = d1, x2 = d2;
        !           328:                q |= bit;
        !           329:                y1 |= 1;                /* now t1, y1 are set in concrete */
        !           330:        }
        !           331:        ODD_DOUBLE;
        !           332:        while ((bit >>= 1) != 0) {
        !           333:                EVEN_DOUBLE;
        !           334:                t2 = y2 | bit;
        !           335:                FPU_SUBS(d2, x2, t2);
        !           336:                FPU_SUBCS(d1, x1, t1);
        !           337:                FPU_SUBC(d0, x0, t0);
        !           338:                if ((int)d0 >= 0) {
        !           339:                        x0 = d0, x1 = d1, x2 = d2;
        !           340:                        q |= bit;
        !           341:                        y2 |= bit << 1;
        !           342:                }
        !           343:                ODD_DOUBLE;
        !           344:        }
        !           345:        x->fp_mant[2] = q;
        !           346: #undef t2
        !           347: 
        !           348:        /*
        !           349:         * The result, which includes guard and round bits, is exact iff
        !           350:         * x is now zero; any nonzero bits in x represent sticky bits.
        !           351:         */
        !           352:        x->fp_sticky = x0 | x1 | x2;
        !           353:        return (x);
        !           354: }

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