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