zig

addf3.zig (6348B) - Raw

---

 const std = @import("std");
 const math = std.math;
 const common = @import("./common.zig");
 const normalize = common.normalize;
 
 /// Ported from:
 ///
 /// https://github.com/llvm/llvm-project/blob/02d85149a05cb1f6dc49f0ba7a2ceca53718ae17/compiler-rt/lib/builtins/fp_add_impl.inc
 pub inline fn addf3(comptime T: type, a: T, b: T) T {
     const bits = @typeInfo(T).float.bits;
     const Z = std.meta.Int(.unsigned, bits);
 
     const typeWidth = bits;
     const significandBits = math.floatMantissaBits(T);
     const fractionalBits = math.floatFractionalBits(T);
     const exponentBits = math.floatExponentBits(T);
 
     const signBit = (@as(Z, 1) << (significandBits + exponentBits));
     const maxExponent = ((1 << exponentBits) - 1);
 
     const integerBit = (@as(Z, 1) << fractionalBits);
     const quietBit = integerBit >> 1;
     const significandMask = (@as(Z, 1) << significandBits) - 1;
 
     const absMask = signBit - 1;
     const qnanRep = @as(Z, @bitCast(math.nan(T))) | quietBit;
 
     var aRep: Z = @bitCast(a);
     var bRep: Z = @bitCast(b);
     const aAbs = aRep & absMask;
     const bAbs = bRep & absMask;
 
     const infRep: Z = @bitCast(math.inf(T));
 
     // Detect if a or b is zero, infinity, or NaN.
     if (aAbs -% @as(Z, 1) >= infRep - @as(Z, 1) or
         bAbs -% @as(Z, 1) >= infRep - @as(Z, 1))
     {
         // NaN + anything = qNaN
         if (aAbs > infRep) return @bitCast(@as(Z, @bitCast(a)) | quietBit);
         // anything + NaN = qNaN
         if (bAbs > infRep) return @bitCast(@as(Z, @bitCast(b)) | quietBit);
 
         if (aAbs == infRep) {
             // +/-infinity + -/+infinity = qNaN
             if ((@as(Z, @bitCast(a)) ^ @as(Z, @bitCast(b))) == signBit) {
                 return @bitCast(qnanRep);
             }
             // +/-infinity + anything remaining = +/- infinity
             else {
                 return a;
             }
         }
 
         // anything remaining + +/-infinity = +/-infinity
         if (bAbs == infRep) return b;
 
         // zero + anything = anything
         if (aAbs == 0) {
             // but we need to get the sign right for zero + zero
             if (bAbs == 0) {
                 return @bitCast(@as(Z, @bitCast(a)) & @as(Z, @bitCast(b)));
             } else {
                 return b;
             }
         }
 
         // anything + zero = anything
         if (bAbs == 0) return a;
     }
 
     // Swap a and b if necessary so that a has the larger absolute value.
     if (bAbs > aAbs) {
         const temp = aRep;
         aRep = bRep;
         bRep = temp;
     }
 
     // Extract the exponent and significand from the (possibly swapped) a and b.
     var aExponent: i32 = @intCast((aRep >> significandBits) & maxExponent);
     var bExponent: i32 = @intCast((bRep >> significandBits) & maxExponent);
     var aSignificand = aRep & significandMask;
     var bSignificand = bRep & significandMask;
 
     // Normalize any denormals, and adjust the exponent accordingly.
     if (aExponent == 0) aExponent = normalize(T, &aSignificand);
     if (bExponent == 0) bExponent = normalize(T, &bSignificand);
 
     // The sign of the result is the sign of the larger operand, a.  If they
     // have opposite signs, we are performing a subtraction; otherwise addition.
     const resultSign = aRep & signBit;
     const subtraction = (aRep ^ bRep) & signBit != 0;
 
     // Shift the significands to give us round, guard and sticky, and or in the
     // implicit significand bit.  (If we fell through from the denormal path it
     // was already set by normalize( ), but setting it twice won't hurt
     // anything.)
     aSignificand = (aSignificand | integerBit) << 3;
     bSignificand = (bSignificand | integerBit) << 3;
 
     // Shift the significand of b by the difference in exponents, with a sticky
     // bottom bit to get rounding correct.
     const @"align": u32 = @intCast(aExponent - bExponent);
     if (@"align" != 0) {
         if (@"align" < typeWidth) {
             const sticky = if (bSignificand << @intCast(typeWidth - @"align") != 0) @as(Z, 1) else 0;
             bSignificand = (bSignificand >> @truncate(@"align")) | sticky;
         } else {
             bSignificand = 1; // sticky; b is known to be non-zero.
         }
     }
     if (subtraction) {
         aSignificand -= bSignificand;
         // If a == -b, return +zero.
         if (aSignificand == 0) return @bitCast(@as(Z, 0));
 
         // If partial cancellation occurred, we need to left-shift the result
         // and adjust the exponent:
         if (aSignificand < integerBit << 3) {
             const shift = @as(i32, @intCast(@clz(aSignificand))) - @as(i32, @intCast(@clz(integerBit << 3)));
             aSignificand <<= @intCast(shift);
             aExponent -= shift;
         }
     } else { // addition
         aSignificand += bSignificand;
 
         // If the addition carried up, we need to right-shift the result and
         // adjust the exponent:
         if (aSignificand & (integerBit << 4) != 0) {
             const sticky = aSignificand & 1;
             aSignificand = aSignificand >> 1 | sticky;
             aExponent += 1;
         }
     }
 
     // If we have overflowed the type, return +/- infinity:
     if (aExponent >= maxExponent) return @bitCast(infRep | resultSign);
 
     if (aExponent <= 0) {
         // Result is denormal; the exponent and round/sticky bits are zero.
         // All we need to do is shift the significand and apply the correct sign.
         aSignificand >>= @intCast(4 - aExponent);
         return @bitCast(resultSign | aSignificand);
     }
 
     // Low three bits are round, guard, and sticky.
     const roundGuardSticky = aSignificand & 0x7;
 
     // Shift the significand into place, and mask off the integer bit, if it's implicit.
     var result = (aSignificand >> 3) & significandMask;
 
     // Insert the exponent and sign.
     result |= @as(Z, @intCast(aExponent)) << significandBits;
     result |= resultSign;
 
     // Final rounding.  The result may overflow to infinity, but that is the
     // correct result in that case.
     if (roundGuardSticky > 0x4) result += 1;
     if (roundGuardSticky == 0x4) result += result & 1;
 
     // Restore any explicit integer bit, if it was rounded off
     if (significandBits != fractionalBits) {
         if ((result >> significandBits) != 0) result |= integerBit;
     }
 
     return @bitCast(result);
 }
 
 test {
     _ = @import("addf3_test.zig");
 }