| ... | @@ -3,8 +3,10 @@ | ... | @@ -3,8 +3,10 @@ |
| 3 | // https://github.com/llvm/llvm-project/commit/d674d96bc56c0f377879d01c9d8dfdaaa7859cdb/compiler-rt/lib/builtins/divsf3.c | 3 | // https://github.com/llvm/llvm-project/commit/d674d96bc56c0f377879d01c9d8dfdaaa7859cdb/compiler-rt/lib/builtins/divsf3.c |
| 4 | | 4 | |
| 5 | const std = @import("std"); | 5 | const std = @import("std"); |
| | 6 | const builtin = @import("builtin"); |
| 6 | | 7 | |
| 7 | pub extern fn __divsf3(a: f32, b: f32) f32 { | 8 | pub extern fn __divsf3(a: f32, b: f32) f32 { |
| | 9 | @setRuntimeSafety(builtin.is_test); |
| 8 | const Z = @IntType(false, f32.bit_count); | 10 | const Z = @IntType(false, f32.bit_count); |
| 9 | | 11 | |
| 10 | const typeWidth = f32.bit_count; | 12 | const typeWidth = f32.bit_count; |
| ... | @@ -37,41 +39,43 @@ pub extern fn __divsf3(a: f32, b: f32) f32 { | ... | @@ -37,41 +39,43 @@ pub extern fn __divsf3(a: f32, b: f32) f32 { |
| 37 | const aAbs: Z = @bitCast(Z, a) & absMask; | 39 | const aAbs: Z = @bitCast(Z, a) & absMask; |
| 38 | const bAbs: Z = @bitCast(Z, b) & absMask; | 40 | const bAbs: Z = @bitCast(Z, b) & absMask; |
| 39 | | 41 | |
| 40 | // NaN * anything = qNaN | 42 | // NaN / anything = qNaN |
| 41 | if (aAbs > infRep) return @bitCast(f32, @bitCast(Z, a) | quietBit); | 43 | if (aAbs > infRep) return @bitCast(f32, @bitCast(Z, a) | quietBit); |
| 42 | // anything * NaN = qNaN | 44 | // anything / NaN = qNaN |
| 43 | if (bAbs > infRep) return @bitCast(f32, @bitCast(Z, b) | quietBit); | 45 | if (bAbs > infRep) return @bitCast(f32, @bitCast(Z, b) | quietBit); |
| 44 | | 46 | |
| 45 | if (aAbs == infRep) { | 47 | if (aAbs == infRep) { |
| 46 | // infinity * non-zero = +/- infinity | 48 | // infinity / infinity = NaN |
| 47 | if (bAbs != 0) { | 49 | if (bAbs == infRep) { |
| 48 | return @bitCast(f32, aAbs | quotientSign); | | |
| 49 | } else { | | |
| 50 | // infinity * zero = NaN | | |
| 51 | return @bitCast(f32, qnanRep); | 50 | return @bitCast(f32, qnanRep); |
| 52 | } | 51 | } |
| | 52 | // infinity / anything else = +/- infinity |
| | 53 | else { |
| | 54 | return @bitCast(f32, aAbs | quotientSign); |
| | 55 | } |
| 53 | } | 56 | } |
| 54 | | 57 | |
| 55 | if (bAbs == infRep) { | 58 | // anything else / infinity = +/- 0 |
| 56 | //? non-zero * infinity = +/- infinity | 59 | if (bAbs == infRep) return @bitCast(f32, quotientSign); |
| 57 | if (aAbs != 0) { | 60 | |
| 58 | return @bitCast(f32, bAbs | quotientSign); | 61 | if (aAbs == 0) { |
| 59 | } else { | 62 | // zero / zero = NaN |
| 60 | // zero * infinity = NaN | 63 | if (bAbs == 0) { |
| 61 | return @bitCast(f32, qnanRep); | 64 | return @bitCast(f32, qnanRep); |
| 62 | } | 65 | } |
| | 66 | // zero / anything else = +/- zero |
| | 67 | else { |
| | 68 | return @bitCast(f32, quotientSign); |
| | 69 | } |
| 63 | } | 70 | } |
| 64 | | 71 | // anything else / zero = +/- infinity |
| 65 | // zero * anything = +/- zero | 72 | if (bAbs == 0) return @bitCast(f32, infRep | quotientSign); |
| 66 | if (aAbs == 0) return @bitCast(f32, quotientSign); | | |
| 67 | // anything * zero = +/- zero | | |
| 68 | if (bAbs == 0) return @bitCast(f32, quotientSign); | | |
| 69 | | 73 | |
| 70 | // one or both of a or b is denormal, the other (if applicable) is a | 74 | // one or both of a or b is denormal, the other (if applicable) is a |
| 71 | // normal number. Renormalize one or both of a and b, and set scale to | 75 | // normal number. Renormalize one or both of a and b, and set scale to |
| 72 | // include the necessary exponent adjustment. | 76 | // include the necessary exponent adjustment. |
| 73 | if (aAbs < implicitBit) scale +%= normalize(f32, &aSignificand); | 77 | if (aAbs < implicitBit) scale +%= normalize(f32, &aSignificand); |
| 74 | if (bAbs < implicitBit) scale +%= normalize(f32, &bSignificand); | 78 | if (bAbs < implicitBit) scale -%= normalize(f32, &bSignificand); |
| 75 | } | 79 | } |
| 76 | | 80 | |
| 77 | // Or in the implicit significand bit. (If we fell through from the | 81 | // Or in the implicit significand bit. (If we fell through from the |
| ... | @@ -85,11 +89,7 @@ pub extern fn __divsf3(a: f32, b: f32) f32 { | ... | @@ -85,11 +89,7 @@ pub extern fn __divsf3(a: f32, b: f32) f32 { |
| 85 | // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax | 89 | // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax |
| 86 | // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This | 90 | // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This |
| 87 | // is accurate to about 3.5 binary digits. | 91 | // is accurate to about 3.5 binary digits. |
| 88 | const q31b = switch (f32) { | 92 | const q31b = bSignificand << 8; |
| 89 | f32 => bSignificand << 8, | | |
| 90 | f64 => bSignificand >> 21, | | |
| 91 | else => @compileError("Type not implemented."), | | |
| 92 | }; | | |
| 93 | var reciprocal = u32(0x7504f333) -% q31b; | 93 | var reciprocal = u32(0x7504f333) -% q31b; |
| 94 | | 94 | |
| 95 | // Now refine the reciprocal estimate using a Newton-Raphson iteration: | 95 | // Now refine the reciprocal estimate using a Newton-Raphson iteration: |
| ... | @@ -186,6 +186,7 @@ pub extern fn __divsf3(a: f32, b: f32) f32 { | ... | @@ -186,6 +186,7 @@ pub extern fn __divsf3(a: f32, b: f32) f32 { |
| 186 | } | 186 | } |
| 187 | | 187 | |
| 188 | fn normalize(comptime T: type, significand: *@IntType(false, T.bit_count)) i32 { | 188 | fn normalize(comptime T: type, significand: *@IntType(false, T.bit_count)) i32 { |
| | 189 | @setRuntimeSafety(builtin.is_test); |
| 189 | const Z = @IntType(false, T.bit_count); | 190 | const Z = @IntType(false, T.bit_count); |
| 190 | const significandBits = std.math.floatMantissaBits(T); | 191 | const significandBits = std.math.floatMantissaBits(T); |
| 191 | const implicitBit = Z(1) << significandBits; | 192 | const implicitBit = Z(1) << significandBits; |