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