From e995e289db40d82bfe388ccfb1088e518cce3b7e Mon Sep 17 00:00:00 2001 From: Taylor Wampler Date: Mon, 23 Dec 2024 19:21:42 -0600 Subject: [PATCH] math::nolibc: acos and asin There are now nolibc definitions for the inverse cosine and inverse sine. More test points were added for acos, asin, and atan in the math_tests module. This was done becuase the nolibc inverse trigonometric functions have various branching conditions depending on the provided input value. Several branches in these functions were neglected. --- lib/std/math/math_nolibc/acos.c3 | 146 +++++++++++++++++++++++++++++++ lib/std/math/math_nolibc/asin.c3 | 136 ++++++++++++++++++++++++++++ test/unit/stdlib/math/math.c3 | 33 ++++++- 3 files changed, 312 insertions(+), 3 deletions(-) create mode 100644 lib/std/math/math_nolibc/acos.c3 create mode 100644 lib/std/math/math_nolibc/asin.c3 diff --git a/lib/std/math/math_nolibc/acos.c3 b/lib/std/math/math_nolibc/acos.c3 new file mode 100644 index 000000000..e638186a2 --- /dev/null +++ b/lib/std/math/math_nolibc/acos.c3 @@ -0,0 +1,146 @@ +module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH)); + +/* origin: FreeBSD /usr/src/lib/msun/src/e_acos.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. + * ==================================================== + */ + +const PIO2_HI @local = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ +const PIO2_LO @local = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */ +const PS0 @local = 1.66666666666666657415e-01; /* 0x3FC55555, 0x55555555 */ +const PS1 @local = -3.25565818622400915405e-01; /* 0xBFD4D612, 0x03EB6F7D */ +const PS2 @local = 2.01212532134862925881e-01; /* 0x3FC9C155, 0x0E884455 */ +const PS3 @local = -4.00555345006794114027e-02; /* 0xBFA48228, 0xB5688F3B */ +const PS4 @local = 7.91534994289814532176e-04; /* 0x3F49EFE0, 0x7501B288 */ +const PS5 @local = 3.47933107596021167570e-05; /* 0x3F023DE1, 0x0DFDF709 */ +const QS1 @local = -2.40339491173441421878e+00; /* 0xC0033A27, 0x1C8A2D4B */ +const QS2 @local = 2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */ +const QS3 @local = -6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */ +const QS4 @local = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + +fn double _r(double z) @local +{ + double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5))))); + double q = 1.0 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); + return p / q; +} + +fn double _acos(double x) @weak @extern("acos") @nostrip +{ + uint hx = x.high_word(); + uint ix = hx & 0x7fffffff; + switch + { + /* |x| >= 1 or nan */ + case ix >= 0x3ff00000: + uint lx = x.low_word(); + if ((ix - 0x3ff00000 | lx) == 0) + { + /* acos(1)=0, acos(-1)=pi */ + if (hx >> 31) return 2. * PIO2_HI + 0x1p-120f; + return 0.; + } + return double.nan; + /* |x| < 0.5 */ + case ix < 0x3fe00000: + /* |x| < 2**-57 */ + if (ix <= 0x3c600000) return PIO2_HI + 0x1p-120f; + return PIO2_HI - (x - (PIO2_LO - x * _r(x * x))); + /* x < -0.5 */ + case (hx >> 31) != 0: + double z = (1. + x) * 0.5; + double s = math::sqrt(z); + double w = _r(z) * s - PIO2_LO; + return 2. * (PIO2_HI - (s + w)); + /* x > 0.5 */ + default: + double z = (1. - x) * 0.5; + double s = math::sqrt(z); + double df @noinit; + { + ulong rep = bitcast(s, ulong); + rep &= 0xffffffff00000000; + df = bitcast(rep, double); + } + double c = (z - df * df) / (s + df); + double w = _r(z) * s + c; + return 2. * (df + w); + } +} + +/* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.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. + * ==================================================== + */ + +const float PIO2_HI_F @local = 1.5707962513e+00; /* 0x3fc90fda */ +const float PIO2_LO_F @local = 7.5497894159e-08; /* 0x33a22168 */ +const float PS0_F @local = 1.6666586697e-01; +const float PS1_F @local = -4.2743422091e-02; +const float PS2_F @local = -8.6563630030e-03; +const float QS1_F @local = -7.0662963390e-01; + +fn float _r_f(float z) @local +{ + float p = z * ( PS0_F + z * (PS1_F + z * PS2_F)); + float q = 1.0f + z * QS1_F; + return p / q; +} + +fn float _acosf(float x) @weak @extern("acosf") @nostrip +{ + uint hx = bitcast(x, uint); + uint ix = hx & 0x7fffffff; + switch + { + /* |x| >= 1 or nan */ + case ix >= 0x3f800000: + if (ix == 0x3f800000) + { + if (hx >> 31) return 2.f * PIO2_HI_F + 0x1p-120f; + return 0; + } + return float.nan; + /* |x| < 0.5 */ + case ix < 0x3f000000: + /* |x| < 2**-26 */ + if (ix <= 0x32800000) return PIO2_HI_F + 0x1p-120f; + return PIO2_HI_F - (x - (PIO2_LO_F - x * _r_f(x * x))); + /* x < -0.5 */ + case (hx >> 31) != 0: + float z = (1.f + x) * 0.5f; + float s = math::sqrt(z); + float w = _r_f(z) * s - PIO2_LO_F; + return 2.f * (PIO2_HI_F - (s + w)); + /* x > 0.5 */ + default: + float z = (1.f - x) * 0.5f; + float s = math::sqrt(z); + float df @noinit; + { + uint rep = bitcast(s, uint); + rep &= 0xfffff000; + df = bitcast(rep, float); + } + float c = (z - df * df) / (s + df); + float w = _r_f(z) * s + c; + return 2.f * (df + w); + } +} \ No newline at end of file diff --git a/lib/std/math/math_nolibc/asin.c3 b/lib/std/math/math_nolibc/asin.c3 new file mode 100644 index 000000000..e0a9f469b --- /dev/null +++ b/lib/std/math/math_nolibc/asin.c3 @@ -0,0 +1,136 @@ +module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH)); + +/* 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. + * ==================================================== + */ + +const PIO2_HI @local = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ +const PIO2_LO @local = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */ +const PS0 @local = 1.66666666666666657415e-01; /* 0x3FC55555, 0x55555555 */ +const PS1 @local = -3.25565818622400915405e-01; /* 0xBFD4D612, 0x03EB6F7D */ +const PS2 @local = 2.01212532134862925881e-01; /* 0x3FC9C155, 0x0E884455 */ +const PS3 @local = -4.00555345006794114027e-02; /* 0xBFA48228, 0xB5688F3B */ +const PS4 @local = 7.91534994289814532176e-04; /* 0x3F49EFE0, 0x7501B288 */ +const PS5 @local = 3.47933107596021167570e-05; /* 0x3F023DE1, 0x0DFDF709 */ +const QS1 @local = -2.40339491173441421878e+00; /* 0xC0033A27, 0x1C8A2D4B */ +const QS2 @local = 2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */ +const QS3 @local = -6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */ +const QS4 @local = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + +fn double _r(double z) @local +{ + double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5))))); + double q = 1.0 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); + return p / q; +} + +fn double _asin(double x) @weak @extern("asin") @nostrip +{ + uint hx = x.high_word(); + uint ix = hx & 0x7fffffff; + switch + { + /* |x| >= 1 or nan */ + case ix >= 0x3ff00000: + uint lx = x.low_word(); + if ((ix-0x3ff00000 | lx) == 0) + { + /* asin(1) = +-pi/2 with inexact */ + return x * PIO2_HI + 0x1p-120f; + } + return double.nan; + /* |x| < 0.5 */ + case 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 */ + default: + double z = (1. - math::abs(x)) * 0.5; + double s = math::sqrt(z); + double r = _r(z); + /* if |x| > 0.975 */ + if (ix >= 0x3fef3333) + { + x = PIO2_HI - (2. * (s + s * r) - PIO2_LO); + } + else + { + double f @noinit; + { + ulong rep = bitcast(s, ulong); + rep &= 0xffffffff00000000; + f = bitcast(rep, double); + } + double 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 != 0) return -x; + return x; + } +} + +/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.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. + * ==================================================== + */ + +const PIO2 @local = 1.570796326794896558e+00; +const float PS0_F @local = 1.6666586697e-01; +const float PS1_F @local = -4.2743422091e-02; +const float PS2_F @local = -8.6563630030e-03; +const float QS1_F @local = -7.0662963390e-01; + +fn float _r_f(float z) @local +{ + float p = z * ( PS0_F + z * (PS1_F + z * PS2_F)); + float q = 1.0f + z * QS1_F; + return p / q; +} + +fn float _asinf(float x) @weak @extern("asinf") @nostrip +{ + uint hx = bitcast(x, uint); + uint ix = hx & 0x7fffffff; + switch + { + /* |x| >= 1 */ + case ix >= 0x3f800000: + if (ix == 0x3f800000) + { + /* asin(+-1) = +-pi/2 with inexact */ + return x * (float)PIO2 + 0x1p-120f; + } + return float.nan; + /* |x| < 0.5 */ + case ix < 0x3f000000: + /* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */ + if (ix < 0x39800000 && ix >= 0x00800000) return x; + return x + x * _r_f(x * x); + /* 1 > |x| >= 0.5 */ + default: + float z = (1.f - math::abs(x)) * 0.5f; + float s = math::sqrt(z); + x = (float)PIO2 - 2.f * (s + s * _r_f(z)); + if (hx >> 31) return -x; + return x; + } +} \ No newline at end of file diff --git a/test/unit/stdlib/math/math.c3 b/test/unit/stdlib/math/math.c3 index d0b8ffe78..3ed416805 100644 --- a/test/unit/stdlib/math/math.c3 +++ b/test/unit/stdlib/math/math.c3 @@ -20,6 +20,8 @@ fn void! test_acos() @test { int [<5>] in = { 231, -231, 1, 0, -1 }; double [<3>] out = { 0., math::PI_2, math::PI }; + double [<6>] in2 = { 0.9, 0.6, 0.1, -0.1, -0.6, -0.9 }; + double [<6>] out2 = { 0.45102681179626236, 0.9272952180016123, 1.4706289056333368, 1.6709637479564565, 2.214297435588181, 2.6905658417935308 }; assert(@typeis(math::acos(in[0]), double)); assert(@typeis(math::acos((float)in[0]), float)); assert(@typeis(math::acos((double)in[0]), double)); @@ -42,6 +44,13 @@ fn void! test_acos() @test x = math::acos((double)in[i]); assert(math::is_approx_rel(x, out[ii], 1e-12), "acos(%f)=%f is not equal to %f", in[i], x, out[ii]); } + for (int i = 0; i < 6; i++) + { + float f = math::acos((float)in2[i]); + assert(math::is_approx(f, (float)out2[i], 1e-6), "acos(%f)=%f is not equal to %f", (float)in2[i], f, (float)out2[i]); + double x = math::acos(in2[i]); + assert(math::is_approx(x, out2[i], 1e-12), "acos(%f)=%f is not equal to %f", in2[i], x, out2[i]); + } } @@ -74,6 +83,8 @@ fn void! test_asin() @test { int [<5>] in = { 231, -231, 1, 0, -1 }; double [<3>] out = { math::PI_2, 0., -math::PI_2 }; + double [<6>] in2 = { 0.98, 0.6, 0.1, -0.1, -0.6, -0.98 }; + double [<6>] out2 = { 1.3704614844717768, 0.6435011087932844, 0.1001674211615598, -0.1001674211615598, -0.6435011087932844, -1.3704614844717768 }; assert(@typeis(math::asin(in[0]), double)); assert(@typeis(math::asin((float)in[0]), float)); assert(@typeis(math::asin((double)in[0]), double)); @@ -93,6 +104,13 @@ fn void! test_asin() @test x = math::asin((double)in[i]); assert(math::is_approx_rel(x, out[ii], 1e-12), "asin(%f)=%f is not equal to %f", in[i], x, out[ii]); } + for (int i = 0; i < 6; i++) + { + float f = math::asin((float)in2[i]); + assert(math::is_approx(f, (float)out2[i], 1e-6), "asin(%f)=%f is not equal to %f", (float)in2[i], f, (float)out2[i]); + double x = math::asin(in2[i]); + assert(math::is_approx(x, out2[i], 1e-12), "asin(%f)=%f is not equal to %f", in2[i], x, out2[i]); + } } fn void! test_asinh() @test @@ -115,12 +133,14 @@ fn void! test_asinh() @test fn void! test_atan() @test { - int [<5>] in = { 231, 1, 0, -1, -231 }; - double [<5>] out = { 1.5664673495078372, math::PI_4, 0., -math::PI_4, -1.5664673495078372 }; + int [<9>] in = { 231, 3, 2, 1, 0, -1, -2, -3, -231 }; + double [<9>] out = { 1.5664673495078372, 1.2490457723982544, 1.1071487177940904, math::PI_4, 0., -math::PI_4, -1.1071487177940904, -1.2490457723982544, -1.5664673495078372 }; + double [<6>] in2 = { 0.6, 0.4, 0.1, -0.1, -0.4, -0.6 }; + double [<6>] out2 = { 0.5404195002705842, 0.3805063771123649, 0.09966865249116204, -0.09966865249116204, -0.3805063771123649, -0.5404195002705842 }; assert(@typeis(math::atan(in[0]), double)); assert(@typeis(math::atan((float)in[0]), float)); assert(@typeis(math::atan((double)in[0]), double)); - for (int i = 0; i < 5; i++) + for (int i = 0; i < 9; i++) { double x = math::atan(in[i]); assert(math::is_approx_rel(x, out[i], 1e-12), "atan(%d)=%f is not equal to %f", in[i], x, out[i]); @@ -129,6 +149,13 @@ fn void! test_atan() @test x = math::atan((double)in[i]); assert(math::is_approx_rel(x, out[i], 1e-12), "atan(%f)=%f is not equal to %f", in[i], x, out[i]); } + for (int i = 0; i < 6; i++) + { + float f = math::atan((float)in2[i]); + assert(math::is_approx(f, (float)out2[i], 1e-6), "atan(%f)=%f is not equal to %f", (float)in2[i], f, (float)out2[i]); + double x = math::atan(in2[i]); + assert(math::is_approx(x, out2[i], 1e-12), "atan(%f)=%f is not equal to %f", in2[i], x, out2[i]); + } } fn void! test_atanh() @test