blob 1e9171c8 (12412B) - Raw
1 // Ported from: 2 // 3 // https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/lib/builtins/fp_mul_impl.inc 4 5 const std = @import("std"); 6 const builtin = @import("builtin"); 7 const compiler_rt = @import("../compiler_rt.zig"); 8 9 pub fn __multf3(a: f128, b: f128) callconv(.C) f128 { 10 return mulXf3(f128, a, b); 11 } 12 pub fn __muldf3(a: f64, b: f64) callconv(.C) f64 { 13 return mulXf3(f64, a, b); 14 } 15 pub fn __mulsf3(a: f32, b: f32) callconv(.C) f32 { 16 return mulXf3(f32, a, b); 17 } 18 19 pub fn __aeabi_fmul(a: f32, b: f32) callconv(.C) f32 { 20 @setRuntimeSafety(false); 21 return @call(.{ .modifier = .always_inline }, __mulsf3, .{ a, b }); 22 } 23 24 pub fn __aeabi_dmul(a: f64, b: f64) callconv(.C) f64 { 25 @setRuntimeSafety(false); 26 return @call(.{ .modifier = .always_inline }, __muldf3, .{ a, b }); 27 } 28 29 fn mulXf3(comptime T: type, a: T, b: T) T { 30 @setRuntimeSafety(builtin.is_test); 31 const typeWidth = @typeInfo(T).Float.bits; 32 const Z = std.meta.Int(.unsigned, typeWidth); 33 34 const significandBits = std.math.floatMantissaBits(T); 35 const exponentBits = std.math.floatExponentBits(T); 36 37 const signBit = (@as(Z, 1) << (significandBits + exponentBits)); 38 const maxExponent = ((1 << exponentBits) - 1); 39 const exponentBias = (maxExponent >> 1); 40 41 const implicitBit = (@as(Z, 1) << significandBits); 42 const quietBit = implicitBit >> 1; 43 const significandMask = implicitBit - 1; 44 45 const absMask = signBit - 1; 46 const exponentMask = absMask ^ significandMask; 47 const qnanRep = exponentMask | quietBit; 48 const infRep = @bitCast(Z, std.math.inf(T)); 49 50 const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent); 51 const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent); 52 const productSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit; 53 54 var aSignificand: Z = @bitCast(Z, a) & significandMask; 55 var bSignificand: Z = @bitCast(Z, b) & significandMask; 56 var scale: i32 = 0; 57 58 // Detect if a or b is zero, denormal, infinity, or NaN. 59 if (aExponent -% 1 >= maxExponent -% 1 or bExponent -% 1 >= maxExponent -% 1) { 60 const aAbs: Z = @bitCast(Z, a) & absMask; 61 const bAbs: Z = @bitCast(Z, b) & absMask; 62 63 // NaN * anything = qNaN 64 if (aAbs > infRep) return @bitCast(T, @bitCast(Z, a) | quietBit); 65 // anything * NaN = qNaN 66 if (bAbs > infRep) return @bitCast(T, @bitCast(Z, b) | quietBit); 67 68 if (aAbs == infRep) { 69 // infinity * non-zero = +/- infinity 70 if (bAbs != 0) { 71 return @bitCast(T, aAbs | productSign); 72 } else { 73 // infinity * zero = NaN 74 return @bitCast(T, qnanRep); 75 } 76 } 77 78 if (bAbs == infRep) { 79 //? non-zero * infinity = +/- infinity 80 if (aAbs != 0) { 81 return @bitCast(T, bAbs | productSign); 82 } else { 83 // zero * infinity = NaN 84 return @bitCast(T, qnanRep); 85 } 86 } 87 88 // zero * anything = +/- zero 89 if (aAbs == 0) return @bitCast(T, productSign); 90 // anything * zero = +/- zero 91 if (bAbs == 0) return @bitCast(T, productSign); 92 93 // 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 95 // include the necessary exponent adjustment. 96 if (aAbs < implicitBit) scale += normalize(T, &aSignificand); 97 if (bAbs < implicitBit) scale += normalize(T, &bSignificand); 98 } 99 100 // 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 102 // won't hurt anything.) 103 aSignificand |= implicitBit; 104 bSignificand |= implicitBit; 105 106 // 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 108 // 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 110 // and bumping the exponent accordingly. 111 var productHi: Z = undefined; 112 var productLo: Z = undefined; 113 wideMultiply(Z, aSignificand, bSignificand << exponentBits, &productHi, &productLo); 114 115 var productExponent: i32 = @bitCast(i32, aExponent +% bExponent) -% exponentBias +% scale; 116 117 // Normalize the significand, adjust exponent if needed. 118 if ((productHi & implicitBit) != 0) { 119 productExponent +%= 1; 120 } else { 121 productHi = (productHi << 1) | (productLo >> (typeWidth - 1)); 122 productLo = productLo << 1; 123 } 124 125 // If we have overflowed the type, return +/- infinity. 126 if (productExponent >= maxExponent) return @bitCast(T, infRep | productSign); 127 128 if (productExponent <= 0) { 129 // Result is denormal before rounding 130 // 131 // If the result is so small that it just underflows to zero, return 132 // a zero of the appropriate sign. Mathematically there is no need to 133 // handle this case separately, but we make it a special case to 134 // simplify the shift logic. 135 const shift: u32 = @truncate(u32, @as(Z, 1) -% @bitCast(u32, productExponent)); 136 if (shift >= typeWidth) return @bitCast(T, productSign); 137 138 // Otherwise, shift the significand of the result so that the round 139 // bit is the high bit of productLo. 140 wideRightShiftWithSticky(Z, &productHi, &productLo, shift); 141 } else { 142 // Result is normal before rounding; insert the exponent. 143 productHi &= significandMask; 144 productHi |= @as(Z, @bitCast(u32, productExponent)) << significandBits; 145 } 146 147 // Insert the sign of the result: 148 productHi |= productSign; 149 150 // 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 152 // default IEEE-754 round-to-nearest, ties-to-even rounding mode. 153 if (productLo > signBit) productHi +%= 1; 154 if (productLo == signBit) productHi +%= productHi & 1; 155 return @bitCast(T, productHi); 156 } 157 158 fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { 159 @setRuntimeSafety(builtin.is_test); 160 switch (Z) { 161 u32 => { 162 // 32x32 --> 64 bit multiply 163 const product = @as(u64, a) * @as(u64, b); 164 hi.* = @truncate(u32, product >> 32); 165 lo.* = @truncate(u32, product); 166 }, 167 u64 => { 168 const S = struct { 169 fn loWord(x: u64) u64 { 170 return @truncate(u32, x); 171 } 172 fn hiWord(x: u64) u64 { 173 return @truncate(u32, x >> 32); 174 } 175 }; 176 // 64x64 -> 128 wide multiply for platforms that don't have such an operation; 177 // many 64-bit platforms have this operation, but they tend to have hardware 178 // floating-point, so we don't bother with a special case for them here. 179 // Each of the component 32x32 -> 64 products 180 const plolo: u64 = S.loWord(a) * S.loWord(b); 181 const plohi: u64 = S.loWord(a) * S.hiWord(b); 182 const philo: u64 = S.hiWord(a) * S.loWord(b); 183 const phihi: u64 = S.hiWord(a) * S.hiWord(b); 184 // Sum terms that contribute to lo in a way that allows us to get the carry 185 const r0: u64 = S.loWord(plolo); 186 const r1: u64 = S.hiWord(plolo) +% S.loWord(plohi) +% S.loWord(philo); 187 lo.* = r0 +% (r1 << 32); 188 // Sum terms contributing to hi with the carry from lo 189 hi.* = S.hiWord(plohi) +% S.hiWord(philo) +% S.hiWord(r1) +% phihi; 190 }, 191 u128 => { 192 const Word_LoMask = @as(u64, 0x00000000ffffffff); 193 const Word_HiMask = @as(u64, 0xffffffff00000000); 194 const Word_FullMask = @as(u64, 0xffffffffffffffff); 195 const S = struct { 196 fn Word_1(x: u128) u64 { 197 return @truncate(u32, x >> 96); 198 } 199 fn Word_2(x: u128) u64 { 200 return @truncate(u32, x >> 64); 201 } 202 fn Word_3(x: u128) u64 { 203 return @truncate(u32, x >> 32); 204 } 205 fn Word_4(x: u128) u64 { 206 return @truncate(u32, x); 207 } 208 }; 209 // 128x128 -> 256 wide multiply for platforms that don't have such an operation; 210 // many 64-bit platforms have this operation, but they tend to have hardware 211 // floating-point, so we don't bother with a special case for them here. 212 213 const product11: u64 = S.Word_1(a) * S.Word_1(b); 214 const product12: u64 = S.Word_1(a) * S.Word_2(b); 215 const product13: u64 = S.Word_1(a) * S.Word_3(b); 216 const product14: u64 = S.Word_1(a) * S.Word_4(b); 217 const product21: u64 = S.Word_2(a) * S.Word_1(b); 218 const product22: u64 = S.Word_2(a) * S.Word_2(b); 219 const product23: u64 = S.Word_2(a) * S.Word_3(b); 220 const product24: u64 = S.Word_2(a) * S.Word_4(b); 221 const product31: u64 = S.Word_3(a) * S.Word_1(b); 222 const product32: u64 = S.Word_3(a) * S.Word_2(b); 223 const product33: u64 = S.Word_3(a) * S.Word_3(b); 224 const product34: u64 = S.Word_3(a) * S.Word_4(b); 225 const product41: u64 = S.Word_4(a) * S.Word_1(b); 226 const product42: u64 = S.Word_4(a) * S.Word_2(b); 227 const product43: u64 = S.Word_4(a) * S.Word_3(b); 228 const product44: u64 = S.Word_4(a) * S.Word_4(b); 229 230 const sum0: u128 = @as(u128, product44); 231 const sum1: u128 = @as(u128, product34) +% 232 @as(u128, product43); 233 const sum2: u128 = @as(u128, product24) +% 234 @as(u128, product33) +% 235 @as(u128, product42); 236 const sum3: u128 = @as(u128, product14) +% 237 @as(u128, product23) +% 238 @as(u128, product32) +% 239 @as(u128, product41); 240 const sum4: u128 = @as(u128, product13) +% 241 @as(u128, product22) +% 242 @as(u128, product31); 243 const sum5: u128 = @as(u128, product12) +% 244 @as(u128, product21); 245 const sum6: u128 = @as(u128, product11); 246 247 const r0: u128 = (sum0 & Word_FullMask) +% 248 ((sum1 & Word_LoMask) << 32); 249 const r1: u128 = (sum0 >> 64) +% 250 ((sum1 >> 32) & Word_FullMask) +% 251 (sum2 & Word_FullMask) +% 252 ((sum3 << 32) & Word_HiMask); 253 254 lo.* = r0 +% (r1 << 64); 255 hi.* = (r1 >> 64) +% 256 (sum1 >> 96) +% 257 (sum2 >> 64) +% 258 (sum3 >> 32) +% 259 sum4 +% 260 (sum5 << 32) +% 261 (sum6 << 64); 262 }, 263 else => @compileError("unsupported"), 264 } 265 } 266 267 fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T).Float.bits)) i32 { 268 @setRuntimeSafety(builtin.is_test); 269 const Z = std.meta.Int(.unsigned, @typeInfo(T).Float.bits); 270 const significandBits = std.math.floatMantissaBits(T); 271 const implicitBit = @as(Z, 1) << significandBits; 272 273 const shift = @clz(Z, significand.*) - @clz(Z, implicitBit); 274 significand.* <<= @intCast(std.math.Log2Int(Z), shift); 275 return @as(i32, 1) - shift; 276 } 277 278 fn wideRightShiftWithSticky(comptime Z: type, hi: *Z, lo: *Z, count: u32) void { 279 @setRuntimeSafety(builtin.is_test); 280 const typeWidth = @typeInfo(Z).Int.bits; 281 const S = std.math.Log2Int(Z); 282 if (count < typeWidth) { 283 const sticky = @boolToInt((lo.* << @intCast(S, typeWidth -% count)) != 0); 284 lo.* = (hi.* << @intCast(S, typeWidth -% count)) | (lo.* >> @intCast(S, count)) | sticky; 285 hi.* = hi.* >> @intCast(S, count); 286 } else if (count < 2 * typeWidth) { 287 const sticky = @boolToInt((hi.* << @intCast(S, 2 * typeWidth -% count) | lo.*) != 0); 288 lo.* = hi.* >> @intCast(S, count -% typeWidth) | sticky; 289 hi.* = 0; 290 } else { 291 const sticky = @boolToInt((hi.* | lo.*) != 0); 292 lo.* = sticky; 293 hi.* = 0; 294 } 295 } 296 297 test { 298 _ = @import("mulXf3_test.zig"); 299 }