From 53bada2a1e60a63aff6a4c42c092fce352608a78 Mon Sep 17 00:00:00 2001 From: Taylor W Date: Sat, 28 Dec 2024 14:13:44 -0600 Subject: [PATCH] math::nolibc: atanh (#1730) * math::nolibc: log1p * math::no_libc: atanh Added atanh nolibc definition and more test points in the math_tests module. --- lib/std/math/math_nolibc/atanh.c3 | 99 +++++++++++++ lib/std/math/math_nolibc/log1p.c3 | 228 ++++++++++++++++++++++++++++++ test/unit/stdlib/math/math.c3 | 8 +- 3 files changed, 331 insertions(+), 4 deletions(-) create mode 100644 lib/std/math/math_nolibc/atanh.c3 create mode 100644 lib/std/math/math_nolibc/log1p.c3 diff --git a/lib/std/math/math_nolibc/atanh.c3 b/lib/std/math/math_nolibc/atanh.c3 new file mode 100644 index 000000000..152e218c0 --- /dev/null +++ b/lib/std/math/math_nolibc/atanh.c3 @@ -0,0 +1,99 @@ +module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH)); + +/* origin: FreeBSD usr/src/lib/msun/src/e_atanh.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. + * ==================================================== + */ + +fn double _atanh(double x) @weak @extern("atanh") @nostrip +{ + double t @noinit; + uint hx = x.high_word(); + uint ix = hx & 0x7fffffff; + uint sign = hx >> 31; + switch + { + /* |x| >= 1 or nan */ + case ix >= 0x3ff00000: + uint lx = x.low_word(); + if ((ix - 0x3ff00000 | lx) == 0) + { + return sign ? -double.inf : double.inf; + } + return double.nan; + /* x<2**-28 */ + case ix < 0x3e300000 && (1e300 + x) > 0.: + return x; + } + { + ulong rep = bitcast(x, ulong); + rep = (rep & 0x00000000ffffffff) | (ulong)ix << 32; + x = bitcast(rep, double); + } + /* |x| < 0.5 */ + if (ix < 0x3fe00000) + { + t = x + x; + t = 0.5 * _log1p(t + t * x / (1. - x)); + } + else + { + t = 0.5 * _log1p((x + x) / (1. - x)); + } + return sign ? -t : t; +} + +/* origin: FreeBSD usr/src/lib/msun/src/e_atanhf.c */ +/* + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +fn float _atanhf(float x) @weak @extern("atanhf") @nostrip +{ + float t @noinit; + uint hx = bitcast(x, uint); + uint ix = hx & 0x7fffffff; + uint sign = hx >> 31; + switch + { + /* |x| >= 1 or nan */ + case ix >= 0x3f800000: + if (ix == 0x3f800000) + { + return sign ? -float.inf : float.inf; + } + return float.nan; + /* x<2**-28 */ + case ix < 0x31800000 && (1e30 + x) > 0.f: + return x; + } + x = bitcast(ix, float); + /* |x| < 0.5 */ + if (ix < 0x3f000000) + { + t = x + x; + t = 0.5f * _log1pf(t + t * x / (1.f - x)); + } + else + { + t = 0.5f * _log1pf((x + x) / (1.f - x)); + } + return sign ? -t : t; +} diff --git a/lib/std/math/math_nolibc/log1p.c3 b/lib/std/math/math_nolibc/log1p.c3 new file mode 100644 index 000000000..9fe0984a5 --- /dev/null +++ b/lib/std/math/math_nolibc/log1p.c3 @@ -0,0 +1,228 @@ +module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH)); + +/* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.c */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* origin: musl libc /src/math/log1p.c */ +/* + * ==================================================== + * Copyright (c) 2005-2020 Rich Felker, et al. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * ==================================================== + */ + +const LN2_HI @local = 6.93147180369123816490e-01; /* 3fe62e42 fee00000 */ +const LN2_LO @local = 1.90821492927058770002e-10; /* 3dea39ef 35793c76 */ +const LG1 @local = 6.666666666666735130e-01; /* 3FE55555 55555593 */ +const LG2 @local = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */ +const LG3 @local = 2.857142874366239149e-01; /* 3FD24924 94229359 */ +const LG4 @local = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */ +const LG5 @local = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */ +const LG6 @local = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */ +const LG7 @local = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + +fn double _log1p(double x) @weak @extern("log1p") @nostrip +{ + uint hx = x.high_word(); + int k = 1; + double c @noinit; + double f @noinit; + switch + { + /* 1+x < sqrt(2)+ */ + case hx < 0x3fda827a || hx >> 31 != 0: + switch + { + /* x <= -1.0 */ + case hx >= 0xbff00000: + if (x == -1) return -double.inf; + return double.nan; + /* |x| < 2**-53 */ + case hx << 1 < 0x3ca00000 << 1: + /* underflow if subnormal */ + if ((hx & 0x7ff00000) == 0) + { + (float)@volatile_load(x); + } + return x; + /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ + case hx <= 0xbfd2bec4: + k = 0; + c = 0; + f = x; + } + case hx >= 0x7ff00000: + return x; + } + if (k) + { + double u = 1 + x; + uint hu = u.high_word(); + hu += 0x3ff00000 - 0x3fe6a09e; + k = (int)(hu >> 20) - 0x3ff; + /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ + if (k < 54) + { + c = (k >= 2) ? 1. - (u - x) : x - (u - 1.); + c /= u; + } + else + { + c = 0; + } + /* reduce u into [sqrt(2)/2, sqrt(2)] */ + hu = (hu & 0x000fffff) + 0x3fe6a09e; + u = bitcast(((ulong)hu << 32) | (bitcast(u, ulong) & 0xffffffff) , double); + f = u - 1.; + } + double hfsq = 0.5 * f * f; + double s = f / (2. + f); + double z = s * s; + double w = z * z; + double t1 = w * (LG2 + w * (LG4 + w * LG6)); + double t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7))); + double r = t1 + t2; + double dk = k; + return s * (hfsq + r) + (dk * LN2_LO + c) - hfsq + f + dk * LN2_HI; +} + +/* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.c */ +/* + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* origin: musl libc /src/math/log1pf.c */ +/* + * ==================================================== + * Copyright (c) 2005-2020 Rich Felker, et al. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * ==================================================== + */ + +const float LN2_HI_F @local = 6.9313812256e-01; /* 0x3f317180 */ +const float LN2_LO_F @local = 9.0580006145e-06; /* 0x3717f7d1 */ +/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ +const float LG1_F @local = 0xaaaaaa.0p-24; /* 0.66666662693 */ +const float LG2_F @local = 0xccce13.0p-25; /* 0.40000972152 */ +const float LG3_F @local = 0x91e9ee.0p-25; /* 0.28498786688 */ +const float LG4_F @local = 0xf89e26.0p-26; /* 0.24279078841 */ + +fn float _log1pf(float x) @weak @extern("log1pf") @nostrip +{ + uint ix = x.word(); + int k = 1; + float c @noinit; + float f @noinit; + switch + { + /* 1+x < sqrt(2)+ */ + case ix < 0x3ed413d0 || ix >> 31 != 0: + switch + { + /* x <= -1.0 */ + case ix >= 0xbf800000: + if (x == -1) return -float.inf; + return float.nan; + /* |x| < 2**-24 */ + case ix << 1 < 0x33800000 << 1: + /* underflow if subnormal */ + if ((ix & 0x7f800000) == 0) + { + float v = @volatile_load(x); + v = v * v; + } + return x; + /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ + case ix <= 0xbe95f619: + k = 0; + c = 0; + f = x; + } + case ix >= 0x7f800000: + return x; + } + if (k) + { + float u = 1 + x; + uint iu = u.word(); + iu += 0x3f800000 - 0x3f3504f3; + k = (int)(iu >> 23) - 0x7f; + /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ + if (k < 25) + { + c = (k >= 2) ? 1.f - (u - x) : x - (u - 1.f); + c /= u; + } + else + { + c = 0; + } + /* reduce u into [sqrt(2)/2, sqrt(2)] */ + iu = (iu & 0x007fffff) + 0x3f3504f3; + f = bitcast(iu, float) - 1.f; + } + float hfsq = 0.5f * f * f; + float s = f / (2.f + f); + float z = s * s; + float w = z * z; + float t1 = w * (LG2_F + w * LG4_F); + float t2 = z * (LG1_F + w * LG3_F); + float r = t1 + t2; + float dk = k; + return s * (hfsq + r) + (dk * LN2_LO_F + c) - hfsq + f + dk * LN2_HI_F; +} \ No newline at end of file diff --git a/test/unit/stdlib/math/math.c3 b/test/unit/stdlib/math/math.c3 index 3ed416805..032bc551e 100644 --- a/test/unit/stdlib/math/math.c3 +++ b/test/unit/stdlib/math/math.c3 @@ -161,14 +161,14 @@ fn void! test_atan() @test fn void! test_atanh() @test { int [<4>] in = { 231, -231, 1, -1 }; - double [<2>] in2 = { 0.5, -0.5 }; - double [<2>] out = { 0.5493061443340548, -0.5493061443340548 }; + double [<6>] in2 = {0.8, 0.5, 0.3, -0.3, -0.5, -0.8 }; + double [<6>] out = { 1.0986122886681098, 0.5493061443340548, 0.30951960420311175, -0.30951960420311175, -0.5493061443340548, -1.0986122886681098 }; assert(@typeis(math::atanh(in[0]), double)); assert(@typeis(math::atanh((float)in[0]), float)); assert(@typeis(math::atanh((double)in[0]), double)); for (int i = 0; i < 2; i++) { - assert(math::is_nan(math::atanh(in[i])), "atanh(%d)=%f is not nan", in[i]); + assert(math::is_nan(math::atanh(in[i])), "atanh(%d) is not nan", in[i]); assert(math::is_nan(math::atanh((float)in[i])), "atanh(%f) is not nan", in[i]); assert(math::is_nan(math::atanh((double)in[i])), "atanh(%f) is not nan", in[i]); } @@ -181,7 +181,7 @@ fn void! test_atanh() @test assert(math::atanh(0) == 0., "atanh(%d) is not equal to %f", 0, 0.); assert(math::atanh(0.f) == 0.f, "atanh(%f) is not equal to %f", 0.f, 0.f); assert(math::atanh(0.) == 0., "atanh(%f) is not equal to %f", 0., 0.); - for (int i = 0; i < 2; i++) + for (int i = 0; i < 6; i++) { float f = math::atanh((float)in2[i]); assert(math::is_approx(f, (float)out[i], 1e-6), "atanh(%f)=%f is not equal to %f", in2[i], f, out[i]);