master
  1//! Ported from musl, which is licensed under the MIT license:
  2//! https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
  3//!
  4//! https://git.musl-libc.org/cgit/musl/tree/src/math/roundf.c
  5//! https://git.musl-libc.org/cgit/musl/tree/src/math/round.c
  6
  7const std = @import("std");
  8const builtin = @import("builtin");
  9const math = std.math;
 10const mem = std.mem;
 11const expect = std.testing.expect;
 12const arch = builtin.cpu.arch;
 13const common = @import("common.zig");
 14
 15pub const panic = common.panic;
 16
 17comptime {
 18    @export(&__roundh, .{ .name = "__roundh", .linkage = common.linkage, .visibility = common.visibility });
 19    @export(&roundf, .{ .name = "roundf", .linkage = common.linkage, .visibility = common.visibility });
 20    @export(&round, .{ .name = "round", .linkage = common.linkage, .visibility = common.visibility });
 21    @export(&__roundx, .{ .name = "__roundx", .linkage = common.linkage, .visibility = common.visibility });
 22    if (common.want_ppc_abi) {
 23        @export(&roundq, .{ .name = "roundf128", .linkage = common.linkage, .visibility = common.visibility });
 24    }
 25    @export(&roundq, .{ .name = "roundq", .linkage = common.linkage, .visibility = common.visibility });
 26    @export(&roundl, .{ .name = "roundl", .linkage = common.linkage, .visibility = common.visibility });
 27}
 28
 29pub fn __roundh(x: f16) callconv(.c) f16 {
 30    // TODO: more efficient implementation
 31    return @floatCast(roundf(x));
 32}
 33
 34pub fn roundf(x_: f32) callconv(.c) f32 {
 35    const f32_toint = 1.0 / math.floatEps(f32);
 36
 37    var x = x_;
 38    const u: u32 = @bitCast(x);
 39    const e = (u >> 23) & 0xFF;
 40    var y: f32 = undefined;
 41
 42    if (e >= 0x7F + 23) {
 43        return x;
 44    }
 45    if (u >> 31 != 0) {
 46        x = -x;
 47    }
 48    if (e < 0x7F - 1) {
 49        if (common.want_float_exceptions) mem.doNotOptimizeAway(x + f32_toint);
 50        return 0 * @as(f32, @bitCast(u));
 51    }
 52
 53    y = x + f32_toint - f32_toint - x;
 54    if (y > 0.5) {
 55        y = y + x - 1;
 56    } else if (y <= -0.5) {
 57        y = y + x + 1;
 58    } else {
 59        y = y + x;
 60    }
 61
 62    if (u >> 31 != 0) {
 63        return -y;
 64    } else {
 65        return y;
 66    }
 67}
 68
 69pub fn round(x_: f64) callconv(.c) f64 {
 70    const f64_toint = 1.0 / math.floatEps(f64);
 71
 72    var x = x_;
 73    const u: u64 = @bitCast(x);
 74    const e = (u >> 52) & 0x7FF;
 75    var y: f64 = undefined;
 76
 77    if (e >= 0x3FF + 52) {
 78        return x;
 79    }
 80    if (u >> 63 != 0) {
 81        x = -x;
 82    }
 83    if (e < 0x3ff - 1) {
 84        if (common.want_float_exceptions) mem.doNotOptimizeAway(x + f64_toint);
 85        return 0 * @as(f64, @bitCast(u));
 86    }
 87
 88    y = x + f64_toint - f64_toint - x;
 89    if (y > 0.5) {
 90        y = y + x - 1;
 91    } else if (y <= -0.5) {
 92        y = y + x + 1;
 93    } else {
 94        y = y + x;
 95    }
 96
 97    if (u >> 63 != 0) {
 98        return -y;
 99    } else {
100        return y;
101    }
102}
103
104pub fn __roundx(x: f80) callconv(.c) f80 {
105    // TODO: more efficient implementation
106    return @floatCast(roundq(x));
107}
108
109pub fn roundq(x_: f128) callconv(.c) f128 {
110    const f128_toint = 1.0 / math.floatEps(f128);
111
112    var x = x_;
113    const u: u128 = @bitCast(x);
114    const e = (u >> 112) & 0x7FFF;
115    var y: f128 = undefined;
116
117    if (e >= 0x3FFF + 112) {
118        return x;
119    }
120    if (u >> 127 != 0) {
121        x = -x;
122    }
123    if (e < 0x3FFF - 1) {
124        if (common.want_float_exceptions) mem.doNotOptimizeAway(x + f128_toint);
125        return 0 * @as(f128, @bitCast(u));
126    }
127
128    y = x + f128_toint - f128_toint - x;
129    if (y > 0.5) {
130        y = y + x - 1;
131    } else if (y <= -0.5) {
132        y = y + x + 1;
133    } else {
134        y = y + x;
135    }
136
137    if (u >> 127 != 0) {
138        return -y;
139    } else {
140        return y;
141    }
142}
143
144pub fn roundl(x: c_longdouble) callconv(.c) c_longdouble {
145    switch (@typeInfo(c_longdouble).float.bits) {
146        16 => return __roundh(x),
147        32 => return roundf(x),
148        64 => return round(x),
149        80 => return __roundx(x),
150        128 => return roundq(x),
151        else => @compileError("unreachable"),
152    }
153}
154
155test "round32" {
156    try expect(roundf(1.3) == 1.0);
157    try expect(roundf(-1.3) == -1.0);
158    try expect(roundf(0.2) == 0.0);
159    try expect(roundf(1.8) == 2.0);
160}
161
162test "round64" {
163    try expect(round(1.3) == 1.0);
164    try expect(round(-1.3) == -1.0);
165    try expect(round(0.2) == 0.0);
166    try expect(round(1.8) == 2.0);
167}
168
169test "round128" {
170    try expect(roundq(1.3) == 1.0);
171    try expect(roundq(-1.3) == -1.0);
172    try expect(roundq(0.2) == 0.0);
173    try expect(roundq(1.8) == 2.0);
174}
175
176test "round32.special" {
177    try expect(roundf(0.0) == 0.0);
178    try expect(roundf(-0.0) == -0.0);
179    try expect(math.isPositiveInf(roundf(math.inf(f32))));
180    try expect(math.isNegativeInf(roundf(-math.inf(f32))));
181    try expect(math.isNan(roundf(math.nan(f32))));
182}
183
184test "round64.special" {
185    try expect(round(0.0) == 0.0);
186    try expect(round(-0.0) == -0.0);
187    try expect(math.isPositiveInf(round(math.inf(f64))));
188    try expect(math.isNegativeInf(round(-math.inf(f64))));
189    try expect(math.isNan(round(math.nan(f64))));
190}
191
192test "round128.special" {
193    try expect(roundq(0.0) == 0.0);
194    try expect(roundq(-0.0) == -0.0);
195    try expect(math.isPositiveInf(roundq(math.inf(f128))));
196    try expect(math.isNegativeInf(roundq(-math.inf(f128))));
197    try expect(math.isNan(roundq(math.nan(f128))));
198}