129 lines · plain
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-exception6//7//===----------------------------------------------------------------------===//8//9// This file implements soft-float multiplication with the IEEE-754 default10// 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 = qNaN33 if (aAbs > infRep)34 return fromRep(toRep(a) | quietBit);35 // anything * NaN = qNaN36 if (bAbs > infRep)37 return fromRep(toRep(b) | quietBit);38 39 if (aAbs == infRep) {40 // infinity * non-zero = +/- infinity41 if (bAbs)42 return fromRep(aAbs | productSign);43 // infinity * zero = NaN44 else45 return fromRep(qnanRep);46 }47 48 if (bAbs == infRep) {49 // non-zero * infinity = +/- infinity50 if (aAbs)51 return fromRep(bAbs | productSign);52 // zero * infinity = NaN53 else54 return fromRep(qnanRep);55 }56 57 // zero * anything = +/- zero58 if (!aAbs)59 return fromRep(productSign);60 // anything * zero = +/- zero61 if (!bAbs)62 return fromRep(productSign);63 64 // One or both of a or b is denormal. The other (if applicable) is a65 // normal number. Renormalize one or both of a and b, and set scale to66 // 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 the74 // denormal path it was already set by normalize( ), but setting it twice75 // won't hurt anything.76 aSignificand |= implicitBit;77 bSignificand |= implicitBit;78 79 // Perform a basic multiplication on the significands. One of them must be80 // 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 else91 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, return101 // zero with the appropriate sign. Mathematically, there is no need to102 // handle this case separately, but we make it a special case to103 // 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 round109 // 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