1const std = @import("std");
2const math = std.math;
3const builtin = @import("builtin");
4const compiler_rt = @import("../compiler_rt.zig");
5const symbol = compiler_rt.symbol;
6
7comptime {
8 symbol(&__mulhf3, "__mulhf3");
9 if (compiler_rt.want_aeabi) {
10 symbol(&__aeabi_fmul, "__aeabi_fmul");
11 symbol(&__aeabi_dmul, "__aeabi_dmul");
12 } else {
13 symbol(&__mulsf3, "__mulsf3");
14 symbol(&__muldf3, "__muldf3");
15 }
16 symbol(&__mulxf3, "__mulxf3");
17 if (compiler_rt.want_ppc_abi) {
18 symbol(&__multf3, "__mulkf3");
19 } else if (compiler_rt.want_sparc64_abi) {
20 symbol(&_Qp_mul, "_Qp_mul");
21 } else if (compiler_rt.want_sparc32_abi) {
22 symbol(&__multf3, "_Q_mul");
23 } else {
24 symbol(&__multf3, "__multf3");
25 }
26}
27
28fn __mulhf3(a: compiler_rt.f16.Abi, b: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
29 return compiler_rt.f16.toAbi(mul_f16(compiler_rt.f16.fromAbi(a), compiler_rt.f16.fromAbi(b)));
30}
31pub fn mul_f16(a: f16, b: f16) f16 {
32 return mulf3(f16, a, b);
33}
34
35fn __mulsf3(a: compiler_rt.f32.Abi, b: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
36 return compiler_rt.f32.toAbi(mul_f32(compiler_rt.f32.fromAbi(a), compiler_rt.f32.fromAbi(b)));
37}
38fn __aeabi_fmul(a: f32, b: f32) callconv(.{ .arm_aapcs = .{} }) f32 {
39 return mul_f32(a, b);
40}
41pub fn mul_f32(a: f32, b: f32) f32 {
42 return mulf3(f32, a, b);
43}
44
45fn __muldf3(a: compiler_rt.f64.Abi, b: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.Abi {
46 return compiler_rt.f64.toAbi(mul_f64(compiler_rt.f64.fromAbi(a), compiler_rt.f64.fromAbi(b)));
47}
48fn __aeabi_dmul(a: f64, b: f64) callconv(.{ .arm_aapcs = .{} }) f64 {
49 return mul_f64(a, b);
50}
51pub fn mul_f64(a: f64, b: f64) f64 {
52 return mulf3(f64, a, b);
53}
54
55fn __mulxf3(a: compiler_rt.f80.Abi, b: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
56 return compiler_rt.f80.toAbi(mul_f80(compiler_rt.f80.fromAbi(a), compiler_rt.f80.fromAbi(b)));
57}
58pub fn mul_f80(a: f80, b: f80) f80 {
59 return mulf3(f80, a, b);
60}
61
62fn __multf3(a: compiler_rt.f128.Abi, b: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.Abi {
63 return compiler_rt.f128.toAbi(mul_f128(compiler_rt.f128.fromAbi(a), compiler_rt.f128.fromAbi(b)));
64}
65fn _Qp_mul(c: *f128, a: *const f128, b: *const f128) callconv(.c) void {
66 c.* = mul_f128(a.*, b.*);
67}
68pub fn mul_f128(a: f128, b: f128) f128 {
69 return mulf3(f128, a, b);
70}
71
72/// Ported from:
73/// https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/lib/builtins/fp_mul_impl.inc
74inline fn mulf3(comptime T: type, a: T, b: T) T {
75 @setRuntimeSafety(compiler_rt.test_safety);
76 const typeWidth = @typeInfo(T).float.bits;
77 const significandBits = math.floatMantissaBits(T);
78 const fractionalBits = math.floatFractionalBits(T);
79 const exponentBits = math.floatExponentBits(T);
80
81 const Z = @Int(.unsigned, typeWidth);
82
83 // ZSignificand is large enough to contain the significand, including an explicit integer bit
84 const ZSignificand = PowerOfTwoSignificandZ(T);
85 const ZSignificandBits = @typeInfo(ZSignificand).int.bits;
86
87 const roundBit = (1 << (ZSignificandBits - 1));
88 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
89 const maxExponent = ((1 << exponentBits) - 1);
90 const exponentBias = (maxExponent >> 1);
91
92 const integerBit = (@as(ZSignificand, 1) << fractionalBits);
93 const quietBit = integerBit >> 1;
94 const significandMask = (@as(Z, 1) << significandBits) - 1;
95
96 const absMask = signBit - 1;
97 const qnanRep = @as(Z, @bitCast(math.nan(T))) | quietBit;
98 const infRep: Z = @bitCast(math.inf(T));
99 const minNormalRep: Z = @bitCast(math.floatMin(T));
100
101 const ZExp = if (typeWidth >= 32) u32 else Z;
102 const aExponent: ZExp = @truncate((@as(Z, @bitCast(a)) >> significandBits) & maxExponent);
103 const bExponent: ZExp = @truncate((@as(Z, @bitCast(b)) >> significandBits) & maxExponent);
104 const productSign: Z = (@as(Z, @bitCast(a)) ^ @as(Z, @bitCast(b))) & signBit;
105
106 var aSignificand: ZSignificand = @intCast(@as(Z, @bitCast(a)) & significandMask);
107 var bSignificand: ZSignificand = @intCast(@as(Z, @bitCast(b)) & significandMask);
108 var scale: i32 = 0;
109
110 // Detect if a or b is zero, denormal, infinity, or NaN.
111 if (aExponent -% 1 >= maxExponent - 1 or bExponent -% 1 >= maxExponent - 1) {
112 const aAbs: Z = @as(Z, @bitCast(a)) & absMask;
113 const bAbs: Z = @as(Z, @bitCast(b)) & absMask;
114
115 // NaN * anything = qNaN
116 if (aAbs > infRep) return @bitCast(@as(Z, @bitCast(a)) | quietBit);
117 // anything * NaN = qNaN
118 if (bAbs > infRep) return @bitCast(@as(Z, @bitCast(b)) | quietBit);
119
120 if (aAbs == infRep) {
121 // infinity * non-zero = +/- infinity
122 if (bAbs != 0) {
123 return @bitCast(aAbs | productSign);
124 } else {
125 // infinity * zero = NaN
126 return @bitCast(qnanRep);
127 }
128 }
129
130 if (bAbs == infRep) {
131 //? non-zero * infinity = +/- infinity
132 if (aAbs != 0) {
133 return @bitCast(bAbs | productSign);
134 } else {
135 // zero * infinity = NaN
136 return @bitCast(qnanRep);
137 }
138 }
139
140 // zero * anything = +/- zero
141 if (aAbs == 0) return @bitCast(productSign);
142 // anything * zero = +/- zero
143 if (bAbs == 0) return @bitCast(productSign);
144
145 // one or both of a or b is denormal, the other (if applicable) is a
146 // normal number. Renormalize one or both of a and b, and set scale to
147 // include the necessary exponent adjustment.
148 if (aAbs < minNormalRep) scale += normalize(T, &aSignificand);
149 if (bAbs < minNormalRep) scale += normalize(T, &bSignificand);
150 }
151
152 // Or in the implicit significand bit. (If we fell through from the
153 // denormal path it was already set by normalize( ), but setting it twice
154 // won't hurt anything.)
155 aSignificand |= integerBit;
156 bSignificand |= integerBit;
157
158 // Get the significand of a*b. Before multiplying the significands, shift
159 // one of them left to left-align it in the field. Thus, the product will
160 // have (exponentBits + 2) integral digits, all but two of which must be
161 // zero. Normalizing this result is just a conditional left-shift by one
162 // and bumping the exponent accordingly.
163 var productHi: ZSignificand = undefined;
164 var productLo: ZSignificand = undefined;
165 const left_align_shift = ZSignificandBits - fractionalBits - 1;
166 compiler_rt.wideMultiply(ZSignificand, aSignificand, bSignificand << left_align_shift, &productHi, &productLo);
167
168 var productExponent: i32 = @as(i32, @intCast(aExponent + bExponent)) - exponentBias + scale;
169
170 // Normalize the significand, adjust exponent if needed.
171 if ((productHi & integerBit) != 0) {
172 productExponent +%= 1;
173 } else {
174 productHi = (productHi << 1) | (productLo >> (ZSignificandBits - 1));
175 productLo = productLo << 1;
176 }
177
178 // If we have overflowed the type, return +/- infinity.
179 if (productExponent >= maxExponent) return @bitCast(infRep | productSign);
180
181 var result: Z = undefined;
182 if (productExponent <= 0) {
183 // Result is denormal before rounding
184 //
185 // If the result is so small that it just underflows to zero, return
186 // a zero of the appropriate sign. Mathematically there is no need to
187 // handle this case separately, but we make it a special case to
188 // simplify the shift logic.
189 const shift: u32 = @truncate(@as(Z, 1) -% @as(u32, @bitCast(productExponent)));
190 if (shift >= ZSignificandBits) return @bitCast(productSign);
191
192 // Otherwise, shift the significand of the result so that the round
193 // bit is the high bit of productLo.
194 const sticky = wideShrWithTruncation(ZSignificand, &productHi, &productLo, shift);
195 productLo |= @intFromBool(sticky);
196 result = productHi;
197
198 // We include the integer bit so that rounding will carry to the exponent,
199 // but it will be removed later if the result is still denormal
200 if (significandBits != fractionalBits) result |= integerBit;
201 } else {
202 // Result is normal before rounding; insert the exponent.
203 result = productHi & significandMask;
204 result |= @as(Z, @intCast(productExponent)) << significandBits;
205 }
206
207 // Final rounding. The final result may overflow to infinity, or underflow
208 // to zero, but those are the correct results in those cases. We use the
209 // default IEEE-754 round-to-nearest, ties-to-even rounding mode.
210 if (productLo > roundBit) result +%= 1;
211 if (productLo == roundBit) result +%= result & 1;
212
213 // Restore any explicit integer bit, if it was rounded off
214 if (significandBits != fractionalBits) {
215 if ((result >> significandBits) != 0) {
216 result |= integerBit;
217 } else {
218 result &= ~integerBit;
219 }
220 }
221
222 // Insert the sign of the result:
223 result |= productSign;
224
225 return @bitCast(result);
226}
227
228/// Returns `true` if the right shift is inexact (i.e. any bit shifted out is non-zero)
229///
230/// This is analogous to an shr version of `@shlWithOverflow`
231fn wideShrWithTruncation(comptime Z: type, hi: *Z, lo: *Z, count: u32) bool {
232 @setRuntimeSafety(compiler_rt.test_safety);
233 const typeWidth = @typeInfo(Z).int.bits;
234 var inexact = false;
235 if (count < typeWidth) {
236 inexact = (lo.* << @intCast(typeWidth -% count)) != 0;
237 lo.* = (hi.* << @intCast(typeWidth -% count)) | (lo.* >> @intCast(count));
238 hi.* = hi.* >> @intCast(count);
239 } else if (count < 2 * typeWidth) {
240 inexact = (hi.* << @intCast(2 * typeWidth -% count) | lo.*) != 0;
241 lo.* = hi.* >> @intCast(count -% typeWidth);
242 hi.* = 0;
243 } else {
244 inexact = (hi.* | lo.*) != 0;
245 lo.* = 0;
246 hi.* = 0;
247 }
248 return inexact;
249}
250
251fn normalize(comptime T: type, significand: *PowerOfTwoSignificandZ(T)) i32 {
252 const Z = PowerOfTwoSignificandZ(T);
253 const integerBit = @as(Z, 1) << math.floatFractionalBits(T);
254
255 const shift = @clz(significand.*) - @clz(integerBit);
256 significand.* <<= @intCast(shift);
257 return @as(i32, 1) - shift;
258}
259
260/// Returns a power-of-two integer type that is large enough to contain
261/// the significand of T, including an explicit integer bit
262fn PowerOfTwoSignificandZ(comptime T: type) type {
263 const bits = math.ceilPowerOfTwoAssert(u16, math.floatFractionalBits(T) + 1);
264 return @Int(.unsigned, bits);
265}
266
267test {
268 _ = @import("mulf3_test.zig");
269}