authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2021-10-05 12:32:26-07:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2021-10-05 12:32:26-07:00
log2fe7b06f3df06b7a442ffc3d9b951c5d52a11a24
treeb94c0f9ede8dd1bd01148f9e70e2f7359a15b040
parent6115cf22404467fd13d0290fc022d51d372d139a

add support for f128 `@mulAdd`

std: add f128 implementations of fma, frexp, and ilogb. Expose `fmal` in zig's freestanding libc. This makes `@mulAdd` work correctly for f128. Fixes a CI regression from yesterday, where I added a usage of f128 `@mulAdd` into the self-hosted compiler.

6 files changed, 324 insertions(+), 37 deletions(-)

lib/std/math.zig+1-2
...@@ -229,8 +229,7 @@ pub const floor = @import("math/floor.zig").floor;...@@ -229,8 +229,7 @@ pub const floor = @import("math/floor.zig").floor;
229pub const trunc = @import("math/trunc.zig").trunc;229pub const trunc = @import("math/trunc.zig").trunc;
230pub const round = @import("math/round.zig").round;230pub const round = @import("math/round.zig").round;
231pub const frexp = @import("math/frexp.zig").frexp;231pub const frexp = @import("math/frexp.zig").frexp;
232pub const frexp32_result = @import("math/frexp.zig").frexp32_result;232pub const Frexp = @import("math/frexp.zig").Frexp;
233pub const frexp64_result = @import("math/frexp.zig").frexp64_result;
234pub const modf = @import("math/modf.zig").modf;233pub const modf = @import("math/modf.zig").modf;
235pub const modf32_result = @import("math/modf.zig").modf32_result;234pub const modf32_result = @import("math/modf.zig").modf32_result;
236pub const modf64_result = @import("math/modf.zig").modf64_result;235pub const modf64_result = @import("math/modf.zig").modf64_result;
lib/std/math/fma.zig+165-4
...@@ -1,6 +1,7 @@...@@ -1,6 +1,7 @@
1// Ported from musl, which is licensed under the MIT license:1// Ported from musl, which is MIT licensed:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/fmal.c
4// https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c5// https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c6// https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c
67
...@@ -13,6 +14,7 @@ pub fn fma(comptime T: type, x: T, y: T, z: T) T {...@@ -13,6 +14,7 @@ pub fn fma(comptime T: type, x: T, y: T, z: T) T {
13 return switch (T) {14 return switch (T) {
14 f32 => fma32(x, y, z),15 f32 => fma32(x, y, z),
15 f64 => fma64(x, y, z),16 f64 => fma64(x, y, z),
17 f128 => fma128(x, y, z),
16 else => @compileError("fma not implemented for " ++ @typeName(T)),18 else => @compileError("fma not implemented for " ++ @typeName(T)),
17 };19 };
18}20}
...@@ -142,12 +144,159 @@ fn add_and_denorm(a: f64, b: f64, scale: i32) f64 {...@@ -142,12 +144,159 @@ fn add_and_denorm(a: f64, b: f64, scale: i32) f64 {
142 return math.scalbn(sum.hi, scale);144 return math.scalbn(sum.hi, scale);
143}145}
144146
145test "math.fma" {147/// A struct that represents a floating-point number with twice the precision
148/// of f128. We maintain the invariant that "hi" stores the high-order
149/// bits of the result.
150const dd128 = struct {
151 hi: f128,
152 lo: f128,
153};
154
155/// Compute a+b exactly, returning the exact result in a struct dd. We assume
156/// that both a and b are finite, but make no assumptions about their relative
157/// magnitudes.
158fn dd_add128(a: f128, b: f128) dd128 {
159 var ret: dd128 = undefined;
160 ret.hi = a + b;
161 const s = ret.hi - a;
162 ret.lo = (a - (ret.hi - s)) + (b - s);
163 return ret;
164}
165
166/// Compute a+b, with a small tweak: The least significant bit of the
167/// result is adjusted into a sticky bit summarizing all the bits that
168/// were lost to rounding. This adjustment negates the effects of double
169/// rounding when the result is added to another number with a higher
170/// exponent. For an explanation of round and sticky bits, see any reference
171/// on FPU design, e.g.,
172///
173/// J. Coonen. An Implementation Guide to a Proposed Standard for
174/// Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980.
175fn add_adjusted128(a: f128, b: f128) f128 {
176 var sum = dd_add128(a, b);
177 if (sum.lo != 0) {
178 var uhii = @bitCast(u128, sum.hi);
179 if (uhii & 1 == 0) {
180 // hibits += copysign(1.0, sum.hi, sum.lo)
181 const uloi = @bitCast(u128, sum.lo);
182 uhii += 1 - ((uhii ^ uloi) >> 126);
183 sum.hi = @bitCast(f128, uhii);
184 }
185 }
186 return sum.hi;
187}
188
189/// Compute ldexp(a+b, scale) with a single rounding error. It is assumed
190/// that the result will be subnormal, and care is taken to ensure that
191/// double rounding does not occur.
192fn add_and_denorm128(a: f128, b: f128, scale: i32) f128 {
193 var sum = dd_add128(a, b);
194 // If we are losing at least two bits of accuracy to denormalization,
195 // then the first lost bit becomes a round bit, and we adjust the
196 // lowest bit of sum.hi to make it a sticky bit summarizing all the
197 // bits in sum.lo. With the sticky bit adjusted, the hardware will
198 // break any ties in the correct direction.
199 //
200 // If we are losing only one bit to denormalization, however, we must
201 // break the ties manually.
202 if (sum.lo != 0) {
203 var uhii = @bitCast(u128, sum.hi);
204 const bits_lost = -@intCast(i32, (uhii >> 112) & 0x7FFF) - scale + 1;
205 if ((bits_lost != 1) == (uhii & 1 != 0)) {
206 const uloi = @bitCast(u128, sum.lo);
207 uhii += 1 - (((uhii ^ uloi) >> 126) & 2);
208 sum.hi = @bitCast(f128, uhii);
209 }
210 }
211 return math.scalbn(sum.hi, scale);
212}
213
214/// Compute a*b exactly, returning the exact result in a struct dd. We assume
215/// that both a and b are normalized, so no underflow or overflow will occur.
216/// The current rounding mode must be round-to-nearest.
217fn dd_mul128(a: f128, b: f128) dd128 {
218 var ret: dd128 = undefined;
219 const split: f128 = 0x1.0p57 + 1.0;
220
221 var p = a * split;
222 var ha = a - p;
223 ha += p;
224 var la = a - ha;
225
226 p = b * split;
227 var hb = b - p;
228 hb += p;
229 var lb = b - hb;
230
231 p = ha * hb;
232 var q = ha * lb + la * hb;
233
234 ret.hi = p + q;
235 ret.lo = p - ret.hi + q + la * lb;
236 return ret;
237}
238
239/// Fused multiply-add: Compute x * y + z with a single rounding error.
240///
241/// We use scaling to avoid overflow/underflow, along with the
242/// canonical precision-doubling technique adapted from:
243///
244/// Dekker, T. A Floating-Point Technique for Extending the
245/// Available Precision. Numer. Math. 18, 224-242 (1971).
246fn fma128(x: f128, y: f128, z: f128) f128 {
247 if (!math.isFinite(x) or !math.isFinite(y)) {
248 return x * y + z;
249 }
250 if (!math.isFinite(z)) {
251 return z;
252 }
253 if (x == 0.0 or y == 0.0) {
254 return x * y + z;
255 }
256 if (z == 0.0) {
257 return x * y;
258 }
259
260 const x1 = math.frexp(x);
261 var ex = x1.exponent;
262 var xs = x1.significand;
263 const x2 = math.frexp(y);
264 var ey = x2.exponent;
265 var ys = x2.significand;
266 const x3 = math.frexp(z);
267 var ez = x3.exponent;
268 var zs = x3.significand;
269
270 var spread = ex + ey - ez;
271 if (spread <= 113 * 2) {
272 zs = math.scalbn(zs, -spread);
273 } else {
274 zs = math.copysign(f128, math.f128_min, zs);
275 }
276
277 const xy = dd_mul128(xs, ys);
278 const r = dd_add128(xy.hi, zs);
279 spread = ex + ey;
280
281 if (r.hi == 0.0) {
282 return xy.hi + zs + math.scalbn(xy.lo, spread);
283 }
284
285 const adj = add_adjusted128(r.lo, xy.lo);
286 if (spread + math.ilogb(r.hi) > -16383) {
287 return math.scalbn(r.hi + adj, spread);
288 } else {
289 return add_and_denorm128(r.hi, adj, spread);
290 }
291}
292
293test "type dispatch" {
146 try expect(fma(f32, 0.0, 1.0, 1.0) == fma32(0.0, 1.0, 1.0));294 try expect(fma(f32, 0.0, 1.0, 1.0) == fma32(0.0, 1.0, 1.0));
147 try expect(fma(f64, 0.0, 1.0, 1.0) == fma64(0.0, 1.0, 1.0));295 try expect(fma(f64, 0.0, 1.0, 1.0) == fma64(0.0, 1.0, 1.0));
296 try expect(fma(f128, 0.0, 1.0, 1.0) == fma128(0.0, 1.0, 1.0));
148}297}
149298
150test "math.fma32" {299test "32" {
151 const epsilon = 0.000001;300 const epsilon = 0.000001;
152301
153 try expect(math.approxEqAbs(f32, fma32(0.0, 5.0, 9.124), 9.124, epsilon));302 try expect(math.approxEqAbs(f32, fma32(0.0, 5.0, 9.124), 9.124, epsilon));
...@@ -159,7 +308,7 @@ test "math.fma32" {...@@ -159,7 +308,7 @@ test "math.fma32" {
159 try expect(math.approxEqAbs(f32, fma32(123123.234375, 5.0, 9.124), 615625.295875, epsilon));308 try expect(math.approxEqAbs(f32, fma32(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
160}309}
161310
162test "math.fma64" {311test "64" {
163 const epsilon = 0.000001;312 const epsilon = 0.000001;
164313
165 try expect(math.approxEqAbs(f64, fma64(0.0, 5.0, 9.124), 9.124, epsilon));314 try expect(math.approxEqAbs(f64, fma64(0.0, 5.0, 9.124), 9.124, epsilon));
...@@ -170,3 +319,15 @@ test "math.fma64" {...@@ -170,3 +319,15 @@ test "math.fma64" {
170 try expect(math.approxEqAbs(f64, fma64(89.123, 5.0, 9.124), 454.739, epsilon));319 try expect(math.approxEqAbs(f64, fma64(89.123, 5.0, 9.124), 454.739, epsilon));
171 try expect(math.approxEqAbs(f64, fma64(123123.234375, 5.0, 9.124), 615625.295875, epsilon));320 try expect(math.approxEqAbs(f64, fma64(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
172}321}
322
323test "128" {
324 const epsilon = 0.000001;
325
326 try expect(math.approxEqAbs(f128, fma128(0.0, 5.0, 9.124), 9.124, epsilon));
327 try expect(math.approxEqAbs(f128, fma128(0.2, 5.0, 9.124), 10.124, epsilon));
328 try expect(math.approxEqAbs(f128, fma128(0.8923, 5.0, 9.124), 13.5855, epsilon));
329 try expect(math.approxEqAbs(f128, fma128(1.5, 5.0, 9.124), 16.624, epsilon));
330 try expect(math.approxEqAbs(f128, fma128(37.45, 5.0, 9.124), 196.374, epsilon));
331 try expect(math.approxEqAbs(f128, fma128(89.123, 5.0, 9.124), 454.739, epsilon));
332 try expect(math.approxEqAbs(f128, fma128(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
333}
lib/std/math/frexp.zig+91-18
...@@ -1,6 +1,7 @@...@@ -1,6 +1,7 @@
1// Ported from musl, which is licensed under the MIT license:1// Ported from musl, which is MIT licensed:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpl.c
4// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpf.c5// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/frexp.c6// https://git.musl-libc.org/cgit/musl/tree/src/math/frexp.c
67
...@@ -8,14 +9,12 @@ const std = @import("../std.zig");...@@ -8,14 +9,12 @@ const std = @import("../std.zig");
8const math = std.math;9const math = std.math;
9const expect = std.testing.expect;10const expect = std.testing.expect;
1011
11fn frexp_result(comptime T: type) type {12pub fn Frexp(comptime T: type) type {
12 return struct {13 return struct {
13 significand: T,14 significand: T,
14 exponent: i32,15 exponent: i32,
15 };16 };
16}17}
17pub const frexp32_result = frexp_result(f32);
18pub const frexp64_result = frexp_result(f64);
1918
20/// Breaks x into a normalized fraction and an integral power of two.19/// Breaks x into a normalized fraction and an integral power of two.
21/// f == frac * 2^exp, with |frac| in the interval [0.5, 1).20/// f == frac * 2^exp, with |frac| in the interval [0.5, 1).
...@@ -24,17 +23,20 @@ pub const frexp64_result = frexp_result(f64);...@@ -24,17 +23,20 @@ pub const frexp64_result = frexp_result(f64);
24/// - frexp(+-0) = +-0, 023/// - frexp(+-0) = +-0, 0
25/// - frexp(+-inf) = +-inf, 024/// - frexp(+-inf) = +-inf, 0
26/// - frexp(nan) = nan, undefined25/// - frexp(nan) = nan, undefined
27pub fn frexp(x: anytype) frexp_result(@TypeOf(x)) {26pub fn frexp(x: anytype) Frexp(@TypeOf(x)) {
28 const T = @TypeOf(x);27 const T = @TypeOf(x);
29 return switch (T) {28 return switch (T) {
30 f32 => frexp32(x),29 f32 => frexp32(x),
31 f64 => frexp64(x),30 f64 => frexp64(x),
31 f128 => frexp128(x),
32 else => @compileError("frexp not implemented for " ++ @typeName(T)),32 else => @compileError("frexp not implemented for " ++ @typeName(T)),
33 };33 };
34}34}
3535
36fn frexp32(x: f32) frexp32_result {36// TODO: unify all these implementations using generics
37 var result: frexp32_result = undefined;37
38fn frexp32(x: f32) Frexp(f32) {
39 var result: Frexp(f32) = undefined;
3840
39 var y = @bitCast(u32, x);41 var y = @bitCast(u32, x);
40 const e = @intCast(i32, y >> 23) & 0xFF;42 const e = @intCast(i32, y >> 23) & 0xFF;
...@@ -70,8 +72,8 @@ fn frexp32(x: f32) frexp32_result {...@@ -70,8 +72,8 @@ fn frexp32(x: f32) frexp32_result {
70 return result;72 return result;
71}73}
7274
73fn frexp64(x: f64) frexp64_result {75fn frexp64(x: f64) Frexp(f64) {
74 var result: frexp64_result = undefined;76 var result: Frexp(f64) = undefined;
7577
76 var y = @bitCast(u64, x);78 var y = @bitCast(u64, x);
77 const e = @intCast(i32, y >> 52) & 0x7FF;79 const e = @intCast(i32, y >> 52) & 0x7FF;
...@@ -107,7 +109,44 @@ fn frexp64(x: f64) frexp64_result {...@@ -107,7 +109,44 @@ fn frexp64(x: f64) frexp64_result {
107 return result;109 return result;
108}110}
109111
110test "math.frexp" {112fn frexp128(x: f128) Frexp(f128) {
113 var result: Frexp(f128) = undefined;
114
115 var y = @bitCast(u128, x);
116 const e = @intCast(i32, y >> 112) & 0x7FFF;
117
118 if (e == 0) {
119 if (x != 0) {
120 // subnormal
121 result = frexp128(x * 0x1.0p120);
122 result.exponent -= 120;
123 } else {
124 // frexp(+-0) = (+-0, 0)
125 result.significand = x;
126 result.exponent = 0;
127 }
128 return result;
129 } else if (e == 0x7FFF) {
130 // frexp(nan) = (nan, undefined)
131 result.significand = x;
132 result.exponent = undefined;
133
134 // frexp(+-inf) = (+-inf, 0)
135 if (math.isInf(x)) {
136 result.exponent = 0;
137 }
138
139 return result;
140 }
141
142 result.exponent = e - 0x3FFE;
143 y &= 0x8000FFFFFFFFFFFFFFFFFFFFFFFFFFFF;
144 y |= 0x3FFE0000000000000000000000000000;
145 result.significand = @bitCast(f128, y);
146 return result;
147}
148
149test "type dispatch" {
111 const a = frexp(@as(f32, 1.3));150 const a = frexp(@as(f32, 1.3));
112 const b = frexp32(1.3);151 const b = frexp32(1.3);
113 try expect(a.significand == b.significand and a.exponent == b.exponent);152 try expect(a.significand == b.significand and a.exponent == b.exponent);
...@@ -115,11 +154,15 @@ test "math.frexp" {...@@ -115,11 +154,15 @@ test "math.frexp" {
115 const c = frexp(@as(f64, 1.3));154 const c = frexp(@as(f64, 1.3));
116 const d = frexp64(1.3);155 const d = frexp64(1.3);
117 try expect(c.significand == d.significand and c.exponent == d.exponent);156 try expect(c.significand == d.significand and c.exponent == d.exponent);
157
158 const e = frexp(@as(f128, 1.3));
159 const f = frexp128(1.3);
160 try expect(e.significand == f.significand and e.exponent == f.exponent);
118}161}
119162
120test "math.frexp32" {163test "32" {
121 const epsilon = 0.000001;164 const epsilon = 0.000001;
122 var r: frexp32_result = undefined;165 var r: Frexp(f32) = undefined;
123166
124 r = frexp32(1.3);167 r = frexp32(1.3);
125 try expect(math.approxEqAbs(f32, r.significand, 0.65, epsilon) and r.exponent == 1);168 try expect(math.approxEqAbs(f32, r.significand, 0.65, epsilon) and r.exponent == 1);
...@@ -128,9 +171,9 @@ test "math.frexp32" {...@@ -128,9 +171,9 @@ test "math.frexp32" {
128 try expect(math.approxEqAbs(f32, r.significand, 0.609558, epsilon) and r.exponent == 7);171 try expect(math.approxEqAbs(f32, r.significand, 0.609558, epsilon) and r.exponent == 7);
129}172}
130173
131test "math.frexp64" {174test "64" {
132 const epsilon = 0.000001;175 const epsilon = 0.000001;
133 var r: frexp64_result = undefined;176 var r: Frexp(f64) = undefined;
134177
135 r = frexp64(1.3);178 r = frexp64(1.3);
136 try expect(math.approxEqAbs(f64, r.significand, 0.65, epsilon) and r.exponent == 1);179 try expect(math.approxEqAbs(f64, r.significand, 0.65, epsilon) and r.exponent == 1);
...@@ -139,8 +182,19 @@ test "math.frexp64" {...@@ -139,8 +182,19 @@ test "math.frexp64" {
139 try expect(math.approxEqAbs(f64, r.significand, 0.609558, epsilon) and r.exponent == 7);182 try expect(math.approxEqAbs(f64, r.significand, 0.609558, epsilon) and r.exponent == 7);
140}183}
141184
142test "math.frexp32.special" {185test "128" {
143 var r: frexp32_result = undefined;186 const epsilon = 0.000001;
187 var r: Frexp(f128) = undefined;
188
189 r = frexp128(1.3);
190 try expect(math.approxEqAbs(f128, r.significand, 0.65, epsilon) and r.exponent == 1);
191
192 r = frexp128(78.0234);
193 try expect(math.approxEqAbs(f128, r.significand, 0.609558, epsilon) and r.exponent == 7);
194}
195
196test "32 special" {
197 var r: Frexp(f32) = undefined;
144198
145 r = frexp32(0.0);199 r = frexp32(0.0);
146 try expect(r.significand == 0.0 and r.exponent == 0);200 try expect(r.significand == 0.0 and r.exponent == 0);
...@@ -158,8 +212,8 @@ test "math.frexp32.special" {...@@ -158,8 +212,8 @@ test "math.frexp32.special" {
158 try expect(math.isNan(r.significand));212 try expect(math.isNan(r.significand));
159}213}
160214
161test "math.frexp64.special" {215test "64 special" {
162 var r: frexp64_result = undefined;216 var r: Frexp(f64) = undefined;
163217
164 r = frexp64(0.0);218 r = frexp64(0.0);
165 try expect(r.significand == 0.0 and r.exponent == 0);219 try expect(r.significand == 0.0 and r.exponent == 0);
...@@ -176,3 +230,22 @@ test "math.frexp64.special" {...@@ -176,3 +230,22 @@ test "math.frexp64.special" {
176 r = frexp64(math.nan(f64));230 r = frexp64(math.nan(f64));
177 try expect(math.isNan(r.significand));231 try expect(math.isNan(r.significand));
178}232}
233
234test "128 special" {
235 var r: Frexp(f128) = undefined;
236
237 r = frexp128(0.0);
238 try expect(r.significand == 0.0 and r.exponent == 0);
239
240 r = frexp128(-0.0);
241 try expect(r.significand == -0.0 and r.exponent == 0);
242
243 r = frexp128(math.inf(f128));
244 try expect(math.isPositiveInf(r.significand) and r.exponent == 0);
245
246 r = frexp128(-math.inf(f128));
247 try expect(math.isNegativeInf(r.significand) and r.exponent == 0);
248
249 r = frexp128(math.nan(f128));
250 try expect(math.isNan(r.significand));
251}
lib/std/math/ilogb.zig+57-6
...@@ -1,6 +1,7 @@...@@ -1,6 +1,7 @@
1// Ported from musl, which is licensed under the MIT license:1// Ported from musl, which is MIT licensed.
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogbl.c
4// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogbf.c5// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogbf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogb.c6// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogb.c
67
...@@ -21,10 +22,13 @@ pub fn ilogb(x: anytype) i32 {...@@ -21,10 +22,13 @@ pub fn ilogb(x: anytype) i32 {
21 return switch (T) {22 return switch (T) {
22 f32 => ilogb32(x),23 f32 => ilogb32(x),
23 f64 => ilogb64(x),24 f64 => ilogb64(x),
25 f128 => ilogb128(x),
24 else => @compileError("ilogb not implemented for " ++ @typeName(T)),26 else => @compileError("ilogb not implemented for " ++ @typeName(T)),
25 };27 };
26}28}
2729
30// TODO: unify these implementations with generics
31
28// NOTE: Should these be exposed publicly?32// NOTE: Should these be exposed publicly?
29const fp_ilogbnan = -1 - @as(i32, maxInt(u32) >> 1);33const fp_ilogbnan = -1 - @as(i32, maxInt(u32) >> 1);
30const fp_ilogb0 = fp_ilogbnan;34const fp_ilogb0 = fp_ilogbnan;
...@@ -100,12 +104,43 @@ fn ilogb64(x: f64) i32 {...@@ -100,12 +104,43 @@ fn ilogb64(x: f64) i32 {
100 return e - 0x3FF;104 return e - 0x3FF;
101}105}
102106
103test "math.ilogb" {107fn ilogb128(x: f128) i32 {
108 var u = @bitCast(u128, x);
109 var e = @intCast(i32, (u >> 112) & 0x7FFF);
110
111 if (math.isNan(x)) {
112 return maxInt(i32);
113 }
114
115 if (e == 0) {
116 u <<= 16;
117 if (u == 0) {
118 math.raiseInvalid();
119 return fp_ilogb0;
120 }
121
122 // subnormal x
123 return ilogb128(x * 0x1p120) - 120;
124 }
125
126 if (e == 0x7FFF) {
127 math.raiseInvalid();
128 if (u << 16 != 0) {
129 return fp_ilogbnan;
130 } else {
131 return maxInt(i32);
132 }
133 }
134
135 return e - 0x3FFF;
136}
137
138test "type dispatch" {
104 try expect(ilogb(@as(f32, 0.2)) == ilogb32(0.2));139 try expect(ilogb(@as(f32, 0.2)) == ilogb32(0.2));
105 try expect(ilogb(@as(f64, 0.2)) == ilogb64(0.2));140 try expect(ilogb(@as(f64, 0.2)) == ilogb64(0.2));
106}141}
107142
108test "math.ilogb32" {143test "32" {
109 try expect(ilogb32(0.0) == fp_ilogb0);144 try expect(ilogb32(0.0) == fp_ilogb0);
110 try expect(ilogb32(0.5) == -1);145 try expect(ilogb32(0.5) == -1);
111 try expect(ilogb32(0.8923) == -1);146 try expect(ilogb32(0.8923) == -1);
...@@ -114,7 +149,7 @@ test "math.ilogb32" {...@@ -114,7 +149,7 @@ test "math.ilogb32" {
114 try expect(ilogb32(2398.23) == 11);149 try expect(ilogb32(2398.23) == 11);
115}150}
116151
117test "math.ilogb64" {152test "64" {
118 try expect(ilogb64(0.0) == fp_ilogb0);153 try expect(ilogb64(0.0) == fp_ilogb0);
119 try expect(ilogb64(0.5) == -1);154 try expect(ilogb64(0.5) == -1);
120 try expect(ilogb64(0.8923) == -1);155 try expect(ilogb64(0.8923) == -1);
...@@ -123,16 +158,32 @@ test "math.ilogb64" {...@@ -123,16 +158,32 @@ test "math.ilogb64" {
123 try expect(ilogb64(2398.23) == 11);158 try expect(ilogb64(2398.23) == 11);
124}159}
125160
126test "math.ilogb32.special" {161test "128" {
162 try expect(ilogb128(0.0) == fp_ilogb0);
163 try expect(ilogb128(0.5) == -1);
164 try expect(ilogb128(0.8923) == -1);
165 try expect(ilogb128(10.0) == 3);
166 try expect(ilogb128(-123984) == 16);
167 try expect(ilogb128(2398.23) == 11);
168}
169
170test "32 special" {
127 try expect(ilogb32(math.inf(f32)) == maxInt(i32));171 try expect(ilogb32(math.inf(f32)) == maxInt(i32));
128 try expect(ilogb32(-math.inf(f32)) == maxInt(i32));172 try expect(ilogb32(-math.inf(f32)) == maxInt(i32));
129 try expect(ilogb32(0.0) == minInt(i32));173 try expect(ilogb32(0.0) == minInt(i32));
130 try expect(ilogb32(math.nan(f32)) == maxInt(i32));174 try expect(ilogb32(math.nan(f32)) == maxInt(i32));
131}175}
132176
133test "math.ilogb64.special" {177test "64 special" {
134 try expect(ilogb64(math.inf(f64)) == maxInt(i32));178 try expect(ilogb64(math.inf(f64)) == maxInt(i32));
135 try expect(ilogb64(-math.inf(f64)) == maxInt(i32));179 try expect(ilogb64(-math.inf(f64)) == maxInt(i32));
136 try expect(ilogb64(0.0) == minInt(i32));180 try expect(ilogb64(0.0) == minInt(i32));
137 try expect(ilogb64(math.nan(f64)) == maxInt(i32));181 try expect(ilogb64(math.nan(f64)) == maxInt(i32));
138}182}
183
184test "128 special" {
185 try expect(ilogb128(math.inf(f128)) == maxInt(i32));
186 try expect(ilogb128(-math.inf(f128)) == maxInt(i32));
187 try expect(ilogb128(0.0) == minInt(i32));
188 try expect(ilogb128(math.nan(f128)) == maxInt(i32));
189}
lib/std/special/c_stage1.zig+4
...@@ -656,6 +656,10 @@ export fn ceil(x: f64) f64 {...@@ -656,6 +656,10 @@ export fn ceil(x: f64) f64 {
656 return math.ceil(x);656 return math.ceil(x);
657}657}
658658
659export fn fmal(a: f128, b: f128, c: f128) f128 {
660 return math.fma(f128, a, b, c);
661}
662
659export fn fma(a: f64, b: f64, c: f64) f64 {663export fn fma(a: f64, b: f64, c: f64) f64 {
660 return math.fma(f64, a, b, c);664 return math.fma(f64, a, b, c);
661}665}
test/behavior/muladd.zig+6-7
...@@ -24,11 +24,10 @@ fn testMulAdd() !void {...@@ -24,11 +24,10 @@ fn testMulAdd() !void {
24 var c: f64 = 6.25;24 var c: f64 = 6.25;
25 try expect(@mulAdd(f64, a, b, c) == 20);25 try expect(@mulAdd(f64, a, b, c) == 20);
26 }26 }
27 // Awaits implementation in libm.zig27 {
28 //{28 var a: f16 = 5.5;
29 // var a: f16 = 5.5;29 var b: f128 = 2.5;
30 // var b: f128 = 2.5;30 var c: f128 = 6.25;
31 // var c: f128 = 6.25;31 try expect(@mulAdd(f128, a, b, c) == 20);
32 //try expect(@mulAdd(f128, a, b, c) == 20);32 }
33 //}
34}33}