adddf3.c revision 229135
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 18ARM_EABI_FNALIAS(dadd, adddf3); 19 20COMPILER_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 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