| ... | ... | @@ -0,0 +1,328 @@ |
| 1 | // Ported from: |
| 2 | // |
| 3 | // https://github.com/llvm/llvm-project/commit/d674d96bc56c0f377879d01c9d8dfdaaa7859cdb/compiler-rt/lib/builtins/divdf3.c |
| 4 | |
| 5 | const std = @import("std"); |
| 6 | const builtin = @import("builtin"); |
| 7 | |
| 8 | pub extern fn __divdf3(a: f64, b: f64) f64 { |
| 9 | @setRuntimeSafety(builtin.is_test); |
| 10 | const Z = @IntType(false, f64.bit_count); |
| 11 | const SignedZ = @IntType(true, f64.bit_count); |
| 12 | |
| 13 | const typeWidth = f64.bit_count; |
| 14 | const significandBits = std.math.floatMantissaBits(f64); |
| 15 | const exponentBits = std.math.floatExponentBits(f64); |
| 16 | |
| 17 | const signBit = (Z(1) << (significandBits + exponentBits)); |
| 18 | const maxExponent = ((1 << exponentBits) - 1); |
| 19 | const exponentBias = (maxExponent >> 1); |
| 20 | |
| 21 | const implicitBit = (Z(1) << significandBits); |
| 22 | const quietBit = implicitBit >> 1; |
| 23 | const significandMask = implicitBit - 1; |
| 24 | |
| 25 | const absMask = signBit - 1; |
| 26 | const exponentMask = absMask ^ significandMask; |
| 27 | const qnanRep = exponentMask | quietBit; |
| 28 | const infRep = @bitCast(Z, std.math.inf(f64)); |
| 29 | |
| 30 | const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent); |
| 31 | const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent); |
| 32 | const quotientSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit; |
| 33 | |
| 34 | var aSignificand: Z = @bitCast(Z, a) & significandMask; |
| 35 | var bSignificand: Z = @bitCast(Z, b) & significandMask; |
| 36 | var scale: i32 = 0; |
| 37 | |
| 38 | // Detect if a or b is zero, denormal, infinity, or NaN. |
| 39 | if (aExponent -% 1 >= maxExponent -% 1 or bExponent -% 1 >= maxExponent -% 1) { |
| 40 | const aAbs: Z = @bitCast(Z, a) & absMask; |
| 41 | const bAbs: Z = @bitCast(Z, b) & absMask; |
| 42 | |
| 43 | // NaN / anything = qNaN |
| 44 | if (aAbs > infRep) return @bitCast(f64, @bitCast(Z, a) | quietBit); |
| 45 | // anything / NaN = qNaN |
| 46 | if (bAbs > infRep) return @bitCast(f64, @bitCast(Z, b) | quietBit); |
| 47 | |
| 48 | if (aAbs == infRep) { |
| 49 | // infinity / infinity = NaN |
| 50 | if (bAbs == infRep) { |
| 51 | return @bitCast(f64, qnanRep); |
| 52 | } |
| 53 | // infinity / anything else = +/- infinity |
| 54 | else { |
| 55 | return @bitCast(f64, aAbs | quotientSign); |
| 56 | } |
| 57 | } |
| 58 | |
| 59 | // anything else / infinity = +/- 0 |
| 60 | if (bAbs == infRep) return @bitCast(f64, quotientSign); |
| 61 | |
| 62 | if (aAbs == 0) { |
| 63 | // zero / zero = NaN |
| 64 | if (bAbs == 0) { |
| 65 | return @bitCast(f64, qnanRep); |
| 66 | } |
| 67 | // zero / anything else = +/- zero |
| 68 | else { |
| 69 | return @bitCast(f64, quotientSign); |
| 70 | } |
| 71 | } |
| 72 | // anything else / zero = +/- infinity |
| 73 | if (bAbs == 0) return @bitCast(f64, infRep | quotientSign); |
| 74 | |
| 75 | // one or both of a or b is denormal, the other (if applicable) is a |
| 76 | // normal number. Renormalize one or both of a and b, and set scale to |
| 77 | // include the necessary exponent adjustment. |
| 78 | if (aAbs < implicitBit) scale +%= normalize(f64, &aSignificand); |
| 79 | if (bAbs < implicitBit) scale -%= normalize(f64, &bSignificand); |
| 80 | } |
| 81 | |
| 82 | // Or in the implicit significand bit. (If we fell through from the |
| 83 | // denormal path it was already set by normalize( ), but setting it twice |
| 84 | // won't hurt anything.) |
| 85 | aSignificand |= implicitBit; |
| 86 | bSignificand |= implicitBit; |
| 87 | var quotientExponent: i32 = @bitCast(i32, aExponent -% bExponent) +% scale; |
| 88 | |
| 89 | // Align the significand of b as a Q31 fixed-point number in the range |
| 90 | // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax |
| 91 | // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This |
| 92 | // is accurate to about 3.5 binary digits. |
| 93 | const q31b: u32 = @truncate(u32, bSignificand >> 21); |
| 94 | var recip32 = u32(0x7504f333) -% q31b; |
| 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, so after three iterations, we have about 28 binary |
| 102 | // digits of accuracy. |
| 103 | var correction32: u32 = undefined; |
| 104 | correction32 = @truncate(u32, ~(u64(recip32) *% q31b >> 32) +% 1); |
| 105 | recip32 = @truncate(u32, u64(recip32) *% correction32 >> 31); |
| 106 | correction32 = @truncate(u32, ~(u64(recip32) *% q31b >> 32) +% 1); |
| 107 | recip32 = @truncate(u32, u64(recip32) *% correction32 >> 31); |
| 108 | correction32 = @truncate(u32, ~(u64(recip32) *% q31b >> 32) +% 1); |
| 109 | recip32 = @truncate(u32, u64(recip32) *% correction32 >> 31); |
| 110 | |
| 111 | // recip32 might have overflowed to exactly zero in the preceding |
| 112 | // computation if the high word of b is exactly 1.0. This would sabotage |
| 113 | // the full-width final stage of the computation that follows, so we adjust |
| 114 | // recip32 downward by one bit. |
| 115 | recip32 -%= 1; |
| 116 | |
| 117 | // We need to perform one more iteration to get us to 56 binary digits; |
| 118 | // The last iteration needs to happen with extra precision. |
| 119 | const q63blo: u32 = @truncate(u32, bSignificand << 11); |
| 120 | var correction: u64 = undefined; |
| 121 | var reciprocal: u64 = undefined; |
| 122 | correction = ~(u64(recip32) *% q31b +% (u64(recip32) *% q63blo >> 32)) +% 1; |
| 123 | const cHi = @truncate(u32, correction >> 32); |
| 124 | const cLo = @truncate(u32, correction); |
| 125 | reciprocal = u64(recip32) *% cHi +% (u64(recip32) *% cLo >> 32); |
| 126 | |
| 127 | // We already adjusted the 32-bit estimate, now we need to adjust the final |
| 128 | // 64-bit reciprocal estimate downward to ensure that it is strictly smaller |
| 129 | // than the infinitely precise exact reciprocal. Because the computation |
| 130 | // of the Newton-Raphson step is truncating at every step, this adjustment |
| 131 | // is small; most of the work is already done. |
| 132 | reciprocal -%= 2; |
| 133 | |
| 134 | // The numerical reciprocal is accurate to within 2^-56, lies in the |
| 135 | // interval [0.5, 1.0), and is strictly smaller than the true reciprocal |
| 136 | // of b. Multiplying a by this reciprocal thus gives a numerical q = a/b |
| 137 | // in Q53 with the following properties: |
| 138 | // |
| 139 | // 1. q < a/b |
| 140 | // 2. q is in the interval [0.5, 2.0) |
| 141 | // 3. the error in q is bounded away from 2^-53 (actually, we have a |
| 142 | // couple of bits to spare, but this is all we need). |
| 143 | |
| 144 | // We need a 64 x 64 multiply high to compute q, which isn't a basic |
| 145 | // operation in C, so we need to be a little bit fussy. |
| 146 | var quotient: Z = undefined; |
| 147 | var quotientLo: Z = undefined; |
| 148 | wideMultiply(Z, aSignificand << 2, reciprocal, "ient, "ientLo); |
| 149 | |
| 150 | // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0). |
| 151 | // In either case, we are going to compute a residual of the form |
| 152 | // |
| 153 | // r = a - q*b |
| 154 | // |
| 155 | // We know from the construction of q that r satisfies: |
| 156 | // |
| 157 | // 0 <= r < ulp(q)*b |
| 158 | // |
| 159 | // if r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we |
| 160 | // already have the correct result. The exact halfway case cannot occur. |
| 161 | // We also take this time to right shift quotient if it falls in the [1,2) |
| 162 | // range and adjust the exponent accordingly. |
| 163 | var residual: Z = undefined; |
| 164 | if (quotient < (implicitBit << 1)) { |
| 165 | residual = (aSignificand << 53) -% quotient *% bSignificand; |
| 166 | quotientExponent -%= 1; |
| 167 | } else { |
| 168 | quotient >>= 1; |
| 169 | residual = (aSignificand << 52) -% quotient *% bSignificand; |
| 170 | } |
| 171 | |
| 172 | const writtenExponent = quotientExponent +% exponentBias; |
| 173 | |
| 174 | if (writtenExponent >= maxExponent) { |
| 175 | // If we have overflowed the exponent, return infinity. |
| 176 | return @bitCast(f64, infRep | quotientSign); |
| 177 | } else if (writtenExponent < 1) { |
| 178 | if (writtenExponent == 0) { |
| 179 | // Check whether the rounded result is normal. |
| 180 | const round = @boolToInt((residual << 1) > bSignificand); |
| 181 | // Clear the implicit bit. |
| 182 | var absResult = quotient & significandMask; |
| 183 | // Round. |
| 184 | absResult += round; |
| 185 | if ((absResult & ~significandMask) != 0) { |
| 186 | // The rounded result is normal; return it. |
| 187 | return @bitCast(f64, absResult | quotientSign); |
| 188 | } |
| 189 | } |
| 190 | // Flush denormals to zero. In the future, it would be nice to add |
| 191 | // code to round them correctly. |
| 192 | return @bitCast(f64, quotientSign); |
| 193 | } else { |
| 194 | const round = @boolToInt((residual << 1) > bSignificand); |
| 195 | // Clear the implicit bit |
| 196 | var absResult = quotient & significandMask; |
| 197 | // Insert the exponent |
| 198 | absResult |= @bitCast(Z, SignedZ(writtenExponent)) << significandBits; |
| 199 | // Round |
| 200 | absResult +%= round; |
| 201 | // Insert the sign and return |
| 202 | return @bitCast(f64, absResult | quotientSign); |
| 203 | } |
| 204 | } |
| 205 | |
| 206 | fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { |
| 207 | @setRuntimeSafety(builtin.is_test); |
| 208 | switch (Z) { |
| 209 | u32 => { |
| 210 | // 32x32 --> 64 bit multiply |
| 211 | const product = u64(a) * u64(b); |
| 212 | hi.* = @truncate(u32, product >> 32); |
| 213 | lo.* = @truncate(u32, product); |
| 214 | }, |
| 215 | u64 => { |
| 216 | const S = struct { |
| 217 | fn loWord(x: u64) u64 { |
| 218 | return @truncate(u32, x); |
| 219 | } |
| 220 | fn hiWord(x: u64) u64 { |
| 221 | return @truncate(u32, x >> 32); |
| 222 | } |
| 223 | }; |
| 224 | // 64x64 -> 128 wide multiply for platforms that don't have such an operation; |
| 225 | // many 64-bit platforms have this operation, but they tend to have hardware |
| 226 | // floating-point, so we don't bother with a special case for them here. |
| 227 | // Each of the component 32x32 -> 64 products |
| 228 | const plolo: u64 = S.loWord(a) * S.loWord(b); |
| 229 | const plohi: u64 = S.loWord(a) * S.hiWord(b); |
| 230 | const philo: u64 = S.hiWord(a) * S.loWord(b); |
| 231 | const phihi: u64 = S.hiWord(a) * S.hiWord(b); |
| 232 | // Sum terms that contribute to lo in a way that allows us to get the carry |
| 233 | const r0: u64 = S.loWord(plolo); |
| 234 | const r1: u64 = S.hiWord(plolo) +% S.loWord(plohi) +% S.loWord(philo); |
| 235 | lo.* = r0 +% (r1 << 32); |
| 236 | // Sum terms contributing to hi with the carry from lo |
| 237 | hi.* = S.hiWord(plohi) +% S.hiWord(philo) +% S.hiWord(r1) +% phihi; |
| 238 | }, |
| 239 | u128 => { |
| 240 | const Word_LoMask = u64(0x00000000ffffffff); |
| 241 | const Word_HiMask = u64(0xffffffff00000000); |
| 242 | const Word_FullMask = u64(0xffffffffffffffff); |
| 243 | const S = struct { |
| 244 | fn Word_1(x: u128) u64 { |
| 245 | return @truncate(u32, x >> 96); |
| 246 | } |
| 247 | fn Word_2(x: u128) u64 { |
| 248 | return @truncate(u32, x >> 64); |
| 249 | } |
| 250 | fn Word_3(x: u128) u64 { |
| 251 | return @truncate(u32, x >> 32); |
| 252 | } |
| 253 | fn Word_4(x: u128) u64 { |
| 254 | return @truncate(u32, x); |
| 255 | } |
| 256 | }; |
| 257 | // 128x128 -> 256 wide multiply for platforms that don't have such an operation; |
| 258 | // many 64-bit platforms have this operation, but they tend to have hardware |
| 259 | // floating-point, so we don't bother with a special case for them here. |
| 260 | |
| 261 | const product11: u64 = S.Word_1(a) * S.Word_1(b); |
| 262 | const product12: u64 = S.Word_1(a) * S.Word_2(b); |
| 263 | const product13: u64 = S.Word_1(a) * S.Word_3(b); |
| 264 | const product14: u64 = S.Word_1(a) * S.Word_4(b); |
| 265 | const product21: u64 = S.Word_2(a) * S.Word_1(b); |
| 266 | const product22: u64 = S.Word_2(a) * S.Word_2(b); |
| 267 | const product23: u64 = S.Word_2(a) * S.Word_3(b); |
| 268 | const product24: u64 = S.Word_2(a) * S.Word_4(b); |
| 269 | const product31: u64 = S.Word_3(a) * S.Word_1(b); |
| 270 | const product32: u64 = S.Word_3(a) * S.Word_2(b); |
| 271 | const product33: u64 = S.Word_3(a) * S.Word_3(b); |
| 272 | const product34: u64 = S.Word_3(a) * S.Word_4(b); |
| 273 | const product41: u64 = S.Word_4(a) * S.Word_1(b); |
| 274 | const product42: u64 = S.Word_4(a) * S.Word_2(b); |
| 275 | const product43: u64 = S.Word_4(a) * S.Word_3(b); |
| 276 | const product44: u64 = S.Word_4(a) * S.Word_4(b); |
| 277 | |
| 278 | const sum0: u128 = u128(product44); |
| 279 | const sum1: u128 = u128(product34) +% |
| 280 | u128(product43); |
| 281 | const sum2: u128 = u128(product24) +% |
| 282 | u128(product33) +% |
| 283 | u128(product42); |
| 284 | const sum3: u128 = u128(product14) +% |
| 285 | u128(product23) +% |
| 286 | u128(product32) +% |
| 287 | u128(product41); |
| 288 | const sum4: u128 = u128(product13) +% |
| 289 | u128(product22) +% |
| 290 | u128(product31); |
| 291 | const sum5: u128 = u128(product12) +% |
| 292 | u128(product21); |
| 293 | const sum6: u128 = u128(product11); |
| 294 | |
| 295 | const r0: u128 = (sum0 & Word_FullMask) +% |
| 296 | ((sum1 & Word_LoMask) << 32); |
| 297 | const r1: u128 = (sum0 >> 64) +% |
| 298 | ((sum1 >> 32) & Word_FullMask) +% |
| 299 | (sum2 & Word_FullMask) +% |
| 300 | ((sum3 << 32) & Word_HiMask); |
| 301 | |
| 302 | lo.* = r0 +% (r1 << 64); |
| 303 | hi.* = (r1 >> 64) +% |
| 304 | (sum1 >> 96) +% |
| 305 | (sum2 >> 64) +% |
| 306 | (sum3 >> 32) +% |
| 307 | sum4 +% |
| 308 | (sum5 << 32) +% |
| 309 | (sum6 << 64); |
| 310 | }, |
| 311 | else => @compileError("unsupported"), |
| 312 | } |
| 313 | } |
| 314 | |
| 315 | fn normalize(comptime T: type, significand: *@IntType(false, T.bit_count)) i32 { |
| 316 | @setRuntimeSafety(builtin.is_test); |
| 317 | const Z = @IntType(false, T.bit_count); |
| 318 | const significandBits = std.math.floatMantissaBits(T); |
| 319 | const implicitBit = Z(1) << significandBits; |
| 320 | |
| 321 | const shift = @clz(significand.*) - @clz(implicitBit); |
| 322 | significand.* <<= @intCast(std.math.Log2Int(Z), shift); |
| 323 | return 1 - shift; |
| 324 | } |
| 325 | |
| 326 | test "import divdf3" { |
| 327 | _ = @import("divdf3_test.zig"); |
| 328 | } |
vegecode <justin.b.alexander1@gmail.com> 2019-04-04 21:46:53-05:00
Andrew Kelley <andrew@ziglang.org> 2019-04-05 18:10:19-04:00
std