1//! Ported from musl, which is MIT licensed:
2//! https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//!
4//! https://git.musl-libc.org/cgit/musl/tree/src/math/fmal.c
5//! https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c
6//! https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c
7
8const std = @import("std");
9const math = std.math;
10const expect = std.testing.expect;
11const compiler_rt = @import("../compiler_rt.zig");
12const symbol = compiler_rt.symbol;
13
14comptime {
15 symbol(&__fmah, "__fmah");
16 symbol(&fmaf, "fmaf");
17 symbol(&fma, "fma");
18 symbol(&__fmax, "__fmax");
19 symbol(&fmaq, "fmaf128");
20 symbol(&fmal, "fmal");
21}
22
23fn __fmah(x: compiler_rt.f16.Abi, y: compiler_rt.f16.Abi, z: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
24 return compiler_rt.f16.toAbi(fma_f16(compiler_rt.f16.fromAbi(x), compiler_rt.f16.fromAbi(y), compiler_rt.f16.fromAbi(z)));
25}
26pub fn fma_f16(x: f16, y: f16, z: f16) f16 {
27 // TODO: more efficient implementation
28 return @floatCast(fma_f32(x, y, z));
29}
30
31fn fmaf(x: compiler_rt.f32.Abi, y: compiler_rt.f32.Abi, z: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
32 return compiler_rt.f32.toAbi(fma_f32(compiler_rt.f32.fromAbi(x), compiler_rt.f32.fromAbi(y), compiler_rt.f32.fromAbi(z)));
33}
34pub fn fma_f32(x: f32, y: f32, z: f32) f32 {
35 const xy = @as(f64, x) * y;
36 const xy_z = xy + z;
37 const u = @as(u64, @bitCast(xy_z));
38 const e = (u >> 52) & 0x7FF;
39
40 if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or (xy_z - xy == z and xy_z - z == xy)) {
41 return @floatCast(xy_z);
42 } else {
43 // TODO: Handle inexact case with double-rounding
44 return @floatCast(xy_z);
45 }
46}
47
48fn fma(x: compiler_rt.f64.Abi, y: compiler_rt.f64.Abi, z: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.Abi {
49 return compiler_rt.f64.toAbi(fma_f64(compiler_rt.f64.fromAbi(x), compiler_rt.f64.fromAbi(y), compiler_rt.f64.fromAbi(z)));
50}
51/// NOTE: Upstream fma.c has been rewritten completely to raise fp exceptions more accurately.
52pub fn fma_f64(x: f64, y: f64, z: f64) f64 {
53 if (!math.isFinite(x) or !math.isFinite(y)) {
54 return x * y + z;
55 }
56 if (!math.isFinite(z)) {
57 return z;
58 }
59 if (x == 0.0 or y == 0.0) {
60 return x * y + z;
61 }
62 if (z == 0.0) {
63 return x * y;
64 }
65
66 const x1 = math.frexp(x);
67 const ex = x1.exponent;
68 const xs = x1.significand;
69 const x2 = math.frexp(y);
70 const ey = x2.exponent;
71 const ys = x2.significand;
72 const x3 = math.frexp(z);
73 const ez = x3.exponent;
74 var zs = x3.significand;
75
76 var spread = ex + ey - ez;
77 if (spread <= 53 * 2) {
78 zs = math.scalbn(zs, -spread);
79 } else {
80 zs = math.copysign(math.floatMin(f64), zs);
81 }
82
83 const xy = dd_mul(xs, ys);
84 const r = dd_add(xy.hi, zs);
85 spread = ex + ey;
86
87 if (r.hi == 0.0) {
88 return xy.hi + zs + math.scalbn(xy.lo, spread);
89 }
90
91 const adj = add_adjusted(r.lo, xy.lo);
92 if (spread + math.ilogb(r.hi) > -1023) {
93 return math.scalbn(r.hi + adj, spread);
94 } else {
95 return add_and_denorm(r.hi, adj, spread);
96 }
97}
98
99fn __fmax(a: compiler_rt.f80.Abi, b: compiler_rt.f80.Abi, c: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
100 return compiler_rt.f80.toAbi(fma_f80(compiler_rt.f80.fromAbi(a), compiler_rt.f80.fromAbi(b), compiler_rt.f80.fromAbi(c)));
101}
102pub fn fma_f80(a: f80, b: f80, c: f80) f80 {
103 // TODO: more efficient implementation
104 return @floatCast(fma_f128(a, b, c));
105}
106
107fn fmaq(x: compiler_rt.f128.Abi, y: compiler_rt.f128.Abi, z: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.Abi {
108 return compiler_rt.f128.toAbi(fma_f128(compiler_rt.f128.fromAbi(x), compiler_rt.f128.fromAbi(y), compiler_rt.f128.fromAbi(z)));
109}
110/// Fused multiply-add: Compute x * y + z with a single rounding error.
111///
112/// We use scaling to avoid overflow/underflow, along with the
113/// canonical precision-doubling technique adapted from:
114///
115/// Dekker, T. A Floating-Point Technique for Extending the
116/// Available Precision. Numer. Math. 18, 224-242 (1971).
117pub fn fma_f128(x: f128, y: f128, z: f128) f128 {
118 if (!math.isFinite(x) or !math.isFinite(y)) {
119 return x * y + z;
120 }
121 if (!math.isFinite(z)) {
122 return z;
123 }
124 if (x == 0.0 or y == 0.0) {
125 return x * y + z;
126 }
127 if (z == 0.0) {
128 return x * y;
129 }
130
131 const x1 = math.frexp(x);
132 const ex = x1.exponent;
133 const xs = x1.significand;
134 const x2 = math.frexp(y);
135 const ey = x2.exponent;
136 const ys = x2.significand;
137 const x3 = math.frexp(z);
138 const ez = x3.exponent;
139 var zs = x3.significand;
140
141 var spread = ex + ey - ez;
142 if (spread <= 113 * 2) {
143 zs = math.scalbn(zs, -spread);
144 } else {
145 zs = math.copysign(math.floatMin(f128), zs);
146 }
147
148 const xy = dd_mul128(xs, ys);
149 const r = dd_add128(xy.hi, zs);
150 spread = ex + ey;
151
152 if (r.hi == 0.0) {
153 return xy.hi + zs + math.scalbn(xy.lo, spread);
154 }
155
156 const adj = add_adjusted128(r.lo, xy.lo);
157 if (spread + math.ilogb(r.hi) > -16383) {
158 return math.scalbn(r.hi + adj, spread);
159 } else {
160 return add_and_denorm128(r.hi, adj, spread);
161 }
162}
163
164pub fn fmal(x: c_longdouble, y: c_longdouble, z: c_longdouble) callconv(.c) c_longdouble {
165 switch (@typeInfo(c_longdouble).float.bits) {
166 64 => return fma_f64(x, y, z),
167 80 => return fma_f80(x, y, z),
168 128 => return fma_f128(x, y, z),
169 else => comptime unreachable,
170 }
171}
172
173const dd = struct {
174 hi: f64,
175 lo: f64,
176};
177
178fn dd_add(a: f64, b: f64) dd {
179 var ret: dd = undefined;
180 ret.hi = a + b;
181 const s = ret.hi - a;
182 ret.lo = (a - (ret.hi - s)) + (b - s);
183 return ret;
184}
185
186fn dd_mul(a: f64, b: f64) dd {
187 var ret: dd = undefined;
188 const split: f64 = 0x1.0p27 + 1.0;
189
190 var p = a * split;
191 var ha = a - p;
192 ha += p;
193 const la = a - ha;
194
195 p = b * split;
196 var hb = b - p;
197 hb += p;
198 const lb = b - hb;
199
200 p = ha * hb;
201 const q = ha * lb + la * hb;
202
203 ret.hi = p + q;
204 ret.lo = p - ret.hi + q + la * lb;
205 return ret;
206}
207
208fn add_adjusted(a: f64, b: f64) f64 {
209 var sum = dd_add(a, b);
210 if (sum.lo != 0) {
211 var uhii: u64 = @bitCast(sum.hi);
212 if (uhii & 1 == 0) {
213 // hibits += copysign(1.0, sum.hi, sum.lo)
214 const uloi: u64 = @bitCast(sum.lo);
215 uhii = uhii + 1 - ((uhii ^ uloi) >> 62);
216 sum.hi = @bitCast(uhii);
217 }
218 }
219 return sum.hi;
220}
221
222fn add_and_denorm(a: f64, b: f64, scale: i32) f64 {
223 var sum = dd_add(a, b);
224 if (sum.lo != 0) {
225 var uhii: u64 = @bitCast(sum.hi);
226 const bits_lost = -@as(i32, @intCast((uhii >> 52) & 0x7FF)) - scale + 1;
227 if ((bits_lost != 1) == (uhii & 1 != 0)) {
228 const uloi: u64 = @bitCast(sum.lo);
229 uhii = uhii + 1 - (((uhii ^ uloi) >> 62) & 2);
230 sum.hi = @bitCast(uhii);
231 }
232 }
233 return math.scalbn(sum.hi, scale);
234}
235
236/// A struct that represents a floating-point number with twice the precision
237/// of f128. We maintain the invariant that "hi" stores the high-order
238/// bits of the result.
239const dd128 = struct {
240 hi: f128,
241 lo: f128,
242};
243
244/// Compute a+b exactly, returning the exact result in a struct dd. We assume
245/// that both a and b are finite, but make no assumptions about their relative
246/// magnitudes.
247fn dd_add128(a: f128, b: f128) dd128 {
248 var ret: dd128 = undefined;
249 ret.hi = a + b;
250 const s = ret.hi - a;
251 ret.lo = (a - (ret.hi - s)) + (b - s);
252 return ret;
253}
254
255/// Compute a+b, with a small tweak: The least significant bit of the
256/// result is adjusted into a sticky bit summarizing all the bits that
257/// were lost to rounding. This adjustment negates the effects of double
258/// rounding when the result is added to another number with a higher
259/// exponent. For an explanation of round and sticky bits, see any reference
260/// on FPU design, e.g.,
261///
262/// J. Coonen. An Implementation Guide to a Proposed Standard for
263/// Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980.
264fn add_adjusted128(a: f128, b: f128) f128 {
265 var sum = dd_add128(a, b);
266 if (sum.lo != 0) {
267 var uhii: u128 = @bitCast(sum.hi);
268 if (uhii & 1 == 0) {
269 // hibits += copysign(1.0, sum.hi, sum.lo)
270 const uloi: u128 = @bitCast(sum.lo);
271 uhii = uhii + 1 - ((uhii ^ uloi) >> 126);
272 sum.hi = @bitCast(uhii);
273 }
274 }
275 return sum.hi;
276}
277
278/// Compute ldexp(a+b, scale) with a single rounding error. It is assumed
279/// that the result will be subnormal, and care is taken to ensure that
280/// double rounding does not occur.
281fn add_and_denorm128(a: f128, b: f128, scale: i32) f128 {
282 var sum = dd_add128(a, b);
283 // If we are losing at least two bits of accuracy to denormalization,
284 // then the first lost bit becomes a round bit, and we adjust the
285 // lowest bit of sum.hi to make it a sticky bit summarizing all the
286 // bits in sum.lo. With the sticky bit adjusted, the hardware will
287 // break any ties in the correct direction.
288 //
289 // If we are losing only one bit to denormalization, however, we must
290 // break the ties manually.
291 if (sum.lo != 0) {
292 var uhii: u128 = @bitCast(sum.hi);
293 const bits_lost = -@as(i32, @intCast((uhii >> 112) & 0x7FFF)) - scale + 1;
294 if ((bits_lost != 1) == (uhii & 1 != 0)) {
295 const uloi: u128 = @bitCast(sum.lo);
296 uhii = uhii + 1 - (((uhii ^ uloi) >> 126) & 2);
297 sum.hi = @bitCast(uhii);
298 }
299 }
300 return math.scalbn(sum.hi, scale);
301}
302
303/// Compute a*b exactly, returning the exact result in a struct dd. We assume
304/// that both a and b are normalized, so no underflow or overflow will occur.
305/// The current rounding mode must be round-to-nearest.
306fn dd_mul128(a: f128, b: f128) dd128 {
307 var ret: dd128 = undefined;
308 const split: f128 = 0x1.0p57 + 1.0;
309
310 var p = a * split;
311 var ha = a - p;
312 ha += p;
313 const la = a - ha;
314
315 p = b * split;
316 var hb = b - p;
317 hb += p;
318 const lb = b - hb;
319
320 p = ha * hb;
321 const q = ha * lb + la * hb;
322
323 ret.hi = p + q;
324 ret.lo = p - ret.hi + q + la * lb;
325 return ret;
326}
327
328test "32" {
329 const epsilon = 0.000001;
330
331 try expect(math.approxEqAbs(f32, fma_f32(0.0, 5.0, 9.124), 9.124, epsilon));
332 try expect(math.approxEqAbs(f32, fma_f32(0.2, 5.0, 9.124), 10.124, epsilon));
333 try expect(math.approxEqAbs(f32, fma_f32(0.8923, 5.0, 9.124), 13.5855, epsilon));
334 try expect(math.approxEqAbs(f32, fma_f32(1.5, 5.0, 9.124), 16.624, epsilon));
335 try expect(math.approxEqAbs(f32, fma_f32(37.45, 5.0, 9.124), 196.374004, epsilon));
336 try expect(math.approxEqAbs(f32, fma_f32(89.123, 5.0, 9.124), 454.739005, epsilon));
337 try expect(math.approxEqAbs(f32, fma_f32(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
338}
339
340test "64" {
341 const epsilon = 0.000001;
342
343 try expect(math.approxEqAbs(f64, fma_f64(0.0, 5.0, 9.124), 9.124, epsilon));
344 try expect(math.approxEqAbs(f64, fma_f64(0.2, 5.0, 9.124), 10.124, epsilon));
345 try expect(math.approxEqAbs(f64, fma_f64(0.8923, 5.0, 9.124), 13.5855, epsilon));
346 try expect(math.approxEqAbs(f64, fma_f64(1.5, 5.0, 9.124), 16.624, epsilon));
347 try expect(math.approxEqAbs(f64, fma_f64(37.45, 5.0, 9.124), 196.374, epsilon));
348 try expect(math.approxEqAbs(f64, fma_f64(89.123, 5.0, 9.124), 454.739, epsilon));
349 try expect(math.approxEqAbs(f64, fma_f64(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
350}
351
352test "128" {
353 const epsilon = 0.000001;
354
355 try expect(math.approxEqAbs(f128, fma_f128(0.0, 5.0, 9.124), 9.124, epsilon));
356 try expect(math.approxEqAbs(f128, fma_f128(0.2, 5.0, 9.124), 10.124, epsilon));
357 try expect(math.approxEqAbs(f128, fma_f128(0.8923, 5.0, 9.124), 13.5855, epsilon));
358 try expect(math.approxEqAbs(f128, fma_f128(1.5, 5.0, 9.124), 16.624, epsilon));
359 try expect(math.approxEqAbs(f128, fma_f128(37.45, 5.0, 9.124), 196.374, epsilon));
360 try expect(math.approxEqAbs(f128, fma_f128(89.123, 5.0, 9.124), 454.739, epsilon));
361 try expect(math.approxEqAbs(f128, fma_f128(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
362}