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