adddf3.c revision 222656
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#include "abi.h" 16 17#define DOUBLE_PRECISION 18#include "fp_lib.h" 19 20ARM_EABI_FNALIAS(dadd, adddf3); 21 22COMPILER_RT_ABI fp_t 23__adddf3(fp_t a, fp_t b) { 24 25 rep_t aRep = toRep(a); 26 rep_t bRep = toRep(b); 27 const rep_t aAbs = aRep & absMask; 28 const rep_t bAbs = bRep & absMask; 29 30 // Detect if a or b is zero, infinity, or NaN. 31 if (aAbs - 1U >= infRep - 1U || bAbs - 1U >= infRep - 1U) { 32 33 // NaN + anything = qNaN 34 if (aAbs > infRep) return fromRep(toRep(a) | quietBit); 35 // anything + NaN = qNaN 36 if (bAbs > infRep) return fromRep(toRep(b) | quietBit); 37 38 if (aAbs == infRep) { 39 // +/-infinity + -/+infinity = qNaN 40 if ((toRep(a) ^ toRep(b)) == signBit) return fromRep(qnanRep); 41 // +/-infinity + anything remaining = +/- infinity 42 else return a; 43 } 44 45 // anything remaining + +/-infinity = +/-infinity 46 if (bAbs == infRep) return b; 47 48 // zero + anything = anything 49 if (!aAbs) { 50 // but we need to get the sign right for zero + zero 51 if (!bAbs) return fromRep(toRep(a) & toRep(b)); 52 else return b; 53 } 54 55 // anything + zero = anything 56 if (!bAbs) 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) aExponent = normalize(&aSignificand); 74 if (bExponent == 0) bExponent = normalize(&bSignificand); 75 76 // The sign of the result is the sign of the larger operand, a. If they 77 // have opposite signs, we are performing a subtraction; otherwise addition. 78 const rep_t resultSign = aRep & signBit; 79 const bool subtraction = (aRep ^ bRep) & signBit; 80 81 // Shift the significands to give us round, guard and sticky, and or in the 82 // implicit significand bit. (If we fell through from the denormal path it 83 // was already set by normalize( ), but setting it twice won't hurt 84 // anything.) 85 aSignificand = (aSignificand | implicitBit) << 3; 86 bSignificand = (bSignificand | implicitBit) << 3; 87 88 // Shift the significand of b by the difference in exponents, with a sticky 89 // bottom bit to get rounding correct. 90 const int align = aExponent - bExponent; 91 if (align) { 92 if (align < typeWidth) { 93 const bool sticky = bSignificand << (typeWidth - align); 94 bSignificand = bSignificand >> align | sticky; 95 } else { 96 bSignificand = 1; // sticky; b is known to be non-zero. 97 } 98 } 99 100 if (subtraction) { 101 aSignificand -= bSignificand; 102 103 // If a == -b, return +zero. 104 if (aSignificand == 0) return fromRep(0); 105 106 // If partial cancellation occured, we need to left-shift the result 107 // and adjust the exponent: 108 if (aSignificand < implicitBit << 3) { 109 const int shift = rep_clz(aSignificand) - rep_clz(implicitBit << 3); 110 aSignificand <<= shift; 111 aExponent -= shift; 112 } 113 } 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) return fromRep(infRep | resultSign); 129 130 if (aExponent <= 0) { 131 // Result is denormal before rounding; the exponent is zero and we 132 // need to shift the significand. 133 const int shift = 1 - aExponent; 134 const bool sticky = aSignificand << (typeWidth - shift); 135 aSignificand = aSignificand >> shift | sticky; 136 aExponent = 0; 137 } 138 139 // Low three bits are round, guard, and sticky. 140 const int roundGuardSticky = aSignificand & 0x7; 141 142 // Shift the significand into place, and mask off the implicit bit. 143 rep_t result = aSignificand >> 3 & significandMask; 144 145 // Insert the exponent and sign. 146 result |= (rep_t)aExponent << significandBits; 147 result |= resultSign; 148 149 // Final rounding. The result may overflow to infinity, but that is the 150 // correct result in that case. 151 if (roundGuardSticky > 0x4) result++; 152 if (roundGuardSticky == 0x4) result += result & 1; 153 return fromRep(result); 154} 155