| ... | ... | @@ -0,0 +1,200 @@ |
| 1 | // Ported from: |
| 2 | // |
| 3 | // https://github.com/llvm/llvm-project/commit/d674d96bc56c0f377879d01c9d8dfdaaa7859cdb/compiler-rt/lib/builtins/divsf3.c |
| 4 | |
| 5 | const std = @import("std"); |
| 6 | |
| 7 | pub extern fn __divsf3(a: f32, b: f32) f32 { |
| 8 | const Z = @IntType(false, f32.bit_count); |
| 9 | |
| 10 | const typeWidth = f32.bit_count; |
| 11 | const significandBits = std.math.floatMantissaBits(f32); |
| 12 | const exponentBits = std.math.floatExponentBits(f32); |
| 13 | |
| 14 | const signBit = (Z(1) << (significandBits + exponentBits)); |
| 15 | const maxExponent = ((1 << exponentBits) - 1); |
| 16 | const exponentBias = (maxExponent >> 1); |
| 17 | |
| 18 | const implicitBit = (Z(1) << significandBits); |
| 19 | const quietBit = implicitBit >> 1; |
| 20 | const significandMask = implicitBit - 1; |
| 21 | |
| 22 | const absMask = signBit - 1; |
| 23 | const exponentMask = absMask ^ significandMask; |
| 24 | const qnanRep = exponentMask | quietBit; |
| 25 | const infRep = @bitCast(Z, std.math.inf(f32)); |
| 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(f32, @bitCast(Z, a) | quietBit); |
| 42 | // anything * NaN = qNaN |
| 43 | if (bAbs > infRep) return @bitCast(f32, @bitCast(Z, b) | quietBit); |
| 44 | |
| 45 | if (aAbs == infRep) { |
| 46 | // infinity * non-zero = +/- infinity |
| 47 | if (bAbs != 0) { |
| 48 | return @bitCast(f32, aAbs | quotientSign); |
| 49 | } else { |
| 50 | // infinity * zero = NaN |
| 51 | return @bitCast(f32, qnanRep); |
| 52 | } |
| 53 | } |
| 54 | |
| 55 | if (bAbs == infRep) { |
| 56 | //? non-zero * infinity = +/- infinity |
| 57 | if (aAbs != 0) { |
| 58 | return @bitCast(f32, bAbs | quotientSign); |
| 59 | } else { |
| 60 | // zero * infinity = NaN |
| 61 | return @bitCast(f32, qnanRep); |
| 62 | } |
| 63 | } |
| 64 | |
| 65 | // zero * anything = +/- zero |
| 66 | if (aAbs == 0) return @bitCast(f32, quotientSign); |
| 67 | // anything * zero = +/- zero |
| 68 | if (bAbs == 0) return @bitCast(f32, quotientSign); |
| 69 | |
| 70 | // 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 |
| 72 | // include the necessary exponent adjustment. |
| 73 | if (aAbs < implicitBit) scale +%= normalize(f32, &aSignificand); |
| 74 | if (bAbs < implicitBit) scale +%= normalize(f32, &bSignificand); |
| 75 | } |
| 76 | |
| 77 | // Or in the implicit significand bit. (If we fell through from the |
| 78 | // denormal path it was already set by normalize( ), but setting it twice |
| 79 | // won't hurt anything.) |
| 80 | aSignificand |= implicitBit; |
| 81 | bSignificand |= implicitBit; |
| 82 | var quotientExponent: i32 = @bitCast(i32, aExponent -% bExponent) +% scale; |
| 83 | |
| 84 | // Align the significand of b as a Q31 fixed-point number in the range |
| 85 | // [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 |
| 87 | // is accurate to about 3.5 binary digits. |
| 88 | const q31b = switch (f32) { |
| 89 | f32 => bSignificand << 8, |
| 90 | f64 => bSignificand >> 21, |
| 91 | else => @compileError("Type not implemented."), |
| 92 | }; |
| 93 | var reciprocal = u32(0x7504f333) -% q31b; |
| 94 | |
| 95 | // Now refine the reciprocal estimate using a Newton-Raphson iteration: |
| 96 | // |
| 97 | // x1 = x0 * (2 - x0 * b) |
| 98 | // |
| 99 | // This doubles the number of correct binary digits in the approximation |
| 100 | // with each iteration, so after three iterations, we have about 28 binary |
| 101 | // digits of accuracy. |
| 102 | var correction: u32 = undefined; |
| 103 | correction = @truncate(u32, ~(u64(reciprocal) *% q31b >> 32) +% 1); |
| 104 | reciprocal = @truncate(u32, u64(reciprocal) *% correction >> 31); |
| 105 | correction = @truncate(u32, ~(u64(reciprocal) *% q31b >> 32) +% 1); |
| 106 | reciprocal = @truncate(u32, u64(reciprocal) *% correction >> 31); |
| 107 | correction = @truncate(u32, ~(u64(reciprocal) *% q31b >> 32) +% 1); |
| 108 | reciprocal = @truncate(u32, u64(reciprocal) *% correction >> 31); |
| 109 | |
| 110 | // Exhaustive testing shows that the error in reciprocal after three steps |
| 111 | // is in the interval [-0x1.f58108p-31, 0x1.d0e48cp-29], in line with our |
| 112 | // expectations. We bump the reciprocal by a tiny value to force the error |
| 113 | // to be strictly positive (in the range [0x1.4fdfp-37,0x1.287246p-29], to |
| 114 | // be specific). This also causes 1/1 to give a sensible approximation |
| 115 | // instead of zero (due to overflow). |
| 116 | reciprocal -%= 2; |
| 117 | |
| 118 | // The numerical reciprocal is accurate to within 2^-28, lies in the |
| 119 | // interval [0x1.000000eep-1, 0x1.fffffffcp-1], and is strictly smaller |
| 120 | // than the true reciprocal of b. Multiplying a by this reciprocal thus |
| 121 | // gives a numerical q = a/b in Q24 with the following properties: |
| 122 | // |
| 123 | // 1. q < a/b |
| 124 | // 2. q is in the interval [0x1.000000eep-1, 0x1.fffffffcp0) |
| 125 | // 3. the error in q is at most 2^-24 + 2^-27 -- the 2^24 term comes |
| 126 | // from the fact that we truncate the product, and the 2^27 term |
| 127 | // is the error in the reciprocal of b scaled by the maximum |
| 128 | // possible value of a. As a consequence of this error bound, |
| 129 | // either q or nextafter(q) is the correctly rounded |
| 130 | var quotient: Z = @truncate(u32, u64(reciprocal) *% (aSignificand << 1) >> 32); |
| 131 | |
| 132 | // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0). |
| 133 | // In either case, we are going to compute a residual of the form |
| 134 | // |
| 135 | // r = a - q*b |
| 136 | // |
| 137 | // We know from the construction of q that r satisfies: |
| 138 | // |
| 139 | // 0 <= r < ulp(q)*b |
| 140 | // |
| 141 | // if r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we |
| 142 | // already have the correct result. The exact halfway case cannot occur. |
| 143 | // We also take this time to right shift quotient if it falls in the [1,2) |
| 144 | // range and adjust the exponent accordingly. |
| 145 | var residual: Z = undefined; |
| 146 | if (quotient < (implicitBit << 1)) { |
| 147 | residual = (aSignificand << 24) -% quotient *% bSignificand; |
| 148 | quotientExponent -%= 1; |
| 149 | } else { |
| 150 | quotient >>= 1; |
| 151 | residual = (aSignificand << 23) -% quotient *% bSignificand; |
| 152 | } |
| 153 | |
| 154 | const writtenExponent = quotientExponent +% exponentBias; |
| 155 | |
| 156 | if (writtenExponent >= maxExponent) { |
| 157 | // If we have overflowed the exponent, return infinity. |
| 158 | return @bitCast(f32, infRep | quotientSign); |
| 159 | } else if (writtenExponent < 1) { |
| 160 | if (writtenExponent == 0) { |
| 161 | // Check whether the rounded result is normal. |
| 162 | const round = @boolToInt((residual << 1) > bSignificand); |
| 163 | // Clear the implicit bit. |
| 164 | var absResult = quotient & significandMask; |
| 165 | // Round. |
| 166 | absResult += round; |
| 167 | if ((absResult & ~significandMask) > 0) { |
| 168 | // The rounded result is normal; return it. |
| 169 | return @bitCast(f32, absResult | quotientSign); |
| 170 | } |
| 171 | } |
| 172 | // Flush denormals to zero. In the future, it would be nice to add |
| 173 | // code to round them correctly. |
| 174 | return @bitCast(f32, quotientSign); |
| 175 | } else { |
| 176 | const round = @boolToInt((residual << 1) > bSignificand); |
| 177 | // Clear the implicit bit |
| 178 | var absResult = quotient & significandMask; |
| 179 | // Insert the exponent |
| 180 | absResult |= @bitCast(Z, writtenExponent) << significandBits; |
| 181 | // Round |
| 182 | absResult +%= round; |
| 183 | // Insert the sign and return |
| 184 | return @bitCast(f32, absResult | quotientSign); |
| 185 | } |
| 186 | } |
| 187 | |
| 188 | fn normalize(comptime T: type, significand: *@IntType(false, T.bit_count)) i32 { |
| 189 | const Z = @IntType(false, T.bit_count); |
| 190 | const significandBits = std.math.floatMantissaBits(T); |
| 191 | const implicitBit = Z(1) << significandBits; |
| 192 | |
| 193 | const shift = @clz(significand.*) - @clz(implicitBit); |
| 194 | significand.* <<= @intCast(std.math.Log2Int(Z), shift); |
| 195 | return 1 - shift; |
| 196 | } |
| 197 | |
| 198 | test "import divsf3" { |
| 199 | _ = @import("divsf3_test.zig"); |
| 200 | } |
vegecode <justin.b.alexander1@gmail.com> 2019-04-03 13:12:31-05:00
Andrew Kelley <andrew@ziglang.org> 2019-04-04 15:38:09-04:00
std