15c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon//===-- lib/mulsf3.c - Single-precision multiplication ------------*- C -*-===// 25c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// 35c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// The LLVM Compiler Infrastructure 45c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// 59ad441ffec97db647fee3725b3424284fb913e14Howard Hinnant// This file is dual licensed under the MIT and the University of Illinois Open 69ad441ffec97db647fee3725b3424284fb913e14Howard Hinnant// Source Licenses. See LICENSE.TXT for details. 75c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// 85c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon//===----------------------------------------------------------------------===// 95c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// 105c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// This file implements single-precision soft-float multiplication 115c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// with the IEEE-754 default rounding (to nearest, ties to even). 125c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon// 135c6d2ecb9c43d8b836b3203a243e24703d473765Stephen Canon//===----------------------------------------------------------------------===// 14e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 15e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon#define SINGLE_PRECISION 16e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon#include "fp_lib.h" 17e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 180193b74976719b8aea4cb8874ba36b75836a8d6eChandler CarruthARM_EABI_FNALIAS(fmul, mulsf3) 1937b97d1cf4501b94347e0b4e880f4b25825a289fAnton Korobeynikov 201c5f89b1dd741135a4007ab577723d422f421eecAnton KorobeynikovCOMPILER_RT_ABI fp_t 211c5f89b1dd741135a4007ab577723d422f421eecAnton Korobeynikov__mulsf3(fp_t a, fp_t b) { 22e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 23e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon const unsigned int aExponent = toRep(a) >> significandBits & maxExponent; 24e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon const unsigned int bExponent = toRep(b) >> significandBits & maxExponent; 25e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit; 26e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 27e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon rep_t aSignificand = toRep(a) & significandMask; 28e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon rep_t bSignificand = toRep(b) & significandMask; 29e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon int scale = 0; 30e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 31e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Detect if a or b is zero, denormal, infinity, or NaN. 32e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) { 33e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 34e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon const rep_t aAbs = toRep(a) & absMask; 35e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon const rep_t bAbs = toRep(b) & absMask; 36e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 37e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // NaN * anything = qNaN 38e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (aAbs > infRep) return fromRep(toRep(a) | quietBit); 39e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // anything * NaN = qNaN 40e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (bAbs > infRep) return fromRep(toRep(b) | quietBit); 41e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 42e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (aAbs == infRep) { 43e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // infinity * non-zero = +/- infinity 44e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (bAbs) return fromRep(aAbs | productSign); 45e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // infinity * zero = NaN 46e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon else return fromRep(qnanRep); 47e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon } 48e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 49e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (bAbs == infRep) { 50e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // non-zero * infinity = +/- infinity 51e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (aAbs) return fromRep(bAbs | productSign); 52e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // zero * infinity = NaN 53e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon else return fromRep(qnanRep); 54e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon } 55e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 56e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // zero * anything = +/- zero 57e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (!aAbs) return fromRep(productSign); 58e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // anything * zero = +/- zero 59e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (!bAbs) return fromRep(productSign); 60e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 61e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // one or both of a or b is denormal, the other (if applicable) is a 62e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // normal number. Renormalize one or both of a and b, and set scale to 63e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // include the necessary exponent adjustment. 64e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (aAbs < implicitBit) scale += normalize(&aSignificand); 65e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (bAbs < implicitBit) scale += normalize(&bSignificand); 66e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon } 67e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 68e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Or in the implicit significand bit. (If we fell through from the 69e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // denormal path it was already set by normalize( ), but setting it twice 70e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // won't hurt anything.) 71e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon aSignificand |= implicitBit; 72e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon bSignificand |= implicitBit; 73e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 74e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Get the significand of a*b. Before multiplying the significands, shift 75e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // one of them left to left-align it in the field. Thus, the product will 76e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // have (exponentBits + 2) integral digits, all but two of which must be 77e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // zero. Normalizing this result is just a conditional left-shift by one 78e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // and bumping the exponent accordingly. 79e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon rep_t productHi, productLo; 80e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon wideMultiply(aSignificand, bSignificand << exponentBits, 81e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon &productHi, &productLo); 82e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 83e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon int productExponent = aExponent + bExponent - exponentBias + scale; 84e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 85e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Normalize the significand, adjust exponent if needed. 86e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (productHi & implicitBit) productExponent++; 87e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon else wideLeftShift(&productHi, &productLo, 1); 88e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 89e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // If we have overflowed the type, return +/- infinity. 90e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (productExponent >= maxExponent) return fromRep(infRep | productSign); 91e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 92e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (productExponent <= 0) { 93e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Result is denormal before rounding, the exponent is zero and we 94e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // need to shift the significand. 950499b8433538a1913ca955566cd34368752bf7afJoerg Sonnenberger wideRightShiftWithSticky(&productHi, &productLo, 1U - (unsigned)productExponent); 96e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon } 97e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 98e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon else { 99e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Result is normal before rounding; insert the exponent. 100e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon productHi &= significandMask; 101e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon productHi |= (rep_t)productExponent << significandBits; 102e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon } 103e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 104e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Insert the sign of the result: 105e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon productHi |= productSign; 106e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon 107e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // Final rounding. The final result may overflow to infinity, or underflow 108e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon // to zero, but those are the correct results in those cases. 109e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (productLo > signBit) productHi++; 110e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon if (productLo == signBit) productHi += productHi & 1; 111e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon return fromRep(productHi); 112e5086322295e5a345af02d09abfcf8ddca2d0897Stephen Canon} 113