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