| ... | @@ -3,12 +3,16 @@ | ... | @@ -3,12 +3,16 @@ |
| 3 | // https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/lib/builtins/fp_mul_impl.inc | 3 | // https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/lib/builtins/fp_mul_impl.inc |
| 4 | | 4 | |
| 5 | const std = @import("std"); | 5 | const std = @import("std"); |
| | 6 | const math = std.math; |
| 6 | const builtin = @import("builtin"); | 7 | const builtin = @import("builtin"); |
| 7 | const compiler_rt = @import("../compiler_rt.zig"); | 8 | const compiler_rt = @import("../compiler_rt.zig"); |
| 8 | | 9 | |
| 9 | pub fn __multf3(a: f128, b: f128) callconv(.C) f128 { | 10 | pub fn __multf3(a: f128, b: f128) callconv(.C) f128 { |
| 10 | return mulXf3(f128, a, b); | 11 | return mulXf3(f128, a, b); |
| 11 | } | 12 | } |
| | 13 | pub fn __mulxf3(a: f80, b: f80) callconv(.C) f80 { |
| | 14 | return mulXf3(f80, a, b); |
| | 15 | } |
| 12 | pub fn __muldf3(a: f64, b: f64) callconv(.C) f64 { | 16 | pub fn __muldf3(a: f64, b: f64) callconv(.C) f64 { |
| 13 | return mulXf3(f64, a, b); | 17 | return mulXf3(f64, a, b); |
| 14 | } | 18 | } |
| ... | @@ -29,30 +33,36 @@ pub fn __aeabi_dmul(a: f64, b: f64) callconv(.C) f64 { | ... | @@ -29,30 +33,36 @@ pub fn __aeabi_dmul(a: f64, b: f64) callconv(.C) f64 { |
| 29 | fn mulXf3(comptime T: type, a: T, b: T) T { | 33 | fn mulXf3(comptime T: type, a: T, b: T) T { |
| 30 | @setRuntimeSafety(builtin.is_test); | 34 | @setRuntimeSafety(builtin.is_test); |
| 31 | const typeWidth = @typeInfo(T).Float.bits; | 35 | const typeWidth = @typeInfo(T).Float.bits; |
| | 36 | const significandBits = math.floatMantissaBits(T); |
| | 37 | const fractionalBits = math.floatFractionalBits(T); |
| | 38 | const exponentBits = math.floatExponentBits(T); |
| | 39 | |
| 32 | const Z = std.meta.Int(.unsigned, typeWidth); | 40 | const Z = std.meta.Int(.unsigned, typeWidth); |
| 33 | | 41 | |
| 34 | const significandBits = std.math.floatMantissaBits(T); | 42 | // ZSignificand is large enough to contain the significand, including an explicit integer bit |
| 35 | const exponentBits = std.math.floatExponentBits(T); | 43 | const ZSignificand = PowerOfTwoSignificandZ(T); |
| | 44 | const ZSignificandBits = @typeInfo(ZSignificand).Int.bits; |
| 36 | | 45 | |
| | 46 | const roundBit = (1 << (ZSignificandBits - 1)); |
| 37 | const signBit = (@as(Z, 1) << (significandBits + exponentBits)); | 47 | const signBit = (@as(Z, 1) << (significandBits + exponentBits)); |
| 38 | const maxExponent = ((1 << exponentBits) - 1); | 48 | const maxExponent = ((1 << exponentBits) - 1); |
| 39 | const exponentBias = (maxExponent >> 1); | 49 | const exponentBias = (maxExponent >> 1); |
| 40 | | 50 | |
| 41 | const implicitBit = (@as(Z, 1) << significandBits); | 51 | const integerBit = (@as(ZSignificand, 1) << fractionalBits); |
| 42 | const quietBit = implicitBit >> 1; | 52 | const quietBit = integerBit >> 1; |
| 43 | const significandMask = implicitBit - 1; | 53 | const significandMask = (@as(Z, 1) << significandBits) - 1; |
| 44 | | 54 | |
| 45 | const absMask = signBit - 1; | 55 | const absMask = signBit - 1; |
| 46 | const exponentMask = absMask ^ significandMask; | 56 | const qnanRep = @bitCast(Z, math.nan(T)) | quietBit; |
| 47 | const qnanRep = exponentMask | quietBit; | 57 | const infRep = @bitCast(Z, math.inf(T)); |
| 48 | const infRep = @bitCast(Z, std.math.inf(T)); | 58 | const minNormalRep = @bitCast(Z, math.floatMin(T)); |
| 49 | | 59 | |
| 50 | const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent); | 60 | const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent); |
| 51 | const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent); | 61 | const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent); |
| 52 | const productSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit; | 62 | const productSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit; |
| 53 | | 63 | |
| 54 | var aSignificand: Z = @bitCast(Z, a) & significandMask; | 64 | var aSignificand: ZSignificand = @intCast(ZSignificand, @bitCast(Z, a) & significandMask); |
| 55 | var bSignificand: Z = @bitCast(Z, b) & significandMask; | 65 | var bSignificand: ZSignificand = @intCast(ZSignificand, @bitCast(Z, b) & significandMask); |
| 56 | var scale: i32 = 0; | 66 | var scale: i32 = 0; |
| 57 | | 67 | |
| 58 | // Detect if a or b is zero, denormal, infinity, or NaN. | 68 | // Detect if a or b is zero, denormal, infinity, or NaN. |
| ... | @@ -93,38 +103,40 @@ fn mulXf3(comptime T: type, a: T, b: T) T { | ... | @@ -93,38 +103,40 @@ fn mulXf3(comptime T: type, a: T, b: T) T { |
| 93 | // one or both of a or b is denormal, the other (if applicable) is a | 103 | // one or both of a or b is denormal, the other (if applicable) is a |
| 94 | // normal number. Renormalize one or both of a and b, and set scale to | 104 | // normal number. Renormalize one or both of a and b, and set scale to |
| 95 | // include the necessary exponent adjustment. | 105 | // include the necessary exponent adjustment. |
| 96 | if (aAbs < implicitBit) scale += normalize(T, &aSignificand); | 106 | if (aAbs < minNormalRep) scale += normalize(T, &aSignificand); |
| 97 | if (bAbs < implicitBit) scale += normalize(T, &bSignificand); | 107 | if (bAbs < minNormalRep) scale += normalize(T, &bSignificand); |
| 98 | } | 108 | } |
| 99 | | 109 | |
| 100 | // Or in the implicit significand bit. (If we fell through from the | 110 | // Or in the implicit significand bit. (If we fell through from the |
| 101 | // denormal path it was already set by normalize( ), but setting it twice | 111 | // denormal path it was already set by normalize( ), but setting it twice |
| 102 | // won't hurt anything.) | 112 | // won't hurt anything.) |
| 103 | aSignificand |= implicitBit; | 113 | aSignificand |= integerBit; |
| 104 | bSignificand |= implicitBit; | 114 | bSignificand |= integerBit; |
| 105 | | 115 | |
| 106 | // Get the significand of a*b. Before multiplying the significands, shift | 116 | // Get the significand of a*b. Before multiplying the significands, shift |
| 107 | // one of them left to left-align it in the field. Thus, the product will | 117 | // one of them left to left-align it in the field. Thus, the product will |
| 108 | // have (exponentBits + 2) integral digits, all but two of which must be | 118 | // have (exponentBits + 2) integral digits, all but two of which must be |
| 109 | // zero. Normalizing this result is just a conditional left-shift by one | 119 | // zero. Normalizing this result is just a conditional left-shift by one |
| 110 | // and bumping the exponent accordingly. | 120 | // and bumping the exponent accordingly. |
| 111 | var productHi: Z = undefined; | 121 | var productHi: ZSignificand = undefined; |
| 112 | var productLo: Z = undefined; | 122 | var productLo: ZSignificand = undefined; |
| 113 | wideMultiply(Z, aSignificand, bSignificand << exponentBits, &productHi, &productLo); | 123 | const left_align_shift = ZSignificandBits - fractionalBits - 1; |
| | 124 | wideMultiply(ZSignificand, aSignificand, bSignificand << left_align_shift, &productHi, &productLo); |
| 114 | | 125 | |
| 115 | var productExponent: i32 = @bitCast(i32, aExponent +% bExponent) -% exponentBias +% scale; | 126 | var productExponent: i32 = @intCast(i32, aExponent + bExponent) - exponentBias + scale; |
| 116 | | 127 | |
| 117 | // Normalize the significand, adjust exponent if needed. | 128 | // Normalize the significand, adjust exponent if needed. |
| 118 | if ((productHi & implicitBit) != 0) { | 129 | if ((productHi & integerBit) != 0) { |
| 119 | productExponent +%= 1; | 130 | productExponent +%= 1; |
| 120 | } else { | 131 | } else { |
| 121 | productHi = (productHi << 1) | (productLo >> (typeWidth - 1)); | 132 | productHi = (productHi << 1) | (productLo >> (ZSignificandBits - 1)); |
| 122 | productLo = productLo << 1; | 133 | productLo = productLo << 1; |
| 123 | } | 134 | } |
| 124 | | 135 | |
| 125 | // If we have overflowed the type, return +/- infinity. | 136 | // If we have overflowed the type, return +/- infinity. |
| 126 | if (productExponent >= maxExponent) return @bitCast(T, infRep | productSign); | 137 | if (productExponent >= maxExponent) return @bitCast(T, infRep | productSign); |
| 127 | | 138 | |
| | 139 | var result: Z = undefined; |
| 128 | if (productExponent <= 0) { | 140 | if (productExponent <= 0) { |
| 129 | // Result is denormal before rounding | 141 | // Result is denormal before rounding |
| 130 | // | 142 | // |
| ... | @@ -133,35 +145,49 @@ fn mulXf3(comptime T: type, a: T, b: T) T { | ... | @@ -133,35 +145,49 @@ fn mulXf3(comptime T: type, a: T, b: T) T { |
| 133 | // handle this case separately, but we make it a special case to | 145 | // handle this case separately, but we make it a special case to |
| 134 | // simplify the shift logic. | 146 | // simplify the shift logic. |
| 135 | const shift: u32 = @truncate(u32, @as(Z, 1) -% @bitCast(u32, productExponent)); | 147 | const shift: u32 = @truncate(u32, @as(Z, 1) -% @bitCast(u32, productExponent)); |
| 136 | if (shift >= typeWidth) return @bitCast(T, productSign); | 148 | if (shift >= ZSignificandBits) return @bitCast(T, productSign); |
| 137 | | 149 | |
| 138 | // Otherwise, shift the significand of the result so that the round | 150 | // Otherwise, shift the significand of the result so that the round |
| 139 | // bit is the high bit of productLo. | 151 | // bit is the high bit of productLo. |
| 140 | wideRightShiftWithSticky(Z, &productHi, &productLo, shift); | 152 | const sticky = wideShrWithTruncation(ZSignificand, &productHi, &productLo, shift); |
| | 153 | productLo |= @boolToInt(sticky); |
| | 154 | result = productHi; |
| 141 | } else { | 155 | } else { |
| 142 | // Result is normal before rounding; insert the exponent. | 156 | // Result is normal before rounding; insert the exponent. |
| 143 | productHi &= significandMask; | 157 | result = productHi & significandMask; |
| 144 | productHi |= @as(Z, @bitCast(u32, productExponent)) << significandBits; | 158 | result |= @intCast(Z, productExponent) << significandBits; |
| 145 | } | 159 | } |
| 146 | | 160 | |
| 147 | // Insert the sign of the result: | | |
| 148 | productHi |= productSign; | | |
| 149 | | | |
| 150 | // Final rounding. The final result may overflow to infinity, or underflow | 161 | // Final rounding. The final result may overflow to infinity, or underflow |
| 151 | // to zero, but those are the correct results in those cases. We use the | 162 | // to zero, but those are the correct results in those cases. We use the |
| 152 | // default IEEE-754 round-to-nearest, ties-to-even rounding mode. | 163 | // default IEEE-754 round-to-nearest, ties-to-even rounding mode. |
| 153 | if (productLo > signBit) productHi +%= 1; | 164 | if (productLo > roundBit) result +%= 1; |
| 154 | if (productLo == signBit) productHi +%= productHi & 1; | 165 | if (productLo == roundBit) result +%= result & 1; |
| 155 | return @bitCast(T, productHi); | 166 | |
| | 167 | // Restore any explicit integer bit, if it was rounded off |
| | 168 | if (significandBits != fractionalBits) { |
| | 169 | if ((result >> significandBits) != 0) result |= integerBit; |
| | 170 | } |
| | 171 | |
| | 172 | // Insert the sign of the result: |
| | 173 | result |= productSign; |
| | 174 | |
| | 175 | return @bitCast(T, result); |
| 156 | } | 176 | } |
| 157 | | 177 | |
| 158 | fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { | 178 | fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { |
| 159 | @setRuntimeSafety(builtin.is_test); | 179 | @setRuntimeSafety(builtin.is_test); |
| 160 | switch (Z) { | 180 | switch (Z) { |
| | 181 | u16 => { |
| | 182 | // 16x16 --> 32 bit multiply |
| | 183 | const product = @as(u32, a) * @as(u32, b); |
| | 184 | hi.* = @intCast(u16, product >> 16); |
| | 185 | lo.* = @truncate(u16, product); |
| | 186 | }, |
| 161 | u32 => { | 187 | u32 => { |
| 162 | // 32x32 --> 64 bit multiply | 188 | // 32x32 --> 64 bit multiply |
| 163 | const product = @as(u64, a) * @as(u64, b); | 189 | const product = @as(u64, a) * @as(u64, b); |
| 164 | hi.* = @truncate(u32, product >> 32); | 190 | hi.* = @intCast(u32, product >> 32); |
| 165 | lo.* = @truncate(u32, product); | 191 | lo.* = @truncate(u32, product); |
| 166 | }, | 192 | }, |
| 167 | u64 => { | 193 | u64 => { |
| ... | @@ -170,7 +196,7 @@ fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { | ... | @@ -170,7 +196,7 @@ fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { |
| 170 | return @truncate(u32, x); | 196 | return @truncate(u32, x); |
| 171 | } | 197 | } |
| 172 | fn hiWord(x: u64) u64 { | 198 | fn hiWord(x: u64) u64 { |
| 173 | return @truncate(u32, x >> 32); | 199 | return @intCast(u32, x >> 32); |
| 174 | } | 200 | } |
| 175 | }; | 201 | }; |
| 176 | // 64x64 -> 128 wide multiply for platforms that don't have such an operation; | 202 | // 64x64 -> 128 wide multiply for platforms that don't have such an operation; |
| ... | @@ -264,34 +290,45 @@ fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { | ... | @@ -264,34 +290,45 @@ fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { |
| 264 | } | 290 | } |
| 265 | } | 291 | } |
| 266 | | 292 | |
| 267 | fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T).Float.bits)) i32 { | 293 | /// Returns a power-of-two integer type that is large enough to contain |
| | 294 | /// the significand of T, including an explicit integer bit |
| | 295 | fn PowerOfTwoSignificandZ(comptime T: type) type { |
| | 296 | const bits = math.ceilPowerOfTwoAssert(u16, math.floatFractionalBits(T) + 1); |
| | 297 | return std.meta.Int(.unsigned, bits); |
| | 298 | } |
| | 299 | |
| | 300 | fn normalize(comptime T: type, significand: *PowerOfTwoSignificandZ(T)) i32 { |
| 268 | @setRuntimeSafety(builtin.is_test); | 301 | @setRuntimeSafety(builtin.is_test); |
| 269 | const Z = std.meta.Int(.unsigned, @typeInfo(T).Float.bits); | 302 | const Z = PowerOfTwoSignificandZ(T); |
| 270 | const significandBits = std.math.floatMantissaBits(T); | 303 | const integerBit = @as(Z, 1) << math.floatFractionalBits(T); |
| 271 | const implicitBit = @as(Z, 1) << significandBits; | | |
| 272 | | 304 | |
| 273 | const shift = @clz(Z, significand.*) - @clz(Z, implicitBit); | 305 | const shift = @clz(Z, significand.*) - @clz(Z, integerBit); |
| 274 | significand.* <<= @intCast(std.math.Log2Int(Z), shift); | 306 | significand.* <<= @intCast(math.Log2Int(Z), shift); |
| 275 | return @as(i32, 1) - shift; | 307 | return @as(i32, 1) - shift; |
| 276 | } | 308 | } |
| 277 | | 309 | |
| 278 | fn wideRightShiftWithSticky(comptime Z: type, hi: *Z, lo: *Z, count: u32) void { | 310 | // Returns `true` if the right shift is inexact (i.e. any bit shifted out is non-zero) |
| | 311 | // |
| | 312 | // This is analogous to an shr version of `@shlWithOverflow` |
| | 313 | fn wideShrWithTruncation(comptime Z: type, hi: *Z, lo: *Z, count: u32) bool { |
| 279 | @setRuntimeSafety(builtin.is_test); | 314 | @setRuntimeSafety(builtin.is_test); |
| 280 | const typeWidth = @typeInfo(Z).Int.bits; | 315 | const typeWidth = @typeInfo(Z).Int.bits; |
| 281 | const S = std.math.Log2Int(Z); | 316 | const S = math.Log2Int(Z); |
| | 317 | var inexact = false; |
| 282 | if (count < typeWidth) { | 318 | if (count < typeWidth) { |
| 283 | const sticky = @boolToInt((lo.* << @intCast(S, typeWidth -% count)) != 0); | 319 | inexact = (lo.* << @intCast(S, typeWidth -% count)) != 0; |
| 284 | lo.* = (hi.* << @intCast(S, typeWidth -% count)) | (lo.* >> @intCast(S, count)) | sticky; | 320 | lo.* = (hi.* << @intCast(S, typeWidth -% count)) | (lo.* >> @intCast(S, count)); |
| 285 | hi.* = hi.* >> @intCast(S, count); | 321 | hi.* = hi.* >> @intCast(S, count); |
| 286 | } else if (count < 2 * typeWidth) { | 322 | } else if (count < 2 * typeWidth) { |
| 287 | const sticky = @boolToInt((hi.* << @intCast(S, 2 * typeWidth -% count) | lo.*) != 0); | 323 | inexact = (hi.* << @intCast(S, 2 * typeWidth -% count) | lo.*) != 0; |
| 288 | lo.* = hi.* >> @intCast(S, count -% typeWidth) | sticky; | 324 | lo.* = hi.* >> @intCast(S, count -% typeWidth); |
| 289 | hi.* = 0; | 325 | hi.* = 0; |
| 290 | } else { | 326 | } else { |
| 291 | const sticky = @boolToInt((hi.* | lo.*) != 0); | 327 | inexact = (hi.* | lo.*) != 0; |
| 292 | lo.* = sticky; | 328 | lo.* = 0; |
| 293 | hi.* = 0; | 329 | hi.* = 0; |
| 294 | } | 330 | } |
| | 331 | return inexact; |
| 295 | } | 332 | } |
| 296 | | 333 | |
| 297 | test { | 334 | test { |
Cody Tapscott <topolarity@tapscott.me> 2022-04-18 20:22:43-07:00