mulsf3.c revision 214152
1214152Sed//===-- lib/mulsf3.c - Single-precision multiplication ------------*- C -*-===// 2214152Sed// 3214152Sed// The LLVM Compiler Infrastructure 4214152Sed// 5214152Sed// This file is distributed under the University of Illinois Open Source 6214152Sed// License. See LICENSE.TXT for details. 7214152Sed// 8214152Sed//===----------------------------------------------------------------------===// 9214152Sed// 10214152Sed// This file implements single-precision soft-float multiplication 11214152Sed// with the IEEE-754 default rounding (to nearest, ties to even). 12214152Sed// 13214152Sed//===----------------------------------------------------------------------===// 14214152Sed 15214152Sed#define SINGLE_PRECISION 16214152Sed#include "fp_lib.h" 17214152Sed 18214152Sedfp_t __mulsf3(fp_t a, fp_t b) { 19214152Sed 20214152Sed const unsigned int aExponent = toRep(a) >> significandBits & maxExponent; 21214152Sed const unsigned int bExponent = toRep(b) >> significandBits & maxExponent; 22214152Sed const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit; 23214152Sed 24214152Sed rep_t aSignificand = toRep(a) & significandMask; 25214152Sed rep_t bSignificand = toRep(b) & significandMask; 26214152Sed int scale = 0; 27214152Sed 28214152Sed // Detect if a or b is zero, denormal, infinity, or NaN. 29214152Sed if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) { 30214152Sed 31214152Sed const rep_t aAbs = toRep(a) & absMask; 32214152Sed const rep_t bAbs = toRep(b) & absMask; 33214152Sed 34214152Sed // NaN * anything = qNaN 35214152Sed if (aAbs > infRep) return fromRep(toRep(a) | quietBit); 36214152Sed // anything * NaN = qNaN 37214152Sed if (bAbs > infRep) return fromRep(toRep(b) | quietBit); 38214152Sed 39214152Sed if (aAbs == infRep) { 40214152Sed // infinity * non-zero = +/- infinity 41214152Sed if (bAbs) return fromRep(aAbs | productSign); 42214152Sed // infinity * zero = NaN 43214152Sed else return fromRep(qnanRep); 44214152Sed } 45214152Sed 46214152Sed if (bAbs == infRep) { 47214152Sed // non-zero * infinity = +/- infinity 48214152Sed if (aAbs) return fromRep(bAbs | productSign); 49214152Sed // zero * infinity = NaN 50214152Sed else return fromRep(qnanRep); 51214152Sed } 52214152Sed 53214152Sed // zero * anything = +/- zero 54214152Sed if (!aAbs) return fromRep(productSign); 55214152Sed // anything * zero = +/- zero 56214152Sed if (!bAbs) return fromRep(productSign); 57214152Sed 58214152Sed // one or both of a or b is denormal, the other (if applicable) is a 59214152Sed // normal number. Renormalize one or both of a and b, and set scale to 60214152Sed // include the necessary exponent adjustment. 61214152Sed if (aAbs < implicitBit) scale += normalize(&aSignificand); 62214152Sed if (bAbs < implicitBit) scale += normalize(&bSignificand); 63214152Sed } 64214152Sed 65214152Sed // Or in the implicit significand bit. (If we fell through from the 66214152Sed // denormal path it was already set by normalize( ), but setting it twice 67214152Sed // won't hurt anything.) 68214152Sed aSignificand |= implicitBit; 69214152Sed bSignificand |= implicitBit; 70214152Sed 71214152Sed // Get the significand of a*b. Before multiplying the significands, shift 72214152Sed // one of them left to left-align it in the field. Thus, the product will 73214152Sed // have (exponentBits + 2) integral digits, all but two of which must be 74214152Sed // zero. Normalizing this result is just a conditional left-shift by one 75214152Sed // and bumping the exponent accordingly. 76214152Sed rep_t productHi, productLo; 77214152Sed wideMultiply(aSignificand, bSignificand << exponentBits, 78214152Sed &productHi, &productLo); 79214152Sed 80214152Sed int productExponent = aExponent + bExponent - exponentBias + scale; 81214152Sed 82214152Sed // Normalize the significand, adjust exponent if needed. 83214152Sed if (productHi & implicitBit) productExponent++; 84214152Sed else wideLeftShift(&productHi, &productLo, 1); 85214152Sed 86214152Sed // If we have overflowed the type, return +/- infinity. 87214152Sed if (productExponent >= maxExponent) return fromRep(infRep | productSign); 88214152Sed 89214152Sed if (productExponent <= 0) { 90214152Sed // Result is denormal before rounding, the exponent is zero and we 91214152Sed // need to shift the significand. 92214152Sed wideRightShiftWithSticky(&productHi, &productLo, 1 - productExponent); 93214152Sed } 94214152Sed 95214152Sed else { 96214152Sed // Result is normal before rounding; insert the exponent. 97214152Sed productHi &= significandMask; 98214152Sed productHi |= (rep_t)productExponent << significandBits; 99214152Sed } 100214152Sed 101214152Sed // Insert the sign of the result: 102214152Sed productHi |= productSign; 103214152Sed 104214152Sed // Final rounding. The final result may overflow to infinity, or underflow 105214152Sed // to zero, but those are the correct results in those cases. 106214152Sed if (productLo > signBit) productHi++; 107214152Sed if (productLo == signBit) productHi += productHi & 1; 108214152Sed return fromRep(productHi); 109214152Sed} 110