Return to fpu_mul.c CVS log

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

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

1.1     ! root        1: /*     $NetBSD: fpu_mul.c,v 1.9 2016/12/06 06:41:14 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_mul.c   8.1 (Berkeley) 6/11/93
        !            41:  */
        !            42: 
        !            43: /*
        !            44:  * Perform an FPU multiply (return x * y).
        !            45:  */
        !            46: 
        !            47: #include "fpu_arith.h"
        !            48: #include "fpu_emulate.h"
        !            49: 
        !            50: /*
        !            51:  * The multiplication algorithm for normal numbers is as follows:
        !            52:  *
        !            53:  * The fraction of the product is built in the usual stepwise fashion.
        !            54:  * Each step consists of shifting the accumulator right one bit
        !            55:  * (maintaining any guard bits) and, if the next bit in y is set,
        !            56:  * adding the multiplicand (x) to the accumulator.  Then, in any case,
        !            57:  * we advance one bit leftward in y.  Algorithmically:
        !            58:  *
        !            59:  *     A = 0;
        !            60:  *     for (bit = 0; bit < FP_NMANT; bit++) {
        !            61:  *             sticky |= A & 1, A >>= 1;
        !            62:  *             if (Y & (1 << bit))
        !            63:  *                     A += X;
        !            64:  *     }
        !            65:  *
        !            66:  * (X and Y here represent the mantissas of x and y respectively.)
        !            67:  * The resultant accumulator (A) is the product's mantissa.  It may
        !            68:  * be as large as 11.11111... in binary and hence may need to be
        !            69:  * shifted right, but at most one bit.
        !            70:  *
        !            71:  * Since we do not have efficient multiword arithmetic, we code the
        !            72:  * accumulator as four separate words, just like any other mantissa.
        !            73:  *
        !            74:  * In the algorithm above, the bits in y are inspected one at a time.
        !            75:  * We will pick them up 32 at a time and then deal with those 32, one
        !            76:  * at a time.  Note, however, that we know several things about y:
        !            77:  *
        !            78:  *    - the guard and round bits at the bottom are sure to be zero;
        !            79:  *
        !            80:  *    - often many low bits are zero (y is often from a single or double
        !            81:  *     precision source);
        !            82:  *
        !            83:  *    - bit FP_NMANT-1 is set, and FP_1*2 fits in a word.
        !            84:  *
        !            85:  * We can also test for 32-zero-bits swiftly.  In this case, the center
        !            86:  * part of the loop---setting sticky, shifting A, and not adding---will
        !            87:  * run 32 times without adding X to A.  We can do a 32-bit shift faster
        !            88:  * by simply moving words.  Since zeros are common, we optimize this case.
        !            89:  * Furthermore, since A is initially zero, we can omit the shift as well
        !            90:  * until we reach a nonzero word.
        !            91:  */
        !            92: struct fpn *
        !            93: fpu_mul(struct fpemu *fe)
        !            94: {
        !            95:        struct fpn *x = &fe->fe_f1, *y = &fe->fe_f2;
        !            96:        uint32_t a2, a1, a0, x2, x1, x0, bit, m;
        !            97:        int sticky;
        !            98:        FPU_DECL_CARRY
        !            99: 
        !           100:        /*
        !           101:         * Put the `heavier' operand on the right (see fpu_emu.h).
        !           102:         * Then we will have one of the following cases, taken in the
        !           103:         * following order:
        !           104:         *
        !           105:         *  - y = NaN.  Implied: if only one is a signalling NaN, y is.
        !           106:         *      The result is y.
        !           107:         *  - y = Inf.  Implied: x != NaN (is 0, number, or Inf: the NaN
        !           108:         *    case was taken care of earlier).
        !           109:         *      If x = 0, the result is NaN.  Otherwise the result
        !           110:         *      is y, with its sign reversed if x is negative.
        !           111:         *  - x = 0.  Implied: y is 0 or number.
        !           112:         *      The result is 0 (with XORed sign as usual).
        !           113:         *  - other.  Implied: both x and y are numbers.
        !           114:         *      The result is x * y (XOR sign, multiply bits, add exponents).
        !           115:         */
        !           116:        ORDER(x, y);
        !           117:        if (ISNAN(y)) {
        !           118:                return (y);
        !           119:        }
        !           120:        if (ISINF(y)) {
        !           121:                if (ISZERO(x))
        !           122:                        return (fpu_newnan(fe));
        !           123:                y->fp_sign ^= x->fp_sign;
        !           124:                return (y);
        !           125:        }
        !           126:        if (ISZERO(x)) {
        !           127:                x->fp_sign ^= y->fp_sign;
        !           128:                return (x);
        !           129:        }
        !           130: 
        !           131:        /*
        !           132:         * Setup.  In the code below, the mask `m' will hold the current
        !           133:         * mantissa byte from y.  The variable `bit' denotes the bit
        !           134:         * within m.  We also define some macros to deal with everything.
        !           135:         */
        !           136:        x2 = x->fp_mant[2];
        !           137:        x1 = x->fp_mant[1];
        !           138:        x0 = x->fp_mant[0];
        !           139:        sticky = a2 = a1 = a0 = 0;
        !           140: 
        !           141: #define        ADD     /* A += X */ \
        !           142:        FPU_ADDS(a2, a2, x2); \
        !           143:        FPU_ADDCS(a1, a1, x1); \
        !           144:        FPU_ADDC(a0, a0, x0)
        !           145: 
        !           146: #define        SHR1    /* A >>= 1, with sticky */ \
        !           147:        sticky |= a2 & 1, \
        !           148:        a2 = (a2 >> 1) | (a1 << 31), a1 = (a1 >> 1) | (a0 << 31), a0 >>= 1
        !           149: 
        !           150: #define        SHR32   /* A >>= 32, with sticky */ \
        !           151:        sticky |= a2, a2 = a1, a1 = a0, a0 = 0
        !           152: 
        !           153: #define        STEP    /* each 1-bit step of the multiplication */ \
        !           154:        SHR1; if (bit & m) { ADD; }; bit <<= 1
        !           155: 
        !           156:        /*
        !           157:         * We are ready to begin.  The multiply loop runs once for each
        !           158:         * of the four 32-bit words.  Some words, however, are special.
        !           159:         * As noted above, the low order bits of Y are often zero.  Even
        !           160:         * if not, the first loop can certainly skip the guard bits.
        !           161:         * The last word of y has its highest 1-bit in position FP_NMANT-1,
        !           162:         * so we stop the loop when we move past that bit.
        !           163:         */
        !           164:        if ((m = y->fp_mant[2]) == 0) {
        !           165:                /* SHR32; */                    /* unneeded since A==0 */
        !           166:        } else {
        !           167:                bit = 1 << FP_NG;
        !           168:                do {
        !           169:                        STEP;
        !           170:                } while (bit != 0);
        !           171:        }
        !           172:        if ((m = y->fp_mant[1]) == 0) {
        !           173:                SHR32;
        !           174:        } else {
        !           175:                bit = 1;
        !           176:                do {
        !           177:                        STEP;
        !           178:                } while (bit != 0);
        !           179:        }
        !           180:        m = y->fp_mant[0];              /* definitely != 0 */
        !           181:        bit = 1;
        !           182:        do {
        !           183:                STEP;
        !           184:        } while (bit <= m);
        !           185: 
        !           186:        /*
        !           187:         * Done with mantissa calculation.  Get exponent and handle
        !           188:         * 11.111...1 case, then put result in place.  We reuse x since
        !           189:         * it already has the right class (FP_NUM).
        !           190:         */
        !           191:        m = x->fp_exp + y->fp_exp;
        !           192:        if (a0 >= FP_2) {
        !           193:                SHR1;
        !           194:                m++;
        !           195:        }
        !           196:        x->fp_sign ^= y->fp_sign;
        !           197:        x->fp_exp = m;
        !           198:        x->fp_sticky = sticky;
        !           199:        x->fp_mant[2] = a2;
        !           200:        x->fp_mant[1] = a1;
        !           201:        x->fp_mant[0] = a0;
        !           202:        return (x);
        !           203: }

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