| ... | ... | @@ -0,0 +1,202 @@ |
| 1 | const std = @import("std"); |
| 2 | const builtin = @import("builtin"); |
| 3 | const normalize = @import("divdf3.zig").normalize; |
| 4 | const wideMultiply = @import("divdf3.zig").wideMultiply; |
| 5 | |
| 6 | pub fn __divxf3(a: f80, b: f80) callconv(.C) f80 { |
| 7 | @setRuntimeSafety(builtin.is_test); |
| 8 | const T = f80; |
| 9 | const Z = std.meta.Int(.unsigned, @bitSizeOf(T)); |
| 10 | |
| 11 | const significandBits = std.math.floatMantissaBits(T); |
| 12 | const fractionalBits = std.math.floatFractionalBits(T); |
| 13 | const exponentBits = std.math.floatExponentBits(T); |
| 14 | |
| 15 | const signBit = (@as(Z, 1) << (significandBits + exponentBits)); |
| 16 | const maxExponent = ((1 << exponentBits) - 1); |
| 17 | const exponentBias = (maxExponent >> 1); |
| 18 | |
| 19 | const integerBit = (@as(Z, 1) << fractionalBits); |
| 20 | const quietBit = integerBit >> 1; |
| 21 | const significandMask = (@as(Z, 1) << significandBits) - 1; |
| 22 | |
| 23 | const absMask = signBit - 1; |
| 24 | const qnanRep = @bitCast(Z, std.math.nan(T)) | quietBit; |
| 25 | const infRep = @bitCast(Z, std.math.inf(T)); |
| 26 | |
| 27 | const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent); |
| 28 | const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent); |
| 29 | const quotientSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit; |
| 30 | |
| 31 | var aSignificand: Z = @bitCast(Z, a) & significandMask; |
| 32 | var bSignificand: Z = @bitCast(Z, b) & significandMask; |
| 33 | var scale: i32 = 0; |
| 34 | |
| 35 | // Detect if a or b is zero, denormal, infinity, or NaN. |
| 36 | if (aExponent -% 1 >= maxExponent - 1 or bExponent -% 1 >= maxExponent - 1) { |
| 37 | const aAbs: Z = @bitCast(Z, a) & absMask; |
| 38 | const bAbs: Z = @bitCast(Z, b) & absMask; |
| 39 | |
| 40 | // NaN / anything = qNaN |
| 41 | if (aAbs > infRep) return @bitCast(T, @bitCast(Z, a) | quietBit); |
| 42 | // anything / NaN = qNaN |
| 43 | if (bAbs > infRep) return @bitCast(T, @bitCast(Z, b) | quietBit); |
| 44 | |
| 45 | if (aAbs == infRep) { |
| 46 | // infinity / infinity = NaN |
| 47 | if (bAbs == infRep) { |
| 48 | return @bitCast(T, qnanRep); |
| 49 | } |
| 50 | // infinity / anything else = +/- infinity |
| 51 | else { |
| 52 | return @bitCast(T, aAbs | quotientSign); |
| 53 | } |
| 54 | } |
| 55 | |
| 56 | // anything else / infinity = +/- 0 |
| 57 | if (bAbs == infRep) return @bitCast(T, quotientSign); |
| 58 | |
| 59 | if (aAbs == 0) { |
| 60 | // zero / zero = NaN |
| 61 | if (bAbs == 0) { |
| 62 | return @bitCast(T, qnanRep); |
| 63 | } |
| 64 | // zero / anything else = +/- zero |
| 65 | else { |
| 66 | return @bitCast(T, quotientSign); |
| 67 | } |
| 68 | } |
| 69 | // anything else / zero = +/- infinity |
| 70 | if (bAbs == 0) return @bitCast(T, infRep | quotientSign); |
| 71 | |
| 72 | // one or both of a or b is denormal, the other (if applicable) is a |
| 73 | // normal number. Renormalize one or both of a and b, and set scale to |
| 74 | // include the necessary exponent adjustment. |
| 75 | if (aAbs < integerBit) scale +%= normalize(T, &aSignificand); |
| 76 | if (bAbs < integerBit) scale -%= normalize(T, &bSignificand); |
| 77 | } |
| 78 | var quotientExponent: i32 = @bitCast(i32, aExponent -% bExponent) +% scale; |
| 79 | |
| 80 | // Align the significand of b as a Q63 fixed-point number in the range |
| 81 | // [1, 2.0) and get a Q64 approximate reciprocal using a small minimax |
| 82 | // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This |
| 83 | // is accurate to about 3.5 binary digits. |
| 84 | const q63b = @intCast(u64, bSignificand); |
| 85 | var recip64 = @as(u64, 0x7504f333F9DE6484) -% q63b; |
| 86 | // 0x7504f333F9DE6484 / 2^64 + 1 = 3/4 + 1/sqrt(2) |
| 87 | |
| 88 | // Now refine the reciprocal estimate using a Newton-Raphson iteration: |
| 89 | // |
| 90 | // x1 = x0 * (2 - x0 * b) |
| 91 | // |
| 92 | // This doubles the number of correct binary digits in the approximation |
| 93 | // with each iteration. |
| 94 | var correction64: u64 = undefined; |
| 95 | correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); |
| 96 | recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); |
| 97 | correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); |
| 98 | recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); |
| 99 | correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); |
| 100 | recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); |
| 101 | correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); |
| 102 | recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); |
| 103 | correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); |
| 104 | recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); |
| 105 | |
| 106 | // The reciprocal may have overflowed to zero if the upper half of b is |
| 107 | // exactly 1.0. This would sabatoge the full-width final stage of the |
| 108 | // computation that follows, so we adjust the reciprocal down by one bit. |
| 109 | recip64 -%= 1; |
| 110 | |
| 111 | // We need to perform one more iteration to get us to 112 binary digits; |
| 112 | // The last iteration needs to happen with extra precision. |
| 113 | |
| 114 | // NOTE: This operation is equivalent to __multi3, which is not implemented |
| 115 | // in some architechures |
| 116 | var reciprocal: u128 = undefined; |
| 117 | var correction: u128 = undefined; |
| 118 | var dummy: u128 = undefined; |
| 119 | wideMultiply(u128, recip64, q63b, &dummy, &correction); |
| 120 | |
| 121 | correction = -%correction; |
| 122 | |
| 123 | const cHi = @truncate(u64, correction >> 64); |
| 124 | const cLo = @truncate(u64, correction); |
| 125 | |
| 126 | var r64cH: u128 = undefined; |
| 127 | var r64cL: u128 = undefined; |
| 128 | wideMultiply(u128, recip64, cHi, &dummy, &r64cH); |
| 129 | wideMultiply(u128, recip64, cLo, &dummy, &r64cL); |
| 130 | |
| 131 | reciprocal = r64cH + (r64cL >> 64); |
| 132 | |
| 133 | // Adjust the final 128-bit reciprocal estimate downward to ensure that it |
| 134 | // is strictly smaller than the infinitely precise exact reciprocal. Because |
| 135 | // the computation of the Newton-Raphson step is truncating at every step, |
| 136 | // this adjustment is small; most of the work is already done. |
| 137 | reciprocal -%= 2; |
| 138 | |
| 139 | // The numerical reciprocal is accurate to within 2^-112, lies in the |
| 140 | // interval [0.5, 1.0), and is strictly smaller than the true reciprocal |
| 141 | // of b. Multiplying a by this reciprocal thus gives a numerical q = a/b |
| 142 | // in Q127 with the following properties: |
| 143 | // |
| 144 | // 1. q < a/b |
| 145 | // 2. q is in the interval [0.5, 2.0) |
| 146 | // 3. The error in q is bounded away from 2^-63 (actually, we have |
| 147 | // many bits to spare, but this is all we need). |
| 148 | |
| 149 | // We need a 128 x 128 multiply high to compute q. |
| 150 | var quotient128: u128 = undefined; |
| 151 | var quotientLo: u128 = undefined; |
| 152 | wideMultiply(u128, aSignificand << 2, reciprocal, "ient128, "ientLo); |
| 153 | |
| 154 | // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0). |
| 155 | // Right shift the quotient if it falls in the [1,2) range and adjust the |
| 156 | // exponent accordingly. |
| 157 | var quotient: u64 = if (quotient128 < (integerBit << 1)) b: { |
| 158 | quotientExponent -= 1; |
| 159 | break :b @intCast(u64, quotient128); |
| 160 | } else @intCast(u64, quotient128 >> 1); |
| 161 | |
| 162 | // We are going to compute a residual of the form |
| 163 | // |
| 164 | // r = a - q*b |
| 165 | // |
| 166 | // We know from the construction of q that r satisfies: |
| 167 | // |
| 168 | // 0 <= r < ulp(q)*b |
| 169 | // |
| 170 | // If r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we |
| 171 | // already have the correct result. The exact halfway case cannot occur. |
| 172 | var residual: u64 = -%(quotient *% q63b); |
| 173 | |
| 174 | const writtenExponent = quotientExponent + exponentBias; |
| 175 | if (writtenExponent >= maxExponent) { |
| 176 | // If we have overflowed the exponent, return infinity. |
| 177 | return @bitCast(T, infRep | quotientSign); |
| 178 | } else if (writtenExponent < 1) { |
| 179 | if (writtenExponent == 0) { |
| 180 | // Check whether the rounded result is normal. |
| 181 | if (residual > (bSignificand >> 1)) { // round |
| 182 | if (quotient == (integerBit - 1)) // If the rounded result is normal, return it |
| 183 | return @bitCast(T, @bitCast(Z, std.math.floatMin(T)) | quotientSign); |
| 184 | } |
| 185 | } |
| 186 | // Flush denormals to zero. In the future, it would be nice to add |
| 187 | // code to round them correctly. |
| 188 | return @bitCast(T, quotientSign); |
| 189 | } else { |
| 190 | const round = @boolToInt(residual > (bSignificand >> 1)); |
| 191 | // Insert the exponent |
| 192 | var absResult = quotient | (@intCast(Z, writtenExponent) << significandBits); |
| 193 | // Round |
| 194 | absResult +%= round; |
| 195 | // Insert the sign and return |
| 196 | return @bitCast(T, absResult | quotientSign | integerBit); |
| 197 | } |
| 198 | } |
| 199 | |
| 200 | test { |
| 201 | _ = @import("divxf3_test.zig"); |
| 202 | } |
Cody Tapscott <topolarity@tapscott.me> 2022-04-25 16:31:43-07:00
lib