mirror of
https://codeberg.org/ziglang/zig.git
synced 2026-04-27 19:09:47 +03:00
Merge pull request 'libzigc: Implement tanhf, tanh, modfl, modff, modf, acoshf, asin' (#31536) from mihael/zig:libzigc/more-math-functions into master
Reviewed-on: https://codeberg.org/ziglang/zig/pulls/31536 Reviewed-by: Andrew Kelley <andrew@ziglang.org>
This commit is contained in:
Andrew Kelley
+133
-17
@@ -2,6 +2,10 @@ const builtin = @import("builtin");
|
||||
|
||||
const std = @import("std");
|
||||
const math = std.math;
|
||||
const expect = std.testing.expect;
|
||||
const expectEqual = std.testing.expectEqual;
|
||||
const expectApproxEqAbs = std.testing.expectApproxEqAbs;
|
||||
const expectApproxEqRel = std.testing.expectApproxEqRel;
|
||||
|
||||
const symbol = @import("../c.zig").symbol;
|
||||
|
||||
@@ -34,14 +38,19 @@ comptime {
|
||||
symbol(&coshf, "coshf");
|
||||
symbol(&hypotf, "hypotf");
|
||||
symbol(&hypotl, "hypotl");
|
||||
symbol(&modff, "modff");
|
||||
symbol(&modfl, "modfl");
|
||||
symbol(&nan, "nan");
|
||||
symbol(&nanf, "nanf");
|
||||
symbol(&nanl, "nanl");
|
||||
symbol(&tanhf, "tanhf");
|
||||
}
|
||||
|
||||
if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
|
||||
symbol(&acos, "acos");
|
||||
symbol(&acosf, "acosf");
|
||||
symbol(&acoshf, "acoshf");
|
||||
symbol(&asin, "asin");
|
||||
symbol(&atan, "atan");
|
||||
symbol(&atanf, "atanf");
|
||||
symbol(&atanl, "atanl");
|
||||
@@ -51,9 +60,11 @@ comptime {
|
||||
symbol(&exp10, "exp10");
|
||||
symbol(&exp10f, "exp10f");
|
||||
symbol(&hypot, "hypot");
|
||||
symbol(&modf, "modf");
|
||||
symbol(&pow, "pow");
|
||||
symbol(&pow10, "pow10");
|
||||
symbol(&pow10f, "pow10f");
|
||||
symbol(&tanh, "tanh");
|
||||
}
|
||||
|
||||
if (builtin.target.isMuslLibC()) {
|
||||
@@ -70,7 +81,15 @@ fn acos(x: f64) callconv(.c) f64 {
|
||||
}
|
||||
|
||||
fn acosf(x: f32) callconv(.c) f32 {
|
||||
return std.math.acos(x);
|
||||
return math.acos(x);
|
||||
}
|
||||
|
||||
fn acoshf(x: f32) callconv(.c) f32 {
|
||||
return math.acosh(x);
|
||||
}
|
||||
|
||||
fn asin(x: f64) callconv(.c) f64 {
|
||||
return math.asin(x);
|
||||
}
|
||||
|
||||
fn atan(x: f64) callconv(.c) f64 {
|
||||
@@ -152,6 +171,95 @@ fn isnanl(x: c_longdouble) callconv(.c) c_int {
|
||||
return if (math.isNan(x)) 1 else 0;
|
||||
}
|
||||
|
||||
fn modfGeneric(comptime T: type, x: T, iptr: *T) T {
|
||||
if (math.isNegativeInf(x)) {
|
||||
iptr.* = -math.inf(T);
|
||||
return -0.0;
|
||||
}
|
||||
|
||||
if (math.isPositiveInf(x)) {
|
||||
iptr.* = math.inf(T);
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
if (math.isNan(x)) {
|
||||
iptr.* = math.nan(T);
|
||||
return math.nan(T);
|
||||
}
|
||||
|
||||
const r = math.modf(x);
|
||||
iptr.* = r.ipart;
|
||||
|
||||
// If the result is a negative zero, we must be explicit about
|
||||
// returning a negative zero.
|
||||
return if (math.isNegativeZero(x) or (x < 0.0 and x == r.ipart)) -0.0 else r.fpart;
|
||||
}
|
||||
|
||||
fn modf(x: f64, iptr: *f64) callconv(.c) f64 {
|
||||
return modfGeneric(f64, x, iptr);
|
||||
}
|
||||
|
||||
fn modff(x: f32, iptr: *f32) callconv(.c) f32 {
|
||||
return modfGeneric(f32, x, iptr);
|
||||
}
|
||||
|
||||
fn modfl(x: c_longdouble, iptr: *c_longdouble) callconv(.c) c_longdouble {
|
||||
return modfGeneric(c_longdouble, x, iptr);
|
||||
}
|
||||
|
||||
fn testModf(comptime T: type) !void {
|
||||
// Choose the appropriate `modf` impl to test based on type
|
||||
const f = switch (T) {
|
||||
f32 => modff,
|
||||
f64 => modf,
|
||||
c_longdouble => modfl,
|
||||
else => @compileError("modf not implemented for " ++ @typeName(T)),
|
||||
};
|
||||
|
||||
var int: T = undefined;
|
||||
const iptr = ∫
|
||||
const eps_val: comptime_float = @max(1e-6, math.floatEps(T));
|
||||
|
||||
const normal_frac = f(@as(T, 1234.567), iptr);
|
||||
// Account for precision error
|
||||
const expected = 1234.567 - @as(T, 1234);
|
||||
try expectApproxEqAbs(expected, normal_frac, eps_val);
|
||||
try expectApproxEqRel(@as(T, 1234.0), iptr.*, eps_val);
|
||||
|
||||
// When `x` is a NaN, NaN is returned and `*iptr` is set to NaN
|
||||
const nan_frac = f(math.nan(T), iptr);
|
||||
try expect(math.isNan(nan_frac));
|
||||
try expect(math.isNan(iptr.*));
|
||||
|
||||
// When `x` is positive infinity, +0 is returned and `*iptr` is set to
|
||||
// positive infinity
|
||||
const pos_zero_frac = f(math.inf(T), iptr);
|
||||
try expect(math.isPositiveZero(pos_zero_frac));
|
||||
try expect(math.isPositiveInf(iptr.*));
|
||||
|
||||
// When `x` is negative infinity, -0 is returned and `*iptr` is set to
|
||||
// negative infinity
|
||||
const neg_zero_frac = f(-math.inf(T), iptr);
|
||||
try expect(math.isNegativeZero(neg_zero_frac));
|
||||
try expect(math.isNegativeInf(iptr.*));
|
||||
|
||||
// Return -0 when `x` is a negative integer
|
||||
const nz_frac = f(@as(T, -1000.0), iptr);
|
||||
try expect(math.isNegativeZero(nz_frac));
|
||||
try expectEqual(@as(T, -1000.0), iptr.*);
|
||||
|
||||
// Return +0 when `x` is a positive integer
|
||||
const pz_frac = f(@as(T, 1000.0), iptr);
|
||||
try expect(math.isPositiveZero(pz_frac));
|
||||
try expectEqual(@as(T, 1000.0), iptr.*);
|
||||
}
|
||||
|
||||
test "modf" {
|
||||
try testModf(f32);
|
||||
try testModf(f64);
|
||||
try testModf(c_longdouble);
|
||||
}
|
||||
|
||||
fn nan(_: [*:0]const c_char) callconv(.c) f64 {
|
||||
return math.nan(f64);
|
||||
}
|
||||
@@ -177,7 +285,7 @@ fn pow10f(x: f32) callconv(.c) f32 {
|
||||
}
|
||||
|
||||
fn rint(x: f64) callconv(.c) f64 {
|
||||
const toint: f64 = 1.0 / @as(f64, std.math.floatEps(f64));
|
||||
const toint: f64 = 1.0 / @as(f64, math.floatEps(f64));
|
||||
const a: u64 = @bitCast(x);
|
||||
const e = a >> 52 & 0x7ff;
|
||||
const s = a >> 63;
|
||||
@@ -199,35 +307,43 @@ fn rint(x: f64) callconv(.c) f64 {
|
||||
|
||||
test "rint" {
|
||||
// Positive numbers round correctly
|
||||
try std.testing.expectEqual(@as(f64, 42.0), rint(42.2));
|
||||
try std.testing.expectEqual(@as(f64, 42.0), rint(41.8));
|
||||
try expectEqual(@as(f64, 42.0), rint(42.2));
|
||||
try expectEqual(@as(f64, 42.0), rint(41.8));
|
||||
|
||||
// Negative numbers round correctly
|
||||
try std.testing.expectEqual(@as(f64, -6.0), rint(-5.9));
|
||||
try std.testing.expectEqual(@as(f64, -6.0), rint(-6.1));
|
||||
try expectEqual(@as(f64, -6.0), rint(-5.9));
|
||||
try expectEqual(@as(f64, -6.0), rint(-6.1));
|
||||
|
||||
// No rounding needed test
|
||||
try std.testing.expectEqual(@as(f64, 5.0), rint(5.0));
|
||||
try std.testing.expectEqual(@as(f64, -10.0), rint(-10.0));
|
||||
try std.testing.expectEqual(@as(f64, 0.0), rint(0.0));
|
||||
try expectEqual(@as(f64, 5.0), rint(5.0));
|
||||
try expectEqual(@as(f64, -10.0), rint(-10.0));
|
||||
try expectEqual(@as(f64, 0.0), rint(0.0));
|
||||
|
||||
// Very large numbers return unchanged
|
||||
const large: f64 = 9007199254740992.0; // 2^53
|
||||
try std.testing.expectEqual(large, rint(large));
|
||||
try std.testing.expectEqual(-large, rint(-large));
|
||||
try expectEqual(large, rint(large));
|
||||
try expectEqual(-large, rint(-large));
|
||||
|
||||
// Small positive numbers round to zero
|
||||
const pos_result = rint(0.3);
|
||||
try std.testing.expectEqual(@as(f64, 0.0), pos_result);
|
||||
try std.testing.expect(@as(u64, @bitCast(pos_result)) == 0);
|
||||
try expectEqual(@as(f64, 0.0), pos_result);
|
||||
try expect(@as(u64, @bitCast(pos_result)) == 0);
|
||||
|
||||
// Small negative numbers round to negative zero
|
||||
const neg_result = rint(-0.3);
|
||||
try std.testing.expectEqual(@as(f64, 0.0), neg_result);
|
||||
try expectEqual(@as(f64, 0.0), neg_result);
|
||||
const bits: u64 = @bitCast(neg_result);
|
||||
try std.testing.expect((bits >> 63) == 1);
|
||||
try expect((bits >> 63) == 1);
|
||||
|
||||
// Exact half rounds to nearest even (banker's rounding)
|
||||
try std.testing.expectEqual(@as(f64, 2.0), rint(2.5));
|
||||
try std.testing.expectEqual(@as(f64, 4.0), rint(3.5));
|
||||
try expectEqual(@as(f64, 2.0), rint(2.5));
|
||||
try expectEqual(@as(f64, 4.0), rint(3.5));
|
||||
}
|
||||
|
||||
fn tanh(x: f64) callconv(.c) f64 {
|
||||
return math.tanh(x);
|
||||
}
|
||||
|
||||
fn tanhf(x: f32) callconv(.c) f32 {
|
||||
return math.tanh(x);
|
||||
}
|
||||
|
||||
Vendored
-42
@@ -1,42 +0,0 @@
|
||||
/**
|
||||
* This file has no copyright assigned and is placed in the Public Domain.
|
||||
* This file is part of the mingw-w64 runtime package.
|
||||
* No warranty is given; refer to the file DISCLAIMER.PD within this package.
|
||||
*/
|
||||
#include <fenv.h>
|
||||
#include <math.h>
|
||||
#include <errno.h>
|
||||
|
||||
float
|
||||
modff (float value, float* iptr)
|
||||
{
|
||||
float int_part = 0.0F;
|
||||
/* truncate */
|
||||
/* truncate */
|
||||
#if (defined(_AMD64_) && !defined(_ARM64EC_)) || (defined(__x86_64__) && !defined(__arm64ec__))
|
||||
asm volatile ("subq $8, %%rsp\n"
|
||||
"fnstcw 4(%%rsp)\n"
|
||||
"movzwl 4(%%rsp), %%eax\n"
|
||||
"orb $12, %%ah\n"
|
||||
"movw %%ax, (%%rsp)\n"
|
||||
"fldcw (%%rsp)\n"
|
||||
"frndint\n"
|
||||
"fldcw 4(%%rsp)\n"
|
||||
"addq $8, %%rsp\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
|
||||
#elif defined(_X86_) || defined(__i386__)
|
||||
asm volatile ("push %%eax\n\tsubl $8, %%esp\n"
|
||||
"fnstcw 4(%%esp)\n"
|
||||
"movzwl 4(%%esp), %%eax\n"
|
||||
"orb $12, %%ah\n"
|
||||
"movw %%ax, (%%esp)\n"
|
||||
"fldcw (%%esp)\n"
|
||||
"frndint\n"
|
||||
"fldcw 4(%%esp)\n"
|
||||
"addl $8, %%esp\n\tpop %%eax\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
|
||||
#else
|
||||
int_part = truncf(value);
|
||||
#endif
|
||||
if (iptr)
|
||||
*iptr = int_part;
|
||||
return (isinf (value) ? 0.0F : value - int_part);
|
||||
}
|
||||
Vendored
-41
@@ -1,41 +0,0 @@
|
||||
/**
|
||||
* This file has no copyright assigned and is placed in the Public Domain.
|
||||
* This file is part of the mingw-w64 runtime package.
|
||||
* No warranty is given; refer to the file DISCLAIMER.PD within this package.
|
||||
*/
|
||||
#include <fenv.h>
|
||||
#include <math.h>
|
||||
#include <errno.h>
|
||||
|
||||
long double
|
||||
modfl (long double value, long double* iptr)
|
||||
{
|
||||
long double int_part = 0.0L;
|
||||
/* truncate */
|
||||
#if (defined(_AMD64_) && !defined(_ARM64EC_)) || (defined(__x86_64__) && !defined(__arm64ec__))
|
||||
asm volatile ("subq $8, %%rsp\n"
|
||||
"fnstcw 4(%%rsp)\n"
|
||||
"movzwl 4(%%rsp), %%eax\n"
|
||||
"orb $12, %%ah\n"
|
||||
"movw %%ax, (%%rsp)\n"
|
||||
"fldcw (%%rsp)\n"
|
||||
"frndint\n"
|
||||
"fldcw 4(%%rsp)\n"
|
||||
"addq $8, %%rsp\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
|
||||
#elif defined(_X86_) || defined(__i386__)
|
||||
asm volatile ("push %%eax\n\tsubl $8, %%esp\n"
|
||||
"fnstcw 4(%%esp)\n"
|
||||
"movzwl 4(%%esp), %%eax\n"
|
||||
"orb $12, %%ah\n"
|
||||
"movw %%ax, (%%esp)\n"
|
||||
"fldcw (%%esp)\n"
|
||||
"frndint\n"
|
||||
"fldcw 4(%%esp)\n"
|
||||
"addl $8, %%esp\n\tpop %%eax\n" : "=t" (int_part) : "0" (value) : "eax"); /* round */
|
||||
#else
|
||||
int_part = truncl(value);
|
||||
#endif
|
||||
if (iptr)
|
||||
*iptr = int_part;
|
||||
return (isinf (value) ? 0.0L : value - int_part);
|
||||
}
|
||||
Vendored
-10
@@ -1,10 +0,0 @@
|
||||
/**
|
||||
* This file has no copyright assigned and is placed in the Public Domain.
|
||||
* This file is part of the mingw-w64 runtime package.
|
||||
* No warranty is given; refer to the file DISCLAIMER.PD within this package.
|
||||
*/
|
||||
#include <math.h>
|
||||
float tanhf (float x)
|
||||
{
|
||||
return (float) tanh (x);
|
||||
}
|
||||
Vendored
-26
@@ -1,26 +0,0 @@
|
||||
#include "libm.h"
|
||||
|
||||
#if FLT_EVAL_METHOD==2
|
||||
#undef sqrtf
|
||||
#define sqrtf sqrtl
|
||||
#elif FLT_EVAL_METHOD==1
|
||||
#undef sqrtf
|
||||
#define sqrtf sqrt
|
||||
#endif
|
||||
|
||||
/* acosh(x) = log(x + sqrt(x*x-1)) */
|
||||
float acoshf(float x)
|
||||
{
|
||||
union {float f; uint32_t i;} u = {x};
|
||||
uint32_t a = u.i & 0x7fffffff;
|
||||
|
||||
if (a < 0x3f800000+(1<<23))
|
||||
/* |x| < 2, invalid if x < 1 */
|
||||
/* up to 2ulp error in [1,1.125] */
|
||||
return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1)));
|
||||
if (u.i < 0x3f800000+(12<<23))
|
||||
/* 2 <= x < 0x1p12 */
|
||||
return logf(2*x - 1/(x+sqrtf(x*x-1)));
|
||||
/* x >= 0x1p12 or x <= -2 or nan */
|
||||
return logf(x) + 0.693147180559945309417232121458176568f;
|
||||
}
|
||||
Vendored
-107
@@ -1,107 +0,0 @@
|
||||
/* origin: FreeBSD /usr/src/lib/msun/src/e_asin.c */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunSoft, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
/* asin(x)
|
||||
* Method :
|
||||
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
|
||||
* we approximate asin(x) on [0,0.5] by
|
||||
* asin(x) = x + x*x^2*R(x^2)
|
||||
* where
|
||||
* R(x^2) is a rational approximation of (asin(x)-x)/x^3
|
||||
* and its remez error is bounded by
|
||||
* |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
|
||||
*
|
||||
* For x in [0.5,1]
|
||||
* asin(x) = pi/2-2*asin(sqrt((1-x)/2))
|
||||
* Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
|
||||
* then for x>0.98
|
||||
* asin(x) = pi/2 - 2*(s+s*z*R(z))
|
||||
* = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
|
||||
* For x<=0.98, let pio4_hi = pio2_hi/2, then
|
||||
* f = hi part of s;
|
||||
* c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
|
||||
* and
|
||||
* asin(x) = pi/2 - 2*(s+s*z*R(z))
|
||||
* = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
|
||||
* = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
|
||||
*
|
||||
* Special cases:
|
||||
* if x is NaN, return x itself;
|
||||
* if |x|>1, return NaN with invalid signal.
|
||||
*
|
||||
*/
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
static const double
|
||||
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
|
||||
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
|
||||
/* coefficients for R(x^2) */
|
||||
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
|
||||
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
|
||||
pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
|
||||
pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
|
||||
pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
|
||||
pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
|
||||
qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
|
||||
qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
|
||||
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
|
||||
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
|
||||
|
||||
static double R(double z)
|
||||
{
|
||||
double_t p, q;
|
||||
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
|
||||
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
|
||||
return p/q;
|
||||
}
|
||||
|
||||
double asin(double x)
|
||||
{
|
||||
double z,r,s;
|
||||
uint32_t hx,ix;
|
||||
|
||||
GET_HIGH_WORD(hx, x);
|
||||
ix = hx & 0x7fffffff;
|
||||
/* |x| >= 1 or nan */
|
||||
if (ix >= 0x3ff00000) {
|
||||
uint32_t lx;
|
||||
GET_LOW_WORD(lx, x);
|
||||
if ((ix-0x3ff00000 | lx) == 0)
|
||||
/* asin(1) = +-pi/2 with inexact */
|
||||
return x*pio2_hi + 0x1p-120f;
|
||||
return 0/(x-x);
|
||||
}
|
||||
/* |x| < 0.5 */
|
||||
if (ix < 0x3fe00000) {
|
||||
/* if 0x1p-1022 <= |x| < 0x1p-26, avoid raising underflow */
|
||||
if (ix < 0x3e500000 && ix >= 0x00100000)
|
||||
return x;
|
||||
return x + x*R(x*x);
|
||||
}
|
||||
/* 1 > |x| >= 0.5 */
|
||||
z = (1 - fabs(x))*0.5;
|
||||
s = sqrt(z);
|
||||
r = R(z);
|
||||
if (ix >= 0x3fef3333) { /* if |x| > 0.975 */
|
||||
x = pio2_hi-(2*(s+s*r)-pio2_lo);
|
||||
} else {
|
||||
double f,c;
|
||||
/* f+c = sqrt(z) */
|
||||
f = s;
|
||||
SET_LOW_WORD(f,0);
|
||||
c = (z-f*f)/(s+f);
|
||||
x = 0.5*pio2_hi - (2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
|
||||
}
|
||||
if (hx >> 31)
|
||||
return -x;
|
||||
return x;
|
||||
}
|
||||
Vendored
-21
@@ -1,21 +0,0 @@
|
||||
.global asin
|
||||
.type asin,@function
|
||||
asin:
|
||||
fldl 4(%esp)
|
||||
mov 8(%esp),%eax
|
||||
add %eax,%eax
|
||||
cmp $0x00200000,%eax
|
||||
jb 1f
|
||||
fld %st(0)
|
||||
fld1
|
||||
fsub %st(0),%st(1)
|
||||
fadd %st(2)
|
||||
fmulp
|
||||
fsqrt
|
||||
fpatan
|
||||
fstpl 4(%esp)
|
||||
fldl 4(%esp)
|
||||
ret
|
||||
# subnormal x, return x with underflow
|
||||
1: fsts 4(%esp)
|
||||
ret
|
||||
Vendored
-34
@@ -1,34 +0,0 @@
|
||||
#include "libm.h"
|
||||
|
||||
double modf(double x, double *iptr)
|
||||
{
|
||||
union {double f; uint64_t i;} u = {x};
|
||||
uint64_t mask;
|
||||
int e = (int)(u.i>>52 & 0x7ff) - 0x3ff;
|
||||
|
||||
/* no fractional part */
|
||||
if (e >= 52) {
|
||||
*iptr = x;
|
||||
if (e == 0x400 && u.i<<12 != 0) /* nan */
|
||||
return x;
|
||||
u.i &= 1ULL<<63;
|
||||
return u.f;
|
||||
}
|
||||
|
||||
/* no integral part*/
|
||||
if (e < 0) {
|
||||
u.i &= 1ULL<<63;
|
||||
*iptr = u.f;
|
||||
return x;
|
||||
}
|
||||
|
||||
mask = -1ULL>>12>>e;
|
||||
if ((u.i & mask) == 0) {
|
||||
*iptr = x;
|
||||
u.i &= 1ULL<<63;
|
||||
return u.f;
|
||||
}
|
||||
u.i &= ~mask;
|
||||
*iptr = u.f;
|
||||
return x - u.f;
|
||||
}
|
||||
Vendored
-34
@@ -1,34 +0,0 @@
|
||||
#include "libm.h"
|
||||
|
||||
float modff(float x, float *iptr)
|
||||
{
|
||||
union {float f; uint32_t i;} u = {x};
|
||||
uint32_t mask;
|
||||
int e = (int)(u.i>>23 & 0xff) - 0x7f;
|
||||
|
||||
/* no fractional part */
|
||||
if (e >= 23) {
|
||||
*iptr = x;
|
||||
if (e == 0x80 && u.i<<9 != 0) { /* nan */
|
||||
return x;
|
||||
}
|
||||
u.i &= 0x80000000;
|
||||
return u.f;
|
||||
}
|
||||
/* no integral part */
|
||||
if (e < 0) {
|
||||
u.i &= 0x80000000;
|
||||
*iptr = u.f;
|
||||
return x;
|
||||
}
|
||||
|
||||
mask = 0x007fffff>>e;
|
||||
if ((u.i & mask) == 0) {
|
||||
*iptr = x;
|
||||
u.i &= 0x80000000;
|
||||
return u.f;
|
||||
}
|
||||
u.i &= ~mask;
|
||||
*iptr = u.f;
|
||||
return x - u.f;
|
||||
}
|
||||
Vendored
-53
@@ -1,53 +0,0 @@
|
||||
#include "libm.h"
|
||||
|
||||
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
||||
long double modfl(long double x, long double *iptr)
|
||||
{
|
||||
double d;
|
||||
long double r;
|
||||
|
||||
r = modf(x, &d);
|
||||
*iptr = d;
|
||||
return r;
|
||||
}
|
||||
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
|
||||
|
||||
static const long double toint = 1/LDBL_EPSILON;
|
||||
|
||||
long double modfl(long double x, long double *iptr)
|
||||
{
|
||||
union ldshape u = {x};
|
||||
int e = (u.i.se & 0x7fff) - 0x3fff;
|
||||
int s = u.i.se >> 15;
|
||||
long double absx;
|
||||
long double y;
|
||||
|
||||
/* no fractional part */
|
||||
if (e >= LDBL_MANT_DIG-1) {
|
||||
*iptr = x;
|
||||
if (isnan(x))
|
||||
return x;
|
||||
return s ? -0.0 : 0.0;
|
||||
}
|
||||
|
||||
/* no integral part*/
|
||||
if (e < 0) {
|
||||
*iptr = s ? -0.0 : 0.0;
|
||||
return x;
|
||||
}
|
||||
|
||||
/* raises spurious inexact */
|
||||
absx = s ? -x : x;
|
||||
y = absx + toint - toint - absx;
|
||||
if (y == 0) {
|
||||
*iptr = x;
|
||||
return s ? -0.0 : 0.0;
|
||||
}
|
||||
if (y > 0)
|
||||
y -= 1;
|
||||
if (s)
|
||||
y = -y;
|
||||
*iptr = x + y;
|
||||
return -y;
|
||||
}
|
||||
#endif
|
||||
Vendored
-45
@@ -1,45 +0,0 @@
|
||||
#include "libm.h"
|
||||
|
||||
/* tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x))
|
||||
* = (exp(2*x) - 1)/(exp(2*x) - 1 + 2)
|
||||
* = (1 - exp(-2*x))/(exp(-2*x) - 1 + 2)
|
||||
*/
|
||||
double tanh(double x)
|
||||
{
|
||||
union {double f; uint64_t i;} u = {.f = x};
|
||||
uint32_t w;
|
||||
int sign;
|
||||
double_t t;
|
||||
|
||||
/* x = |x| */
|
||||
sign = u.i >> 63;
|
||||
u.i &= (uint64_t)-1/2;
|
||||
x = u.f;
|
||||
w = u.i >> 32;
|
||||
|
||||
if (w > 0x3fe193ea) {
|
||||
/* |x| > log(3)/2 ~= 0.5493 or nan */
|
||||
if (w > 0x40340000) {
|
||||
/* |x| > 20 or nan */
|
||||
/* note: this branch avoids raising overflow */
|
||||
t = 1 - 0/x;
|
||||
} else {
|
||||
t = expm1(2*x);
|
||||
t = 1 - 2/(t+2);
|
||||
}
|
||||
} else if (w > 0x3fd058ae) {
|
||||
/* |x| > log(5/3)/2 ~= 0.2554 */
|
||||
t = expm1(2*x);
|
||||
t = t/(t+2);
|
||||
} else if (w >= 0x00100000) {
|
||||
/* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */
|
||||
t = expm1(-2*x);
|
||||
t = -t/(t+2);
|
||||
} else {
|
||||
/* |x| is subnormal */
|
||||
/* note: the branch above would not raise underflow in [0x1p-1023,0x1p-1022) */
|
||||
FORCE_EVAL((float)x);
|
||||
t = x;
|
||||
}
|
||||
return sign ? -t : t;
|
||||
}
|
||||
Vendored
-39
@@ -1,39 +0,0 @@
|
||||
#include "libm.h"
|
||||
|
||||
float tanhf(float x)
|
||||
{
|
||||
union {float f; uint32_t i;} u = {.f = x};
|
||||
uint32_t w;
|
||||
int sign;
|
||||
float t;
|
||||
|
||||
/* x = |x| */
|
||||
sign = u.i >> 31;
|
||||
u.i &= 0x7fffffff;
|
||||
x = u.f;
|
||||
w = u.i;
|
||||
|
||||
if (w > 0x3f0c9f54) {
|
||||
/* |x| > log(3)/2 ~= 0.5493 or nan */
|
||||
if (w > 0x41200000) {
|
||||
/* |x| > 10 */
|
||||
t = 1 + 0/x;
|
||||
} else {
|
||||
t = expm1f(2*x);
|
||||
t = 1 - 2/(t+2);
|
||||
}
|
||||
} else if (w > 0x3e82c578) {
|
||||
/* |x| > log(5/3)/2 ~= 0.2554 */
|
||||
t = expm1f(2*x);
|
||||
t = t/(t+2);
|
||||
} else if (w >= 0x00800000) {
|
||||
/* |x| >= 0x1p-126 */
|
||||
t = expm1f(-2*x);
|
||||
t = -t/(t+2);
|
||||
} else {
|
||||
/* |x| is subnormal */
|
||||
FORCE_EVAL(x*x);
|
||||
t = x;
|
||||
}
|
||||
return sign ? -t : t;
|
||||
}
|
||||
@@ -619,7 +619,6 @@ const mingw32_generic_src = [_][]const u8{
|
||||
"math" ++ path.sep_str ++ "lgamma.c",
|
||||
"math" ++ path.sep_str ++ "lgammaf.c",
|
||||
"math" ++ path.sep_str ++ "lgammal.c",
|
||||
"math" ++ path.sep_str ++ "modfl.c",
|
||||
"math" ++ path.sep_str ++ "powi.c",
|
||||
"math" ++ path.sep_str ++ "powif.c",
|
||||
"math" ++ path.sep_str ++ "powil.c",
|
||||
@@ -977,10 +976,8 @@ const mingw32_x86_src = [_][]const u8{
|
||||
|
||||
const mingw32_x86_32_src = [_][]const u8{
|
||||
// ucrtbase
|
||||
"math" ++ path.sep_str ++ "modff.c",
|
||||
"math" ++ path.sep_str ++ "powf.c",
|
||||
"math" ++ path.sep_str ++ "sinhf.c",
|
||||
"math" ++ path.sep_str ++ "tanhf.c",
|
||||
"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "acosf.c",
|
||||
"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "asinf.c",
|
||||
"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atan2f.c",
|
||||
|
||||
@@ -796,12 +796,10 @@ const src_files = [_][]const u8{
|
||||
"musl/src/math/aarch64/nearbyintf.c",
|
||||
"musl/src/math/aarch64/rintf.c",
|
||||
"musl/src/math/acosh.c",
|
||||
"musl/src/math/acoshf.c",
|
||||
"musl/src/math/acoshl.c",
|
||||
"musl/src/math/acosl.c",
|
||||
"musl/src/math/arm/fma.c",
|
||||
"musl/src/math/arm/fmaf.c",
|
||||
"musl/src/math/asin.c",
|
||||
"musl/src/math/asinf.c",
|
||||
"musl/src/math/asinh.c",
|
||||
"musl/src/math/asinhf.c",
|
||||
@@ -849,7 +847,6 @@ const src_files = [_][]const u8{
|
||||
"musl/src/math/i386/acosl.s",
|
||||
"musl/src/math/i386/asinf.s",
|
||||
"musl/src/math/i386/asinl.s",
|
||||
"musl/src/math/i386/asin.s",
|
||||
"musl/src/math/i386/atan2f.s",
|
||||
"musl/src/math/i386/atan2l.s",
|
||||
"musl/src/math/i386/atan2.s",
|
||||
@@ -937,9 +934,6 @@ const src_files = [_][]const u8{
|
||||
"musl/src/math/__math_uflowf.c",
|
||||
"musl/src/math/__math_xflow.c",
|
||||
"musl/src/math/__math_xflowf.c",
|
||||
"musl/src/math/modf.c",
|
||||
"musl/src/math/modff.c",
|
||||
"musl/src/math/modfl.c",
|
||||
"musl/src/math/nearbyint.c",
|
||||
"musl/src/math/nearbyintf.c",
|
||||
"musl/src/math/nearbyintl.c",
|
||||
@@ -1009,8 +1003,6 @@ const src_files = [_][]const u8{
|
||||
"musl/src/math/sinl.c",
|
||||
"musl/src/math/__tan.c",
|
||||
"musl/src/math/__tandf.c",
|
||||
"musl/src/math/tanh.c",
|
||||
"musl/src/math/tanhf.c",
|
||||
"musl/src/math/tanhl.c",
|
||||
"musl/src/math/__tanl.c",
|
||||
"musl/src/math/tanl.c",
|
||||
|
||||
@@ -665,10 +665,8 @@ const libc_top_half_src_files = [_][]const u8{
|
||||
"musl/src/locale/wcscoll.c",
|
||||
"musl/src/locale/wcsxfrm.c",
|
||||
"musl/src/math/acosh.c",
|
||||
"musl/src/math/acoshf.c",
|
||||
"musl/src/math/acoshl.c",
|
||||
"musl/src/math/acosl.c",
|
||||
"musl/src/math/asin.c",
|
||||
"musl/src/math/asinf.c",
|
||||
"musl/src/math/asinh.c",
|
||||
"musl/src/math/asinhf.c",
|
||||
@@ -757,9 +755,6 @@ const libc_top_half_src_files = [_][]const u8{
|
||||
"musl/src/math/__math_uflowf.c",
|
||||
"musl/src/math/__math_xflow.c",
|
||||
"musl/src/math/__math_xflowf.c",
|
||||
"musl/src/math/modf.c",
|
||||
"musl/src/math/modff.c",
|
||||
"musl/src/math/modfl.c",
|
||||
"musl/src/math/nearbyintl.c",
|
||||
"musl/src/math/nextafter.c",
|
||||
"musl/src/math/nextafterf.c",
|
||||
@@ -798,8 +793,6 @@ const libc_top_half_src_files = [_][]const u8{
|
||||
"musl/src/math/sinl.c",
|
||||
"musl/src/math/__tan.c",
|
||||
"musl/src/math/__tandf.c",
|
||||
"musl/src/math/tanh.c",
|
||||
"musl/src/math/tanhf.c",
|
||||
"musl/src/math/tanhl.c",
|
||||
"musl/src/math/__tanl.c",
|
||||
"musl/src/math/tanl.c",
|
||||
|
||||
Reference in New Issue
Block a user