master
1const std = @import("std");
2const math = std.math;
3const common = @import("./common.zig");
4const normalize = common.normalize;
5
6/// Ported from:
7///
8/// https://github.com/llvm/llvm-project/blob/02d85149a05cb1f6dc49f0ba7a2ceca53718ae17/compiler-rt/lib/builtins/fp_add_impl.inc
9pub inline fn addf3(comptime T: type, a: T, b: T) T {
10 const bits = @typeInfo(T).float.bits;
11 const Z = std.meta.Int(.unsigned, bits);
12
13 const typeWidth = bits;
14 const significandBits = math.floatMantissaBits(T);
15 const fractionalBits = math.floatFractionalBits(T);
16 const exponentBits = math.floatExponentBits(T);
17
18 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
19 const maxExponent = ((1 << exponentBits) - 1);
20
21 const integerBit = (@as(Z, 1) << fractionalBits);
22 const quietBit = integerBit >> 1;
23 const significandMask = (@as(Z, 1) << significandBits) - 1;
24
25 const absMask = signBit - 1;
26 const qnanRep = @as(Z, @bitCast(math.nan(T))) | quietBit;
27
28 var aRep: Z = @bitCast(a);
29 var bRep: Z = @bitCast(b);
30 const aAbs = aRep & absMask;
31 const bAbs = bRep & absMask;
32
33 const infRep: Z = @bitCast(math.inf(T));
34
35 // Detect if a or b is zero, infinity, or NaN.
36 if (aAbs -% @as(Z, 1) >= infRep - @as(Z, 1) or
37 bAbs -% @as(Z, 1) >= infRep - @as(Z, 1))
38 {
39 // NaN + anything = qNaN
40 if (aAbs > infRep) return @bitCast(@as(Z, @bitCast(a)) | quietBit);
41 // anything + NaN = qNaN
42 if (bAbs > infRep) return @bitCast(@as(Z, @bitCast(b)) | quietBit);
43
44 if (aAbs == infRep) {
45 // +/-infinity + -/+infinity = qNaN
46 if ((@as(Z, @bitCast(a)) ^ @as(Z, @bitCast(b))) == signBit) {
47 return @bitCast(qnanRep);
48 }
49 // +/-infinity + anything remaining = +/- infinity
50 else {
51 return a;
52 }
53 }
54
55 // anything remaining + +/-infinity = +/-infinity
56 if (bAbs == infRep) return b;
57
58 // zero + anything = anything
59 if (aAbs == 0) {
60 // but we need to get the sign right for zero + zero
61 if (bAbs == 0) {
62 return @bitCast(@as(Z, @bitCast(a)) & @as(Z, @bitCast(b)));
63 } else {
64 return b;
65 }
66 }
67
68 // anything + zero = anything
69 if (bAbs == 0) return a;
70 }
71
72 // Swap a and b if necessary so that a has the larger absolute value.
73 if (bAbs > aAbs) {
74 const temp = aRep;
75 aRep = bRep;
76 bRep = temp;
77 }
78
79 // Extract the exponent and significand from the (possibly swapped) a and b.
80 var aExponent: i32 = @intCast((aRep >> significandBits) & maxExponent);
81 var bExponent: i32 = @intCast((bRep >> significandBits) & maxExponent);
82 var aSignificand = aRep & significandMask;
83 var bSignificand = bRep & significandMask;
84
85 // Normalize any denormals, and adjust the exponent accordingly.
86 if (aExponent == 0) aExponent = normalize(T, &aSignificand);
87 if (bExponent == 0) bExponent = normalize(T, &bSignificand);
88
89 // The sign of the result is the sign of the larger operand, a. If they
90 // have opposite signs, we are performing a subtraction; otherwise addition.
91 const resultSign = aRep & signBit;
92 const subtraction = (aRep ^ bRep) & signBit != 0;
93
94 // Shift the significands to give us round, guard and sticky, and or in the
95 // implicit significand bit. (If we fell through from the denormal path it
96 // was already set by normalize( ), but setting it twice won't hurt
97 // anything.)
98 aSignificand = (aSignificand | integerBit) << 3;
99 bSignificand = (bSignificand | integerBit) << 3;
100
101 // Shift the significand of b by the difference in exponents, with a sticky
102 // bottom bit to get rounding correct.
103 const @"align": u32 = @intCast(aExponent - bExponent);
104 if (@"align" != 0) {
105 if (@"align" < typeWidth) {
106 const sticky = if (bSignificand << @intCast(typeWidth - @"align") != 0) @as(Z, 1) else 0;
107 bSignificand = (bSignificand >> @truncate(@"align")) | sticky;
108 } else {
109 bSignificand = 1; // sticky; b is known to be non-zero.
110 }
111 }
112 if (subtraction) {
113 aSignificand -= bSignificand;
114 // If a == -b, return +zero.
115 if (aSignificand == 0) return @bitCast(@as(Z, 0));
116
117 // If partial cancellation occurred, we need to left-shift the result
118 // and adjust the exponent:
119 if (aSignificand < integerBit << 3) {
120 const shift = @as(i32, @intCast(@clz(aSignificand))) - @as(i32, @intCast(@clz(integerBit << 3)));
121 aSignificand <<= @intCast(shift);
122 aExponent -= shift;
123 }
124 } else { // addition
125 aSignificand += bSignificand;
126
127 // If the addition carried up, we need to right-shift the result and
128 // adjust the exponent:
129 if (aSignificand & (integerBit << 4) != 0) {
130 const sticky = aSignificand & 1;
131 aSignificand = aSignificand >> 1 | sticky;
132 aExponent += 1;
133 }
134 }
135
136 // If we have overflowed the type, return +/- infinity:
137 if (aExponent >= maxExponent) return @bitCast(infRep | resultSign);
138
139 if (aExponent <= 0) {
140 // Result is denormal; the exponent and round/sticky bits are zero.
141 // All we need to do is shift the significand and apply the correct sign.
142 aSignificand >>= @intCast(4 - aExponent);
143 return @bitCast(resultSign | aSignificand);
144 }
145
146 // Low three bits are round, guard, and sticky.
147 const roundGuardSticky = aSignificand & 0x7;
148
149 // Shift the significand into place, and mask off the integer bit, if it's implicit.
150 var result = (aSignificand >> 3) & significandMask;
151
152 // Insert the exponent and sign.
153 result |= @as(Z, @intCast(aExponent)) << significandBits;
154 result |= resultSign;
155
156 // Final rounding. The result may overflow to infinity, but that is the
157 // correct result in that case.
158 if (roundGuardSticky > 0x4) result += 1;
159 if (roundGuardSticky == 0x4) result += result & 1;
160
161 // Restore any explicit integer bit, if it was rounded off
162 if (significandBits != fractionalBits) {
163 if ((result >> significandBits) != 0) result |= integerBit;
164 }
165
166 return @bitCast(result);
167}
168
169test {
170 _ = @import("addf3_test.zig");
171}