1//===---- lib/fp_mul_impl.inc - floating point multiplication -----*- C -*-===// 2// 3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. 4// See https://llvm.org/LICENSE.txt for license information. 5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception 6// 7//===----------------------------------------------------------------------===// 8// 9// This file implements soft-float multiplication with the IEEE-754 default 10// rounding (to nearest, ties to even). 11// 12//===----------------------------------------------------------------------===// 13 14#include "fp_lib.h" 15 16static __inline fp_t __mulXf3__(fp_t a, fp_t b) { 17 const unsigned int aExponent = toRep(a) >> significandBits & maxExponent; 18 const unsigned int bExponent = toRep(b) >> significandBits & maxExponent; 19 const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit; 20 21 rep_t aSignificand = toRep(a) & significandMask; 22 rep_t bSignificand = toRep(b) & significandMask; 23 int scale = 0; 24 25 // Detect if a or b is zero, denormal, infinity, or NaN. 26 if (aExponent - 1U >= maxExponent - 1U || 27 bExponent - 1U >= maxExponent - 1U) { 28 29 const rep_t aAbs = toRep(a) & absMask; 30 const rep_t bAbs = toRep(b) & absMask; 31 32 // NaN * anything = qNaN 33 if (aAbs > infRep) 34 return fromRep(toRep(a) | quietBit); 35 // anything * NaN = qNaN 36 if (bAbs > infRep) 37 return fromRep(toRep(b) | quietBit); 38 39 if (aAbs == infRep) { 40 // infinity * non-zero = +/- infinity 41 if (bAbs) 42 return fromRep(aAbs | productSign); 43 // infinity * zero = NaN 44 else 45 return fromRep(qnanRep); 46 } 47 48 if (bAbs == infRep) { 49 // non-zero * infinity = +/- infinity 50 if (aAbs) 51 return fromRep(bAbs | productSign); 52 // zero * infinity = NaN 53 else 54 return fromRep(qnanRep); 55 } 56 57 // zero * anything = +/- zero 58 if (!aAbs) 59 return fromRep(productSign); 60 // anything * zero = +/- zero 61 if (!bAbs) 62 return fromRep(productSign); 63 64 // One or both of a or b is denormal. The other (if applicable) is a 65 // normal number. Renormalize one or both of a and b, and set scale to 66 // include the necessary exponent adjustment. 67 if (aAbs < implicitBit) 68 scale += normalize(&aSignificand); 69 if (bAbs < implicitBit) 70 scale += normalize(&bSignificand); 71 } 72 73 // Set the implicit significand bit. If we fell through from the 74 // denormal path it was already set by normalize( ), but setting it twice 75 // won't hurt anything. 76 aSignificand |= implicitBit; 77 bSignificand |= implicitBit; 78 79 // Perform a basic multiplication on the significands. One of them must be 80 // shifted beforehand to be aligned with the exponent. 81 rep_t productHi, productLo; 82 wideMultiply(aSignificand, bSignificand << exponentBits, &productHi, 83 &productLo); 84 85 int productExponent = aExponent + bExponent - exponentBias + scale; 86 87 // Normalize the significand and adjust the exponent if needed. 88 if (productHi & implicitBit) 89 productExponent++; 90 else 91 wideLeftShift(&productHi, &productLo, 1); 92 93 // If we have overflowed the type, return +/- infinity. 94 if (productExponent >= maxExponent) 95 return fromRep(infRep | productSign); 96 97 if (productExponent <= 0) { 98 // The result is denormal before rounding. 99 // 100 // If the result is so small that it just underflows to zero, return 101 // zero with the appropriate sign. Mathematically, there is no need to 102 // handle this case separately, but we make it a special case to 103 // simplify the shift logic. 104 const unsigned int shift = REP_C(1) - (unsigned int)productExponent; 105 if (shift >= typeWidth) 106 return fromRep(productSign); 107 108 // Otherwise, shift the significand of the result so that the round 109 // bit is the high bit of productLo. 110 wideRightShiftWithSticky(&productHi, &productLo, shift); 111 } else { 112 // The result is normal before rounding. Insert the exponent. 113 productHi &= significandMask; 114 productHi |= (rep_t)productExponent << significandBits; 115 } 116 117 // Insert the sign of the result. 118 productHi |= productSign; 119 120 // Perform the final rounding. The final result may overflow to infinity, 121 // or underflow to zero, but those are the correct results in those cases. 122 // We use the default IEEE-754 round-to-nearest, ties-to-even rounding mode. 123 if (productLo > signBit) 124 productHi++; 125 if (productLo == signBit) 126 productHi += productHi & 1; 127 return fromRep(productHi); 128} 129