contrib/llvm-project/compiler-rt/lib/builtins/fp_add_impl.inc

   1 //===----- lib/fp_add_impl.inc - floaing point addition -----------*- 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 addition with the IEEE-754 default rounding
  10 // (to nearest, ties to even).
  11 //
  12 //===----------------------------------------------------------------------===//
  13
  14 #include "fp_lib.h"
  15
  16 static __inline fp_t __addXf3__(fp_t a, fp_t b) {
  17   rep_t aRep = toRep(a);
  18   rep_t bRep = toRep(b);
  19   const rep_t aAbs = aRep & absMask;
  20   const rep_t bAbs = bRep & absMask;
  21
  22   // Detect if a or b is zero, infinity, or NaN.
  23   if (aAbs - REP_C(1) >= infRep - REP_C(1) ||
  24       bAbs - REP_C(1) >= infRep - REP_C(1)) {
  25     // NaN + anything = qNaN
  26     if (aAbs > infRep)
  27       return fromRep(toRep(a) | quietBit);
  28     // anything + NaN = qNaN
  29     if (bAbs > infRep)
  30       return fromRep(toRep(b) | quietBit);
  31
  32     if (aAbs == infRep) {
  33       // +/-infinity + -/+infinity = qNaN
  34       if ((toRep(a) ^ toRep(b)) == signBit)
  35         return fromRep(qnanRep);
  36       // +/-infinity + anything remaining = +/- infinity
  37       else
  38         return a;
  39     }
  40
  41     // anything remaining + +/-infinity = +/-infinity
  42     if (bAbs == infRep)
  43       return b;
  44
  45     // zero + anything = anything
  46     if (!aAbs) {
  47       // We need to get the sign right for zero + zero.
  48       if (!bAbs)
  49         return fromRep(toRep(a) & toRep(b));
  50       else
  51         return b;
  52     }
  53
  54     // anything + zero = anything
  55     if (!bAbs)
  56       return a;
  57   }
  58
  59   // Swap a and b if necessary so that a has the larger absolute value.
  60   if (bAbs > aAbs) {
  61     const rep_t temp = aRep;
  62     aRep = bRep;
  63     bRep = temp;
  64   }
  65
  66   // Extract the exponent and significand from the (possibly swapped) a and b.
  67   int aExponent = aRep >> significandBits & maxExponent;
  68   int bExponent = bRep >> significandBits & maxExponent;
  69   rep_t aSignificand = aRep & significandMask;
  70   rep_t bSignificand = bRep & significandMask;
  71
  72   // Normalize any denormals, and adjust the exponent accordingly.
  73   if (aExponent == 0)
  74     aExponent = normalize(&aSignificand);
  75   if (bExponent == 0)
  76     bExponent = normalize(&bSignificand);
  77
  78   // The sign of the result is the sign of the larger operand, a.  If they
  79   // have opposite signs, we are performing a subtraction.  Otherwise, we
  80   // perform addition.
  81   const rep_t resultSign = aRep & signBit;
  82   const bool subtraction = (aRep ^ bRep) & signBit;
  83
  84   // Shift the significands to give us round, guard and sticky, and set the
  85   // implicit significand bit.  If we fell through from the denormal path it
  86   // was already set by normalize( ), but setting it twice won't hurt
  87   // anything.
  88   aSignificand = (aSignificand | implicitBit) << 3;
  89   bSignificand = (bSignificand | implicitBit) << 3;
  90
  91   // Shift the significand of b by the difference in exponents, with a sticky
  92   // bottom bit to get rounding correct.
  93   const unsigned int align = aExponent - bExponent;
  94   if (align) {
  95     if (align < typeWidth) {
  96       const bool sticky = bSignificand << (typeWidth - align);
  97       bSignificand = bSignificand >> align | sticky;
  98     } else {
  99       bSignificand = 1; // Set the sticky bit.  b is known to be non-zero.
 100     }
 101   }
 102   if (subtraction) {
 103     aSignificand -= bSignificand;
 104     // If a == -b, return +zero.
 105     if (aSignificand == 0)
 106       return fromRep(0);
 107
 108     // If partial cancellation occured, we need to left-shift the result
 109     // and adjust the exponent.
 110     if (aSignificand < implicitBit << 3) {
 111       const int shift = rep_clz(aSignificand) - rep_clz(implicitBit << 3);
 112       aSignificand <<= shift;
 113       aExponent -= shift;
 114     }
 115   } else /* addition */ {
 116     aSignificand += bSignificand;
 117
 118     // If the addition carried up, we need to right-shift the result and
 119     // adjust the exponent.
 120     if (aSignificand & implicitBit << 4) {
 121       const bool sticky = aSignificand & 1;
 122       aSignificand = aSignificand >> 1 | sticky;
 123       aExponent += 1;
 124     }
 125   }
 126
 127   // If we have overflowed the type, return +/- infinity.
 128   if (aExponent >= maxExponent)
 129     return fromRep(infRep | resultSign);
 130
 131   if (aExponent <= 0) {
 132     // The result is denormal before rounding.  The exponent is zero and we
 133     // need to shift the significand.
 134     const int shift = 1 - aExponent;
 135     const bool sticky = aSignificand << (typeWidth - shift);
 136     aSignificand = aSignificand >> shift | sticky;
 137     aExponent = 0;
 138   }
 139
 140   // Low three bits are round, guard, and sticky.
 141   const int roundGuardSticky = aSignificand & 0x7;
 142
 143   // Shift the significand into place, and mask off the implicit bit.
 144   rep_t result = aSignificand >> 3 & significandMask;
 145
 146   // Insert the exponent and sign.
 147   result |= (rep_t)aExponent << significandBits;
 148   result |= resultSign;
 149
 150   // Perform the final rounding.  The result may overflow to infinity, but
 151   // that is the correct result in that case.
 152   if (roundGuardSticky > 0x4)
 153     result++;
 154   if (roundGuardSticky == 0x4)
 155     result += result & 1;
 156   return fromRep(result);
 157 }