From 3027d67037c0488ec6c125b90b76e687422399fa Mon Sep 17 00:00:00 2001 From: Justine Tunney Date: Tue, 12 Jul 2022 15:49:11 -0700 Subject: [PATCH] Import more Musl math --- libc/tinymath/{acosf.S => acosf-tiny.S} | 3 + libc/tinymath/acosf.c | 111 ++++++++++++ libc/tinymath/{asinf.S => asinf-tiny.S} | 3 + libc/tinymath/asinf.c | 102 +++++++++++ libc/tinymath/{atan.S => atan-tiny.S} | 3 + libc/tinymath/atan.c | 156 ++++++++++++++++ libc/tinymath/{atanf.S => atanf-tiny.S} | 3 + libc/tinymath/atanf.c | 134 ++++++++++++++ libc/tinymath/cos.c | 6 +- libc/tinymath/cosf.S | 28 --- libc/tinymath/cosf.c | 121 +++++++++++++ libc/tinymath/kernel.internal.h | 7 +- libc/tinymath/ldexp.S | 41 ----- libc/tinymath/ldexp.c | 70 ++++++++ libc/tinymath/ldexpf.S | 41 ----- libc/tinymath/ldexpf.c | 68 +++++++ libc/tinymath/ldexpl.S | 1 - libc/tinymath/{log1p.S => log1p-tiny.S} | 3 + libc/tinymath/log1p.c | 164 +++++++++++++++++ libc/tinymath/pow.c | 1 + libc/tinymath/{powf.S => powf-tiny.S} | 3 + libc/tinymath/powf.c | 226 ++++++++++++++++++++++++ libc/tinymath/powf_data.c | 67 +++++++ libc/tinymath/powf_data.internal.h | 25 +++ libc/tinymath/powl.c | 1 + libc/tinymath/rempio2.c | 2 +- libc/tinymath/rempio2f.c | 124 +++++++++++++ libc/tinymath/scalbln.c | 26 +++ libc/tinymath/scalblnf.c | 26 +++ libc/tinymath/scalbn.c | 26 +++ libc/tinymath/scalbnf.c | 26 +++ libc/tinymath/sin.c | 6 +- libc/tinymath/sincos.c | 4 + libc/tinymath/sinf.S | 28 --- libc/tinymath/sinf.c | 119 +++++++++++++ libc/tinymath/tan.c | 3 + test/libc/tinymath/acos_test.c | 13 +- test/libc/tinymath/asin_test.c | 13 +- test/libc/tinymath/atan_test.c | 10 ++ test/libc/tinymath/cos_test.c | 10 ++ test/libc/tinymath/ldexp_test.c | 11 +- test/libc/tinymath/log1p_test.c | 10 ++ test/libc/tinymath/powl_test.c | 2 +- test/libc/tinymath/sin_test.c | 11 +- test/libc/tinymath/sincos_test.c | 40 ++++- test/libc/tinymath/tan_test.c | 10 ++ test/tool/net/ljson_test.lua | 17 +- tool/net/ljson.c | 4 +- 48 files changed, 1749 insertions(+), 180 deletions(-) rename libc/tinymath/{acosf.S => acosf-tiny.S} (98%) create mode 100644 libc/tinymath/acosf.c rename libc/tinymath/{asinf.S => asinf-tiny.S} (98%) create mode 100644 libc/tinymath/asinf.c rename libc/tinymath/{atan.S => atan-tiny.S} (98%) create mode 100644 libc/tinymath/atan.c rename libc/tinymath/{atanf.S => atanf-tiny.S} (98%) create mode 100644 libc/tinymath/atanf.c delete mode 100644 libc/tinymath/cosf.S create mode 100644 libc/tinymath/cosf.c delete mode 100644 libc/tinymath/ldexp.S create mode 100644 libc/tinymath/ldexp.c delete mode 100644 libc/tinymath/ldexpf.S create mode 100644 libc/tinymath/ldexpf.c rename libc/tinymath/{log1p.S => log1p-tiny.S} (98%) create mode 100644 libc/tinymath/log1p.c rename libc/tinymath/{powf.S => powf-tiny.S} (98%) create mode 100644 libc/tinymath/powf.c create mode 100644 libc/tinymath/powf_data.c create mode 100644 libc/tinymath/powf_data.internal.h create mode 100644 libc/tinymath/rempio2f.c create mode 100644 libc/tinymath/scalbln.c create mode 100644 libc/tinymath/scalblnf.c create mode 100644 libc/tinymath/scalbn.c create mode 100644 libc/tinymath/scalbnf.c delete mode 100644 libc/tinymath/sinf.S create mode 100644 libc/tinymath/sinf.c diff --git a/libc/tinymath/acosf.S b/libc/tinymath/acosf-tiny.S similarity index 98% rename from libc/tinymath/acosf.S rename to libc/tinymath/acosf-tiny.S index 144c6743d..deb9f6f1a 100644 --- a/libc/tinymath/acosf.S +++ b/libc/tinymath/acosf-tiny.S @@ -17,6 +17,7 @@ │ PERFORMANCE OF THIS SOFTWARE. │ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/macros.internal.h" +#ifdef TINY // Returns arc cosine of 𝑥. // @@ -25,3 +26,5 @@ acosf: ezlea acosl,ax jmp _f2ld2 .endfn acosf,globl + +#endif /* TINY */ diff --git a/libc/tinymath/acosf.c b/libc/tinymath/acosf.c new file mode 100644 index 000000000..b80eb1431 --- /dev/null +++ b/libc/tinymath/acosf.c @@ -0,0 +1,111 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/complex.internal.h" +#ifndef TINY + +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* 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. + * ==================================================== + */ + +static const float +pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */ +pio2_lo = 7.5497894159e-08, /* 0x33a22168 */ +pS0 = 1.6666586697e-01, +pS1 = -4.2743422091e-02, +pS2 = -8.6563630030e-03, +qS1 = -7.0662963390e-01; + +static float R(float z) +{ + float_t p, q; + p = z*(pS0+z*(pS1+z*pS2)); + q = 1.0f+z*qS1; + return p/q; +} + +/** + * Returns arc cosine of 𝑥. + */ +float acosf(float x) +{ + float z,w,s,c,df; + uint32_t hx,ix; + + GET_FLOAT_WORD(hx, x); + ix = hx & 0x7fffffff; + /* |x| >= 1 or nan */ + if (ix >= 0x3f800000) { + if (ix == 0x3f800000) { + if (hx >> 31) + return 2*pio2_hi + 0x1p-120f; + return 0; + } + return 0/(x-x); + } + /* |x| < 0.5 */ + if (ix < 0x3f000000) { + if (ix <= 0x32800000) /* |x| < 2**-26 */ + return pio2_hi + 0x1p-120f; + return pio2_hi - (x - (pio2_lo-x*R(x*x))); + } + /* x < -0.5 */ + if (hx >> 31) { + z = (1+x)*0.5f; + s = sqrtf(z); + w = R(z)*s-pio2_lo; + return 2*(pio2_hi - (s+w)); + } + /* x > 0.5 */ + z = (1-x)*0.5f; + s = sqrtf(z); + GET_FLOAT_WORD(hx,s); + SET_FLOAT_WORD(df,hx&0xfffff000); + c = (z-df*df)/(s+df); + w = R(z)*s+c; + return 2*(df+w); +} + +#endif /* TINY */ diff --git a/libc/tinymath/asinf.S b/libc/tinymath/asinf-tiny.S similarity index 98% rename from libc/tinymath/asinf.S rename to libc/tinymath/asinf-tiny.S index 3f65532f8..4605276f6 100644 --- a/libc/tinymath/asinf.S +++ b/libc/tinymath/asinf-tiny.S @@ -17,6 +17,7 @@ │ PERFORMANCE OF THIS SOFTWARE. │ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/macros.internal.h" +#ifdef TINY // Returns arc sine of 𝑥. // @@ -25,3 +26,5 @@ asinf: ezlea asinl,ax jmp _f2ld2 .endfn asinf,globl + +#endif /* TINY */ diff --git a/libc/tinymath/asinf.c b/libc/tinymath/asinf.c new file mode 100644 index 000000000..48627b488 --- /dev/null +++ b/libc/tinymath/asinf.c @@ -0,0 +1,102 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/complex.internal.h" +#ifndef TINY + +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* 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. + * ==================================================== + */ + +static const double +pio2 = 1.570796326794896558e+00; + +static const float +/* coefficients for R(x^2) */ +pS0 = 1.6666586697e-01, +pS1 = -4.2743422091e-02, +pS2 = -8.6563630030e-03, +qS1 = -7.0662963390e-01; + +static float R(float z) +{ + float_t p, q; + p = z*(pS0+z*(pS1+z*pS2)); + q = 1.0f+z*qS1; + return p/q; +} + +/** + * Returns arc sine of 𝑥. + */ +float asinf(float x) +{ + double s; + float z; + uint32_t hx,ix; + + GET_FLOAT_WORD(hx, x); + ix = hx & 0x7fffffff; + if (ix >= 0x3f800000) { /* |x| >= 1 */ + if (ix == 0x3f800000) /* |x| == 1 */ + return x*pio2 + 0x1p-120f; /* asin(+-1) = +-pi/2 with inexact */ + return 0/(x-x); /* asin(|x|>1) is NaN */ + } + if (ix < 0x3f000000) { /* |x| < 0.5 */ + /* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */ + if (ix < 0x39800000 && ix >= 0x00800000) + return x; + return x + x*R(x*x); + } + /* 1 > |x| >= 0.5 */ + z = (1 - fabsf(x))*0.5f; + s = sqrt(z); + x = pio2 - 2*(s+s*R(z)); + if (hx >> 31) + return -x; + return x; +} + +#endif /* TINY */ diff --git a/libc/tinymath/atan.S b/libc/tinymath/atan-tiny.S similarity index 98% rename from libc/tinymath/atan.S rename to libc/tinymath/atan-tiny.S index d50db1438..f622e1cca 100644 --- a/libc/tinymath/atan.S +++ b/libc/tinymath/atan-tiny.S @@ -17,6 +17,7 @@ │ PERFORMANCE OF THIS SOFTWARE. │ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/macros.internal.h" +#ifdef TINY // Returns arc tangent of 𝑥. // @@ -25,3 +26,5 @@ atan: ezlea atanl,ax jmp _d2ld2 .endfn atan,globl + +#endif /* TINY */ diff --git a/libc/tinymath/atan.c b/libc/tinymath/atan.c new file mode 100644 index 000000000..f630eb6e6 --- /dev/null +++ b/libc/tinymath/atan.c @@ -0,0 +1,156 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/complex.internal.h" +#ifndef TINY + +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* origin: FreeBSD /usr/src/lib/msun/src/s_atan.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. + * ==================================================== + */ +/* atan(x) + * Method + * 1. Reduce x to positive by atan(x) = -atan(-x). + * 2. According to the integer k=4t+0.25 chopped, t=x, the argument + * is further reduced to one of the following intervals and the + * arctangent of t is evaluated by the corresponding formula: + * + * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) + * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) + * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) + * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) + * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +static const double atanhi[] = { + 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ + 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ + 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ + 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ +}; + +static const double atanlo[] = { + 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ + 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ + 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ + 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ +}; + +static const double aT[] = { + 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ + -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ + 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ + -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ + 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ + -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ + 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ + -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ + 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ + -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ + 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ +}; + +/** + * Returns arc tangent of 𝑥. + * @note should take ~18ns + */ +double atan(double x) +{ + double_t w,s1,s2,z; + uint32_t ix,sign; + int id; + + GET_HIGH_WORD(ix, x); + sign = ix >> 31; + ix &= 0x7fffffff; + if (ix >= 0x44100000) { /* if |x| >= 2^66 */ + if (isnan(x)) + return x; + z = atanhi[3] + 0x1p-120f; + return sign ? -z : z; + } + if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ + if (ix < 0x3e400000) { /* |x| < 2^-27 */ + if (ix < 0x00100000) + /* raise underflow for subnormal x */ + FORCE_EVAL((float)x); + return x; + } + id = -1; + } else { + x = fabs(x); + if (ix < 0x3ff30000) { /* |x| < 1.1875 */ + if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */ + id = 0; + x = (2.0*x-1.0)/(2.0+x); + } else { /* 11/16 <= |x| < 19/16 */ + id = 1; + x = (x-1.0)/(x+1.0); + } + } else { + if (ix < 0x40038000) { /* |x| < 2.4375 */ + id = 2; + x = (x-1.5)/(1.0+1.5*x); + } else { /* 2.4375 <= |x| < 2^66 */ + id = 3; + x = -1.0/x; + } + } + } + /* end of argument reduction */ + z = x*x; + w = z*z; + /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ + s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); + s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); + if (id < 0) + return x - x*(s1+s2); + z = atanhi[id] - (x*(s1+s2) - atanlo[id] - x); + return sign ? -z : z; +} + +#endif /* TINY */ diff --git a/libc/tinymath/atanf.S b/libc/tinymath/atanf-tiny.S similarity index 98% rename from libc/tinymath/atanf.S rename to libc/tinymath/atanf-tiny.S index 174916519..9f0671dde 100644 --- a/libc/tinymath/atanf.S +++ b/libc/tinymath/atanf-tiny.S @@ -17,6 +17,7 @@ │ PERFORMANCE OF THIS SOFTWARE. │ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/macros.internal.h" +#ifdef TINY // Returns arc tangent of 𝑥. // @@ -25,3 +26,5 @@ atanf: ezlea atanl,ax jmp _f2ld2 .endfn atanf,globl + +#endif /* TINY */ diff --git a/libc/tinymath/atanf.c b/libc/tinymath/atanf.c new file mode 100644 index 000000000..889e65385 --- /dev/null +++ b/libc/tinymath/atanf.c @@ -0,0 +1,134 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/complex.internal.h" +#ifndef TINY + +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.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. + * ==================================================== + */ + +static const float atanhi[] = { + 4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */ + 7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */ + 9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */ + 1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */ +}; + +static const float atanlo[] = { + 5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */ + 3.7748947079e-08, /* atan(1.0)lo 0x33222168 */ + 3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */ + 7.5497894159e-08, /* atan(inf)lo 0x33a22168 */ +}; + +static const float aT[] = { + 3.3333328366e-01, + -1.9999158382e-01, + 1.4253635705e-01, + -1.0648017377e-01, + 6.1687607318e-02, +}; + +/** + * Returns arc tangent of 𝑥. + * @note should take ~12ns + */ +float atanf(float x) +{ + float_t w,s1,s2,z; + uint32_t ix,sign; + int id; + + GET_FLOAT_WORD(ix, x); + sign = ix>>31; + ix &= 0x7fffffff; + if (ix >= 0x4c800000) { /* if |x| >= 2**26 */ + if (isnan(x)) + return x; + z = atanhi[3] + 0x1p-120f; + return sign ? -z : z; + } + if (ix < 0x3ee00000) { /* |x| < 0.4375 */ + if (ix < 0x39800000) { /* |x| < 2**-12 */ + if (ix < 0x00800000) + /* raise underflow for subnormal x */ + FORCE_EVAL(x*x); + return x; + } + id = -1; + } else { + x = fabsf(x); + if (ix < 0x3f980000) { /* |x| < 1.1875 */ + if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */ + id = 0; + x = (2.0f*x - 1.0f)/(2.0f + x); + } else { /* 11/16 <= |x| < 19/16 */ + id = 1; + x = (x - 1.0f)/(x + 1.0f); + } + } else { + if (ix < 0x401c0000) { /* |x| < 2.4375 */ + id = 2; + x = (x - 1.5f)/(1.0f + 1.5f*x); + } else { /* 2.4375 <= |x| < 2**26 */ + id = 3; + x = -1.0f/x; + } + } + } + /* end of argument reduction */ + z = x*x; + w = z*z; + /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ + s1 = z*(aT[0]+w*(aT[2]+w*aT[4])); + s2 = w*(aT[1]+w*aT[3]); + if (id < 0) + return x - x*(s1+s2); + z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); + return sign ? -z : z; +} + +#endif /* TINY */ diff --git a/libc/tinymath/cos.c b/libc/tinymath/cos.c index 2c8c967dd..44ec3a96c 100644 --- a/libc/tinymath/cos.c +++ b/libc/tinymath/cos.c @@ -36,8 +36,8 @@ asm(".ident\t\"\\n\\n\ Musl libc (MIT License)\\n\ Copyright 2005-2014 Rich Felker, et. al.\""); asm(".include \"libc/disclaimer.inc\""); - /* clang-format off */ + /* origin: FreeBSD /usr/src/lib/msun/src/s_cos.c */ /* * ==================================================== @@ -83,6 +83,10 @@ asm(".include \"libc/disclaimer.inc\""); #define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i #define gethighw(hi,d) (hi) = asuint64(d) >> 32 +/** + * Returns cosine of 𝑥. + * @note should take ~5ns + */ double cos(double x) { double y[2]; diff --git a/libc/tinymath/cosf.S b/libc/tinymath/cosf.S deleted file mode 100644 index 8dacf7b85..000000000 --- a/libc/tinymath/cosf.S +++ /dev/null @@ -1,28 +0,0 @@ -/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│ -│vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi│ -╞══════════════════════════════════════════════════════════════════════════════╡ -│ Copyright 2020 Justine Alexandra Roberts Tunney │ -│ │ -│ Permission to use, copy, modify, and/or distribute this software for │ -│ any purpose with or without fee is hereby granted, provided that the │ -│ above copyright notice and this permission notice appear in all copies. │ -│ │ -│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ -│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ -│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ -│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ -│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ -│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ -│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ -│ PERFORMANCE OF THIS SOFTWARE. │ -╚─────────────────────────────────────────────────────────────────────────────*/ -#include "libc/macros.internal.h" - -// Returns cosine of 𝑥. -// -// @param 𝑥 is float scalar in low quarter of %xmm0 -// @return float scalar in low quarter of %xmm0 -// @domain -(3π/8) < 𝑥 < 3π/8 for best accuracy -cosf: ezlea cosl,ax - jmp _f2ld2 - .endfn cosf,globl diff --git a/libc/tinymath/cosf.c b/libc/tinymath/cosf.c new file mode 100644 index 000000000..a4e370d16 --- /dev/null +++ b/libc/tinymath/cosf.c @@ -0,0 +1,121 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/complex.internal.h" +#include "libc/tinymath/feval.internal.h" +#include "libc/tinymath/kernel.internal.h" + +asm(".ident\t\"\\n\\n\ +fdlibm (fdlibm license)\\n\ +Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\""); +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* origin: FreeBSD /usr/src/lib/msun/src/s_cosf.c */ +/* + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Optimized by Bruce D. Evans. + */ +/* + * ==================================================== + * 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. + * ==================================================== + */ + +/* Small multiples of pi/2 rounded to double precision. */ +static const double +c1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */ +c2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */ +c3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */ +c4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */ + +/** + * Returns cosine of 𝑥. + * @note should take about 5ns + */ +float cosf(float x) +{ + double y; + uint32_t ix; + unsigned n, sign; + + GET_FLOAT_WORD(ix, x); + sign = ix >> 31; + ix &= 0x7fffffff; + + if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */ + if (ix < 0x39800000) { /* |x| < 2**-12 */ + /* raise inexact if x != 0 */ + FORCE_EVAL(x + 0x1p120f); + return 1.0f; + } + return __cosdf(x); + } + if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */ + if (ix > 0x4016cbe3) /* |x| ~> 3*pi/4 */ + return -__cosdf(sign ? x+c2pio2 : x-c2pio2); + else { + if (sign) + return __sindf(x + c1pio2); + else + return __sindf(c1pio2 - x); + } + } + if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */ + if (ix > 0x40afeddf) /* |x| ~> 7*pi/4 */ + return __cosdf(sign ? x+c4pio2 : x-c4pio2); + else { + if (sign) + return __sindf(-x - c3pio2); + else + return __sindf(x - c3pio2); + } + } + + /* cos(Inf or NaN) is NaN */ + if (ix >= 0x7f800000) + return x-x; + + /* general argument reduction needed */ + n = __rem_pio2f(x,&y); + switch (n&3) { + case 0: return __cosdf(y); + case 1: return __sindf(-y); + case 2: return -__cosdf(y); + default: + return __sindf(y); + } +} diff --git a/libc/tinymath/kernel.internal.h b/libc/tinymath/kernel.internal.h index 565373d0f..9263e4814 100644 --- a/libc/tinymath/kernel.internal.h +++ b/libc/tinymath/kernel.internal.h @@ -5,14 +5,15 @@ COSMOPOLITAN_C_START_ extern int __signgam; +double __cos(double, double) hidden; +double __sin(double, double, int) hidden; +double __tan(double, double, int) hidden; float __cosdf(double) hidden; float __sindf(double) hidden; float __tandf(double, int) hidden; -double __sin(double, double, int) hidden; -double __tan(double, double, int) hidden; -double __cos(double, double) hidden; int __rem_pio2(double, double *) hidden; int __rem_pio2_large(double *, double *, int, int, int) hidden; +int __rem_pio2f(float, double *) hidden; COSMOPOLITAN_C_END_ #endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */ diff --git a/libc/tinymath/ldexp.S b/libc/tinymath/ldexp.S deleted file mode 100644 index 7fd154075..000000000 --- a/libc/tinymath/ldexp.S +++ /dev/null @@ -1,41 +0,0 @@ -/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│ -│vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi│ -╞══════════════════════════════════════════════════════════════════════════════╡ -│ Copyright 2020 Justine Alexandra Roberts Tunney │ -│ │ -│ Permission to use, copy, modify, and/or distribute this software for │ -│ any purpose with or without fee is hereby granted, provided that the │ -│ above copyright notice and this permission notice appear in all copies. │ -│ │ -│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ -│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ -│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ -│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ -│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ -│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ -│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ -│ PERFORMANCE OF THIS SOFTWARE. │ -╚─────────────────────────────────────────────────────────────────────────────*/ -#include "libc/macros.internal.h" - -// Returns 𝑥 × 2ʸ. -// -// @param 𝑥 is double passed in %xmm0 -// @param 𝑦 is exponent via %edi -// @return double in %xmm0 -ldexp: push %rbp - mov %rsp,%rbp - .profilable - push %rdi - fildl (%rsp) - movsd %xmm0,(%rsp) - fldl (%rsp) - fscale - fstp %st(1) - fstpl (%rsp) - movsd (%rsp),%xmm0 - leave - ret - .endfn ldexp,globl - .alias ldexp,scalbn - .alias ldexp,scalbln diff --git a/libc/tinymath/ldexp.c b/libc/tinymath/ldexp.c new file mode 100644 index 000000000..3e4756e6e --- /dev/null +++ b/libc/tinymath/ldexp.c @@ -0,0 +1,70 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/feval.internal.h" +#include "libc/tinymath/kernel.internal.h" + +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/** + * Returns 𝑥 × 2ʸ. + */ +double ldexp(double x, int n) +{ + union {double f; uint64_t i;} u; + double_t y = x; + + if (n > 1023) { + y *= 0x1p1023; + n -= 1023; + if (n > 1023) { + y *= 0x1p1023; + n -= 1023; + if (n > 1023) + n = 1023; + } + } else if (n < -1022) { + /* make sure final n < -53 to avoid double + rounding in the subnormal range */ + y *= 0x1p-1022 * 0x1p53; + n += 1022 - 53; + if (n < -1022) { + y *= 0x1p-1022 * 0x1p53; + n += 1022 - 53; + if (n < -1022) + n = -1022; + } + } + u.i = (uint64_t)(0x3ff+n)<<52; + x = y * u.f; + return x; +} diff --git a/libc/tinymath/ldexpf.S b/libc/tinymath/ldexpf.S deleted file mode 100644 index 485c5e3c7..000000000 --- a/libc/tinymath/ldexpf.S +++ /dev/null @@ -1,41 +0,0 @@ -/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│ -│vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi│ -╞══════════════════════════════════════════════════════════════════════════════╡ -│ Copyright 2020 Justine Alexandra Roberts Tunney │ -│ │ -│ Permission to use, copy, modify, and/or distribute this software for │ -│ any purpose with or without fee is hereby granted, provided that the │ -│ above copyright notice and this permission notice appear in all copies. │ -│ │ -│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ -│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ -│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ -│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ -│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ -│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ -│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ -│ PERFORMANCE OF THIS SOFTWARE. │ -╚─────────────────────────────────────────────────────────────────────────────*/ -#include "libc/macros.internal.h" - -// Returns 𝑥 × 2ʸ. -// -// @param 𝑥 is float passed in %xmm0 -// @param 𝑦 is exponent via %edi -// @return float in %xmm0 -ldexpf: push %rbp - mov %rsp,%rbp - .profilable - push %rdi - fildl (%rsp) - movss %xmm0,(%rsp) - flds (%rsp) - fscale - fstp %st(1) - fstps (%rsp) - movss (%rsp),%xmm0 - leave - ret - .endfn ldexpf,globl - .alias ldexpf,scalbnf - .alias ldexpf,scalblnf diff --git a/libc/tinymath/ldexpf.c b/libc/tinymath/ldexpf.c new file mode 100644 index 000000000..de2e21aca --- /dev/null +++ b/libc/tinymath/ldexpf.c @@ -0,0 +1,68 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/feval.internal.h" +#include "libc/tinymath/kernel.internal.h" + +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/** + * Returns 𝑥 × 2ʸ. + */ +float ldexpf(float x, int n) +{ + union {float f; uint32_t i;} u; + float_t y = x; + + if (n > 127) { + y *= 0x1p127f; + n -= 127; + if (n > 127) { + y *= 0x1p127f; + n -= 127; + if (n > 127) + n = 127; + } + } else if (n < -126) { + y *= 0x1p-126f * 0x1p24f; + n += 126 - 24; + if (n < -126) { + y *= 0x1p-126f * 0x1p24f; + n += 126 - 24; + if (n < -126) + n = -126; + } + } + u.i = (uint32_t)(0x7f+n)<<23; + x = y * u.f; + return x; +} diff --git a/libc/tinymath/ldexpl.S b/libc/tinymath/ldexpl.S index 7aa6ae2f1..ec5991e12 100644 --- a/libc/tinymath/ldexpl.S +++ b/libc/tinymath/ldexpl.S @@ -16,7 +16,6 @@ │ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ │ PERFORMANCE OF THIS SOFTWARE. │ ╚─────────────────────────────────────────────────────────────────────────────*/ -#include "libc/runtime/pc.internal.h" #include "libc/macros.internal.h" // Returns 𝑥 × 2ʸ. diff --git a/libc/tinymath/log1p.S b/libc/tinymath/log1p-tiny.S similarity index 98% rename from libc/tinymath/log1p.S rename to libc/tinymath/log1p-tiny.S index f8f4e8866..4ecb21371 100644 --- a/libc/tinymath/log1p.S +++ b/libc/tinymath/log1p-tiny.S @@ -17,6 +17,7 @@ │ PERFORMANCE OF THIS SOFTWARE. │ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/macros.internal.h" +#ifdef TINY // Returns log(𝟷+𝑥). // @@ -57,3 +58,5 @@ log1p: push %rbp .long 0x95f61998 .long 0x3ffd .long 0 + +#endif /* TINY */ diff --git a/libc/tinymath/log1p.c b/libc/tinymath/log1p.c new file mode 100644 index 000000000..02f47b40c --- /dev/null +++ b/libc/tinymath/log1p.c @@ -0,0 +1,164 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/bits/likely.h" +#include "libc/math.h" +#include "libc/tinymath/internal.h" +#include "libc/tinymath/log_data.internal.h" +#ifndef TINY + +asm(".ident\t\"\\n\\n\ +Double-precision math functions (MIT License)\\n\ +Copyright 2018 ARM Limited\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* 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. + * ==================================================== + */ +/* double log1p(double x) + * Return the natural logarithm of 1+x. + * + * Method : + * 1. Argument Reduction: find k and f such that + * 1+x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * Note. If k=0, then f=x is exact. However, if k!=0, then f + * may not be representable exactly. In that case, a correction + * term is need. Let u=1+x rounded. Let c = (1+x)-u, then + * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u), + * and add back the correction term c/u. + * (Note: when x > 2**53, one can simply return log(x)) + * + * 2. Approximation of log(1+f): See log.c + * + * 3. Finally, log1p(x) = k*ln2 + log(1+f) + c/u. See log.c + * + * Special cases: + * log1p(x) is NaN with signal if x < -1 (including -INF) ; + * log1p(+INF) is +INF; log1p(-1) is -INF with signal; + * log1p(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + * + * Note: Assuming log() return accurate answer, the following + * algorithm can be used to compute log1p(x) to within a few ULP: + * + * u = 1+x; + * if(u==1.0) return x ; else + * return log(u)*(x/(u-1.0)); + * + * See HP-15C Advanced Functions Handbook, p.193. + */ + +static const double +ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ +ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ +Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ +Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ +Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ +Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ +Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ +Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ +Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + +/** + * Returns log(𝟷+𝑥). + */ +double log1p(double x) +{ + union {double f; uint64_t i;} u = {x}; + double_t hfsq,f,c,s,z,R,w,t1,t2,dk; + uint32_t hx,hu; + int k; + + hx = u.i>>32; + k = 1; + if (hx < 0x3fda827a || hx>>31) { /* 1+x < sqrt(2)+ */ + if (hx >= 0xbff00000) { /* x <= -1.0 */ + if (x == -1) + return x/0.0; /* log1p(-1) = -inf */ + return (x-x)/0.0; /* log1p(x<-1) = NaN */ + } + if (hx<<1 < 0x3ca00000<<1) { /* |x| < 2**-53 */ + /* underflow if subnormal */ + if ((hx&0x7ff00000) == 0) + FORCE_EVAL((float)x); + return x; + } + if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ + k = 0; + c = 0; + f = x; + } + } else if (hx >= 0x7ff00000) + return x; + if (k) { + u.f = 1 + x; + hu = u.i>>32; + 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.f-x) : x-(u.f-1); + c /= u.f; + } else + c = 0; + /* reduce u into [sqrt(2)/2, sqrt(2)] */ + hu = (hu&0x000fffff) + 0x3fe6a09e; + u.i = (uint64_t)hu<<32 | (u.i&0xffffffff); + f = u.f - 1; + } + hfsq = 0.5*f*f; + s = f/(2.0+f); + z = s*s; + w = z*z; + t1 = w*(Lg2+w*(Lg4+w*Lg6)); + t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + R = t2 + t1; + dk = k; + return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi; +} + +#endif /* TINY */ diff --git a/libc/tinymath/pow.c b/libc/tinymath/pow.c index 4b6a83d90..898fcbc22 100644 --- a/libc/tinymath/pow.c +++ b/libc/tinymath/pow.c @@ -288,6 +288,7 @@ static inline int zeroinfnan(uint64_t i) /** * Returns 𝑥^𝑦. + * @note should take ~18ns */ double pow(double x, double y) { diff --git a/libc/tinymath/powf.S b/libc/tinymath/powf-tiny.S similarity index 98% rename from libc/tinymath/powf.S rename to libc/tinymath/powf-tiny.S index 7ff1ac8b8..c65670167 100644 --- a/libc/tinymath/powf.S +++ b/libc/tinymath/powf-tiny.S @@ -17,6 +17,7 @@ │ PERFORMANCE OF THIS SOFTWARE. │ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/macros.internal.h" +#ifdef TINY // Returns 𝑥^𝑦. // @@ -26,3 +27,5 @@ powf: ezlea powl,ax jmp _f2ld2 .endfn powf,globl + +#endif /* TINY */ diff --git a/libc/tinymath/powf.c b/libc/tinymath/powf.c new file mode 100644 index 000000000..0d8dd07a0 --- /dev/null +++ b/libc/tinymath/powf.c @@ -0,0 +1,226 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/bits/likely.h" +#include "libc/math.h" +#include "libc/tinymath/exp2f_data.internal.h" +#include "libc/tinymath/exp_data.internal.h" +#include "libc/tinymath/internal.h" +#include "libc/tinymath/powf_data.internal.h" +#ifndef TINY + +asm(".ident\t\"\\n\\n\ +Double-precision math functions (MIT License)\\n\ +Copyright 2018 ARM Limited\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* + * Copyright (c) 2017-2018, Arm Limited. + * SPDX-License-Identifier: MIT + */ + +/* +POWF_LOG2_POLY_ORDER = 5 +EXP2F_TABLE_BITS = 5 + +ULP error: 0.82 (~ 0.5 + relerr*2^24) +relerr: 1.27 * 2^-26 (Relative error ~= 128*Ln2*relerr_log2 + relerr_exp2) +relerr_log2: 1.83 * 2^-33 (Relative error of logx.) +relerr_exp2: 1.69 * 2^-34 (Relative error of exp2(ylogx).) +*/ + +#define N (1 << POWF_LOG2_TABLE_BITS) +#define T __powf_log2_data.tab +#define A __powf_log2_data.poly +#define OFF 0x3f330000 + +/* Subnormal input is normalized so ix has negative biased exponent. + Output is multiplied by N (POWF_SCALE) if TOINT_INTRINICS is set. */ +static inline double_t log2_inline(uint32_t ix) +{ + double_t z, r, r2, r4, p, q, y, y0, invc, logc; + uint32_t iz, top, tmp; + int k, i; + + /* x = 2^k z; where z is in range [OFF,2*OFF] and exact. + The range is split into N subintervals. + The ith subinterval contains z and c is near its center. */ + tmp = ix - OFF; + i = (tmp >> (23 - POWF_LOG2_TABLE_BITS)) % N; + top = tmp & 0xff800000; + iz = ix - top; + k = (int32_t)top >> (23 - POWF_SCALE_BITS); /* arithmetic shift */ + invc = T[i].invc; + logc = T[i].logc; + z = (double_t)asfloat(iz); + + /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ + r = z * invc - 1; + y0 = logc + (double_t)k; + + /* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ + r2 = r * r; + y = A[0] * r + A[1]; + p = A[2] * r + A[3]; + r4 = r2 * r2; + q = A[4] * r + y0; + q = p * r2 + q; + y = y * r4 + q; + return y; +} + +#undef N +#undef T +#define N (1 << EXP2F_TABLE_BITS) +#define T __exp2f_data.tab +#define SIGN_BIAS (1 << (EXP2F_TABLE_BITS + 11)) + +/* The output of log2 and thus the input of exp2 is either scaled by N + (in case of fast toint intrinsics) or not. The unscaled xd must be + in [-1021,1023], sign_bias sets the sign of the result. */ +static inline float exp2_inline(double_t xd, uint32_t sign_bias) +{ + uint64_t ki, ski, t; + double_t kd, z, r, r2, y, s; + +#if TOINT_INTRINSICS +#define C __exp2f_data.poly_scaled + /* N*x = k + r with r in [-1/2, 1/2] */ + kd = roundtoint(xd); /* k */ + ki = converttoint(xd); +#else +#define C __exp2f_data.poly +#define SHIFT __exp2f_data.shift_scaled + /* x = k/N + r with r in [-1/(2N), 1/(2N)] */ + kd = eval_as_double(xd + SHIFT); + ki = asuint64(kd); + kd -= SHIFT; /* k/N */ +#endif + r = xd - kd; + + /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + ski = ki + sign_bias; + t += ski << (52 - EXP2F_TABLE_BITS); + s = asdouble(t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return eval_as_float(y); +} + +/* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is + the bit representation of a non-zero finite floating-point value. */ +static inline int checkint(uint32_t iy) +{ + int e = iy >> 23 & 0xff; + if (e < 0x7f) + return 0; + if (e > 0x7f + 23) + return 2; + if (iy & ((1 << (0x7f + 23 - e)) - 1)) + return 0; + if (iy & (1 << (0x7f + 23 - e))) + return 1; + return 2; +} + +static inline int zeroinfnan(uint32_t ix) +{ + return 2 * ix - 1 >= 2u * 0x7f800000 - 1; +} + +/** + * Returns 𝑥^𝑦. + * @note should take ~16ns + */ +float powf(float x, float y) +{ + uint32_t sign_bias = 0; + uint32_t ix, iy; + + ix = asuint(x); + iy = asuint(y); + if (UNLIKELY(ix - 0x00800000 >= 0x7f800000 - 0x00800000 || + zeroinfnan(iy))) { + /* Either (x < 0x1p-126 or inf or nan) or (y is 0 or inf or nan). */ + if (UNLIKELY(zeroinfnan(iy))) { + if (2 * iy == 0) + return issignalingf_inline(x) ? x + y : 1.0f; + if (ix == 0x3f800000) + return issignalingf_inline(y) ? x + y : 1.0f; + if (2 * ix > 2u * 0x7f800000 || + 2 * iy > 2u * 0x7f800000) + return x + y; + if (2 * ix == 2 * 0x3f800000) + return 1.0f; + if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000)) + return 0.0f; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ + return y * y; + } + if (UNLIKELY(zeroinfnan(ix))) { + float_t x2 = x * x; + if (ix & 0x80000000 && checkint(iy) == 1) + x2 = -x2; + /* Without the barrier some versions of clang hoist the 1/x2 and + thus division by zero exception can be signaled spuriously. */ + return iy & 0x80000000 ? fp_barrierf(1 / x2) : x2; + } + /* x and y are non-zero finite. */ + if (ix & 0x80000000) { + /* Finite x < 0. */ + int yint = checkint(iy); + if (yint == 0) + return __math_invalidf(x); + if (yint == 1) + sign_bias = SIGN_BIAS; + ix &= 0x7fffffff; + } + if (ix < 0x00800000) { + /* Normalize subnormal x so exponent becomes negative. */ + ix = asuint(x * 0x1p23f); + ix &= 0x7fffffff; + ix -= 23 << 23; + } + } + double_t logx = log2_inline(ix); + double_t ylogx = y * logx; /* cannot overflow, y is single prec. */ + if (UNLIKELY((asuint64(ylogx) >> 47 & 0xffff) >= + asuint64(126.0 * POWF_SCALE) >> 47)) { + /* |y*log(x)| >= 126. */ + if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE) + return __math_oflowf(sign_bias); + if (ylogx <= -150.0 * POWF_SCALE) + return __math_uflowf(sign_bias); + } + return exp2_inline(ylogx, sign_bias); +} + +#endif /* TINY */ diff --git a/libc/tinymath/powf_data.c b/libc/tinymath/powf_data.c new file mode 100644 index 000000000..c6376dbdf --- /dev/null +++ b/libc/tinymath/powf_data.c @@ -0,0 +1,67 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/tinymath/powf_data.internal.h" + +asm(".ident\t\"\\n\\n\ +Double-precision math functions (MIT License)\\n\ +Copyright 2018 ARM Limited\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* + * Data definition for powf. + * + * Copyright (c) 2017-2018, Arm Limited. + * SPDX-License-Identifier: MIT + */ + +const struct powf_log2_data __powf_log2_data = { + .tab = { + { 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 * POWF_SCALE }, + { 0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2 * POWF_SCALE }, + { 0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2 * POWF_SCALE }, + { 0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2 * POWF_SCALE }, + { 0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2 * POWF_SCALE }, + { 0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3 * POWF_SCALE }, + { 0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3 * POWF_SCALE }, + { 0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4 * POWF_SCALE }, + { 0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5 * POWF_SCALE }, + { 0x1p+0, 0x0p+0 * POWF_SCALE }, + { 0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4 * POWF_SCALE }, + { 0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3 * POWF_SCALE }, + { 0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3 * POWF_SCALE }, + { 0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2 * POWF_SCALE }, + { 0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2 * POWF_SCALE }, + { 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 * POWF_SCALE }, + }, + .poly = { + 0x1.27616c9496e0bp-2 * POWF_SCALE, -0x1.71969a075c67ap-2 * POWF_SCALE, + 0x1.ec70a6ca7baddp-2 * POWF_SCALE, -0x1.7154748bef6c8p-1 * POWF_SCALE, + 0x1.71547652ab82bp0 * POWF_SCALE, + } +}; diff --git a/libc/tinymath/powf_data.internal.h b/libc/tinymath/powf_data.internal.h new file mode 100644 index 000000000..1b0bf87f9 --- /dev/null +++ b/libc/tinymath/powf_data.internal.h @@ -0,0 +1,25 @@ +#ifndef COSMOPOLITAN_LIBC_TINYMATH_POWF_DATA_INTERNAL_H_ +#define COSMOPOLITAN_LIBC_TINYMATH_POWF_DATA_INTERNAL_H_ + +#define POWF_LOG2_TABLE_BITS 4 +#define POWF_LOG2_POLY_ORDER 5 +#if TOINT_INTRINSICS +#define POWF_SCALE_BITS EXP2F_TABLE_BITS +#else +#define POWF_SCALE_BITS 0 +#endif +#define POWF_SCALE ((double)(1 << POWF_SCALE_BITS)) + +#if !(__ASSEMBLER__ + __LINKER__ + 0) +COSMOPOLITAN_C_START_ + +extern hidden const struct powf_log2_data { + struct { + double invc, logc; + } tab[1 << POWF_LOG2_TABLE_BITS]; + double poly[POWF_LOG2_POLY_ORDER]; +} __powf_log2_data; + +COSMOPOLITAN_C_END_ +#endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */ +#endif /* COSMOPOLITAN_LIBC_TINYMATH_POWF_DATA_INTERNAL_H_ */ diff --git a/libc/tinymath/powl.c b/libc/tinymath/powl.c index e2a58a23b..40a73e455 100644 --- a/libc/tinymath/powl.c +++ b/libc/tinymath/powl.c @@ -21,6 +21,7 @@ /** * Returns 𝑥^𝑦. + * @note should take ~56ns */ long double powl(long double x, long double y) { long double t, u; diff --git a/libc/tinymath/rempio2.c b/libc/tinymath/rempio2.c index 6f257a92a..a14fda98a 100644 --- a/libc/tinymath/rempio2.c +++ b/libc/tinymath/rempio2.c @@ -36,8 +36,8 @@ asm(".ident\t\"\\n\\n\ Musl libc (MIT License)\\n\ Copyright 2005-2014 Rich Felker, et. al.\""); asm(".include \"libc/disclaimer.inc\""); - /* clang-format off */ + /* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c */ /* * ==================================================== diff --git a/libc/tinymath/rempio2f.c b/libc/tinymath/rempio2f.c new file mode 100644 index 000000000..c99d77d81 --- /dev/null +++ b/libc/tinymath/rempio2f.c @@ -0,0 +1,124 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/bits/likely.h" +#include "libc/math.h" +#include "libc/tinymath/kernel.internal.h" + +asm(".ident\t\"\\n\\n\ +fdlibm (fdlibm license)\\n\ +Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\""); +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2f.c */ +/* + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Debugged and optimized by Bruce D. Evans. + */ +/* + * ==================================================== + * 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. + * ==================================================== + */ +/* __rem_pio2f(x,y) + * + * return the remainder of x rem pi/2 in *y + * use double precision for everything except passing x + * use __rem_pio2_large() for large x + */ + +#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 +#define EPS DBL_EPSILON +#elif FLT_EVAL_METHOD==2 +#define EPS LDBL_EPSILON +#endif + +/* + * invpio2: 53 bits of 2/pi + * pio2_1: first 25 bits of pi/2 + * pio2_1t: pi/2 - pio2_1 + */ +static const double +toint = 1.5/EPS, +pio4 = 0x1.921fb6p-1, +invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ +pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ +pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ + +int __rem_pio2f(float x, double *y) +{ + union {float f; uint32_t i;} u = {x}; + double tx[1],ty[1]; + double_t fn; + uint32_t ix; + int n, sign, e0; + + ix = u.i & 0x7fffffff; + /* 25+53 bit pi is good enough for medium size */ + if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ + /* Use a specialized rint() to get fn. */ + fn = (double_t)x*invpio2 + toint - toint; + n = (int32_t)fn; + *y = x - fn*pio2_1 - fn*pio2_1t; + /* Matters with directed rounding. */ + if (UNLIKELY(*y < -pio4)) { + n--; + fn--; + *y = x - fn*pio2_1 - fn*pio2_1t; + } else if (UNLIKELY(*y > pio4)) { + n++; + fn++; + *y = x - fn*pio2_1 - fn*pio2_1t; + } + return n; + } + if(ix>=0x7f800000) { /* x is inf or NaN */ + *y = x-x; + return 0; + } + /* scale x into [2^23, 2^24-1] */ + sign = u.i>>31; + e0 = (ix>>23) - (0x7f+23); /* e0 = ilogb(|x|)-23, positive */ + u.i = ix - (e0<<23); + tx[0] = u.f; + n = __rem_pio2_large(tx,ty,e0,1,0); + if (sign) { + *y = -ty[0]; + return -n; + } + *y = ty[0]; + return n; +} diff --git a/libc/tinymath/scalbln.c b/libc/tinymath/scalbln.c new file mode 100644 index 000000000..a783482e2 --- /dev/null +++ b/libc/tinymath/scalbln.c @@ -0,0 +1,26 @@ +/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│ +│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│ +╞══════════════════════════════════════════════════════════════════════════════╡ +│ Copyright 2022 Justine Alexandra Roberts Tunney │ +│ │ +│ Permission to use, copy, modify, and/or distribute this software for │ +│ any purpose with or without fee is hereby granted, provided that the │ +│ above copyright notice and this permission notice appear in all copies. │ +│ │ +│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ +│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ +│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ +│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ +│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ +│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ +│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ +│ PERFORMANCE OF THIS SOFTWARE. │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" + +/** + * Returns 𝑥 × 2ʸ. + */ +double scalbln(double x, long n) { + return ldexp(x, n); +} diff --git a/libc/tinymath/scalblnf.c b/libc/tinymath/scalblnf.c new file mode 100644 index 000000000..dfa88ffaa --- /dev/null +++ b/libc/tinymath/scalblnf.c @@ -0,0 +1,26 @@ +/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│ +│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│ +╞══════════════════════════════════════════════════════════════════════════════╡ +│ Copyright 2022 Justine Alexandra Roberts Tunney │ +│ │ +│ Permission to use, copy, modify, and/or distribute this software for │ +│ any purpose with or without fee is hereby granted, provided that the │ +│ above copyright notice and this permission notice appear in all copies. │ +│ │ +│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ +│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ +│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ +│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ +│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ +│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ +│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ +│ PERFORMANCE OF THIS SOFTWARE. │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" + +/** + * Returns 𝑥 × 2ʸ. + */ +float scalblnf(float x, long n) { + return ldexpf(x, n); +} diff --git a/libc/tinymath/scalbn.c b/libc/tinymath/scalbn.c new file mode 100644 index 000000000..0e28095d5 --- /dev/null +++ b/libc/tinymath/scalbn.c @@ -0,0 +1,26 @@ +/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│ +│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│ +╞══════════════════════════════════════════════════════════════════════════════╡ +│ Copyright 2022 Justine Alexandra Roberts Tunney │ +│ │ +│ Permission to use, copy, modify, and/or distribute this software for │ +│ any purpose with or without fee is hereby granted, provided that the │ +│ above copyright notice and this permission notice appear in all copies. │ +│ │ +│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ +│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ +│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ +│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ +│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ +│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ +│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ +│ PERFORMANCE OF THIS SOFTWARE. │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" + +/** + * Returns 𝑥 × 2ʸ. + */ +double scalbn(double x, int n) { + return ldexp(x, n); +} diff --git a/libc/tinymath/scalbnf.c b/libc/tinymath/scalbnf.c new file mode 100644 index 000000000..86e627507 --- /dev/null +++ b/libc/tinymath/scalbnf.c @@ -0,0 +1,26 @@ +/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│ +│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│ +╞══════════════════════════════════════════════════════════════════════════════╡ +│ Copyright 2022 Justine Alexandra Roberts Tunney │ +│ │ +│ Permission to use, copy, modify, and/or distribute this software for │ +│ any purpose with or without fee is hereby granted, provided that the │ +│ above copyright notice and this permission notice appear in all copies. │ +│ │ +│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ +│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ +│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ +│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ +│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ +│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ +│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ +│ PERFORMANCE OF THIS SOFTWARE. │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" + +/** + * Returns 𝑥 × 2ʸ. + */ +float scalbnf(float x, int n) { + return ldexpf(x, n); +} diff --git a/libc/tinymath/sin.c b/libc/tinymath/sin.c index 26a7c5643..e29ebd555 100644 --- a/libc/tinymath/sin.c +++ b/libc/tinymath/sin.c @@ -36,8 +36,8 @@ asm(".ident\t\"\\n\\n\ Musl libc (MIT License)\\n\ Copyright 2005-2014 Rich Felker, et. al.\""); asm(".include \"libc/disclaimer.inc\""); - /* clang-format off */ + /* origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */ /* * ==================================================== @@ -83,6 +83,10 @@ asm(".include \"libc/disclaimer.inc\""); #define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i #define gethighw(hi,d) (hi) = asuint64(d) >> 32 +/** + * Returns sine of 𝑥. + * @note should take ~5ns + */ double sin(double x) { double y[2]; diff --git a/libc/tinymath/sincos.c b/libc/tinymath/sincos.c index 06cb4d91f..32c0264cb 100644 --- a/libc/tinymath/sincos.c +++ b/libc/tinymath/sincos.c @@ -54,6 +54,10 @@ asm(".include \"libc/disclaimer.inc\""); #define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i #define gethighw(hi,d) (hi) = asuint64(d) >> 32 +/** + * Returns sine and cosine of 𝑥. + * @note should take ~10ns + */ void sincos(double x, double *sin, double *cos) { double y[2], s, c; diff --git a/libc/tinymath/sinf.S b/libc/tinymath/sinf.S deleted file mode 100644 index 4c4214f59..000000000 --- a/libc/tinymath/sinf.S +++ /dev/null @@ -1,28 +0,0 @@ -/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│ -│vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi│ -╞══════════════════════════════════════════════════════════════════════════════╡ -│ Copyright 2020 Justine Alexandra Roberts Tunney │ -│ │ -│ Permission to use, copy, modify, and/or distribute this software for │ -│ any purpose with or without fee is hereby granted, provided that the │ -│ above copyright notice and this permission notice appear in all copies. │ -│ │ -│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │ -│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │ -│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │ -│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │ -│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │ -│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │ -│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │ -│ PERFORMANCE OF THIS SOFTWARE. │ -╚─────────────────────────────────────────────────────────────────────────────*/ -#include "libc/macros.internal.h" - -// Returns sine of 𝑥. -// -// @param 𝑥 is float scalar in low quarter of %xmm0 -// @return float scalar in low quarter of %xmm0 -// @domain -(3π/8) < 𝑥 < 3π/8 for best accuracy -sinf: ezlea sinl,ax - jmp _f2ld2 - .endfn sinf,globl diff --git a/libc/tinymath/sinf.c b/libc/tinymath/sinf.c new file mode 100644 index 000000000..2eff94dfc --- /dev/null +++ b/libc/tinymath/sinf.c @@ -0,0 +1,119 @@ +/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ +│vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ +╚──────────────────────────────────────────────────────────────────────────────╝ +│ │ +│ Musl Libc │ +│ Copyright © 2005-2014 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#include "libc/math.h" +#include "libc/tinymath/complex.internal.h" +#include "libc/tinymath/feval.internal.h" +#include "libc/tinymath/kernel.internal.h" + +asm(".ident\t\"\\n\\n\ +fdlibm (fdlibm license)\\n\ +Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\""); +asm(".ident\t\"\\n\\n\ +Musl libc (MIT License)\\n\ +Copyright 2005-2014 Rich Felker, et. al.\""); +asm(".include \"libc/disclaimer.inc\""); +/* clang-format off */ + +/* origin: FreeBSD /usr/src/lib/msun/src/s_sinf.c */ +/* + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Optimized by Bruce D. Evans. + */ +/* + * ==================================================== + * 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. + * ==================================================== + */ + +/* Small multiples of pi/2 rounded to double precision. */ +static const double +s1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */ +s2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */ +s3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */ +s4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */ + +/** + * Returns sine of 𝑥. + * @note should take about 5ns + */ +float sinf(float x) +{ + double y; + uint32_t ix; + int n, sign; + + GET_FLOAT_WORD(ix, x); + sign = ix >> 31; + ix &= 0x7fffffff; + + if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */ + if (ix < 0x39800000) { /* |x| < 2**-12 */ + /* raise inexact if x!=0 and underflow if subnormal */ + FORCE_EVAL(ix < 0x00800000 ? x/0x1p120f : x+0x1p120f); + return x; + } + return __sindf(x); + } + if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */ + if (ix <= 0x4016cbe3) { /* |x| ~<= 3pi/4 */ + if (sign) + return -__cosdf(x + s1pio2); + else + return __cosdf(x - s1pio2); + } + return __sindf(sign ? -(x + s2pio2) : -(x - s2pio2)); + } + if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */ + if (ix <= 0x40afeddf) { /* |x| ~<= 7*pi/4 */ + if (sign) + return __cosdf(x + s3pio2); + else + return -__cosdf(x - s3pio2); + } + return __sindf(sign ? x + s4pio2 : x - s4pio2); + } + + /* sin(Inf or NaN) is NaN */ + if (ix >= 0x7f800000) + return x - x; + + /* general argument reduction needed */ + n = __rem_pio2f(x, &y); + switch (n&3) { + case 0: return __sindf(y); + case 1: return __cosdf(y); + case 2: return __sindf(-y); + default: + return -__cosdf(y); + } +} diff --git a/libc/tinymath/tan.c b/libc/tinymath/tan.c index 350424d54..dd1163c49 100644 --- a/libc/tinymath/tan.c +++ b/libc/tinymath/tan.c @@ -83,6 +83,9 @@ asm(".include \"libc/disclaimer.inc\""); * TRIG(x) returns trig(x) nearly rounded */ +/** + * Returns tangent of x. + */ double tan(double x) { double y[2]; diff --git a/test/libc/tinymath/acos_test.c b/test/libc/tinymath/acos_test.c index d22cb7deb..0a8b22a71 100644 --- a/test/libc/tinymath/acos_test.c +++ b/test/libc/tinymath/acos_test.c @@ -44,10 +44,11 @@ TEST(acos, test) { EXPECT_TRUE(isnan(acos(__DBL_MAX__))); } -BENCH(acos, bench) { - EZBENCH2("acos(+0)", donothing, acos(0)); - EZBENCH2("acos(-0)", donothing, acos(-0.)); - EZBENCH2("acos(NAN)", donothing, acos(NAN)); - EZBENCH2("acos(INFINITY)", donothing, acos(INFINITY)); - EZBENCH_C("acos", _real1(vigna()), acos(_real1(vigna()))); +BENCH(acosl, bench) { + double _acos(double) asm("acos"); + float _acosf(float) asm("acosf"); + long double _acosl(long double) asm("acosl"); + EZBENCH2("-acos", donothing, _acos(.7)); /* ~17ns */ + EZBENCH2("-acosf", donothing, _acosf(.7)); /* ~11ns */ + EZBENCH2("-acosl", donothing, _acosl(.7)); /* ~40ns */ } diff --git a/test/libc/tinymath/asin_test.c b/test/libc/tinymath/asin_test.c index ed35e42ef..778d8afef 100644 --- a/test/libc/tinymath/asin_test.c +++ b/test/libc/tinymath/asin_test.c @@ -44,10 +44,11 @@ TEST(asin, test) { EXPECT_TRUE(isnan(asin(__DBL_MAX__))); } -BENCH(asin, bench) { - EZBENCH2("asin(+0)", donothing, asin(0)); - EZBENCH2("asin(-0)", donothing, asin(-0.)); - EZBENCH2("asin(NAN)", donothing, asin(NAN)); - EZBENCH2("asin(INFINITY)", donothing, asin(INFINITY)); - EZBENCH_C("asin", _real1(lemur64()), asin(_real1(lemur64()))); +BENCH(asinl, bench) { + double _asin(double) asm("asin"); + float _asinf(float) asm("asinf"); + long double _asinl(long double) asm("asinl"); + EZBENCH2("-asin", donothing, _asin(.7)); /* ~16ns */ + EZBENCH2("-asinf", donothing, _asinf(.7)); /* ~12ns */ + EZBENCH2("-asinl", donothing, _asinl(.7)); /* ~39ns */ } diff --git a/test/libc/tinymath/atan_test.c b/test/libc/tinymath/atan_test.c index 34686b484..dbd435c4d 100644 --- a/test/libc/tinymath/atan_test.c +++ b/test/libc/tinymath/atan_test.c @@ -18,6 +18,7 @@ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/math.h" #include "libc/runtime/gc.internal.h" +#include "libc/testlib/ezbench.h" #include "libc/testlib/testlib.h" #include "libc/x/x.h" @@ -38,3 +39,12 @@ TEST(atan, test) { gc(xasprintf("%.15g", atan(__DBL_MIN__)))); EXPECT_STREQ("1.5707963267949", gc(xasprintf("%.15g", atan(__DBL_MAX__)))); } + +BENCH(atanl, bench) { + double _atan(double) asm("atan"); + float _atanf(float) asm("atanf"); + long double _atanl(long double) asm("atanl"); + EZBENCH2("-atan", donothing, _atan(.7)); /* ~18ns */ + EZBENCH2("-atanf", donothing, _atanf(.7)); /* ~12ns */ + EZBENCH2("-atanl", donothing, _atanl(.7)); /* ~34ns */ +} diff --git a/test/libc/tinymath/cos_test.c b/test/libc/tinymath/cos_test.c index cb2566c53..2fc1f13f9 100644 --- a/test/libc/tinymath/cos_test.c +++ b/test/libc/tinymath/cos_test.c @@ -18,6 +18,7 @@ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/math.h" #include "libc/runtime/gc.internal.h" +#include "libc/testlib/ezbench.h" #include "libc/testlib/testlib.h" #include "libc/x/x.h" @@ -43,3 +44,12 @@ TEST(cos, test) { gc(xasprintf("%.15g", cos(-1.0000000000000002)))); EXPECT_STREQ("1", gc(xasprintf("%.15g", cos(-2.1073424255447e-08)))); } + +BENCH(cos, bench) { + double _cos(double) asm("cos"); + float _cosf(float) asm("cosf"); + long double _cosl(long double) asm("cosl"); + EZBENCH2("cos", donothing, _cos(.7)); /* ~6ns */ + EZBENCH2("cosf", donothing, _cosf(.7)); /* ~5ns */ + EZBENCH2("cosl", donothing, _cosl(.7)); /* ~28ns */ +} diff --git a/test/libc/tinymath/ldexp_test.c b/test/libc/tinymath/ldexp_test.c index 382b8a3af..eb291d559 100644 --- a/test/libc/tinymath/ldexp_test.c +++ b/test/libc/tinymath/ldexp_test.c @@ -20,6 +20,7 @@ #include "libc/rand/rand.h" #include "libc/runtime/gc.internal.h" #include "libc/stdio/stdio.h" +#include "libc/testlib/ezbench.h" #include "libc/testlib/testlib.h" #include "libc/x/x.h" @@ -110,7 +111,11 @@ TEST(ldexp, stuff) { ASSERT_STREQ("100.48", gc(xasprintf("%.2Lf", scalblnl(pi, twopow)))); } -TEST(exp10, test) { - ASSERT_EQ(100, (int)exp10(2)); - ASSERT_STREQ("100.000000", gc(xasprintf("%Lf", exp10l(2)))); +BENCH(ldexpl, bench) { + double _ldexp(double, int) asm("ldexp"); + float _ldexpf(float, int) asm("ldexpf"); + long double _ldexpl(long double, int) asm("ldexpl"); + EZBENCH2("ldexp", donothing, _ldexp(.7, 3)); /* ~1ns */ + EZBENCH2("ldexpf", donothing, _ldexpf(.7, 3)); /* ~1ns */ + EZBENCH2("ldexpl", donothing, _ldexpl(.7, 3)); /* ~7ns */ } diff --git a/test/libc/tinymath/log1p_test.c b/test/libc/tinymath/log1p_test.c index 2833e112b..0d1cbb429 100644 --- a/test/libc/tinymath/log1p_test.c +++ b/test/libc/tinymath/log1p_test.c @@ -18,6 +18,7 @@ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/math.h" #include "libc/runtime/gc.internal.h" +#include "libc/testlib/ezbench.h" #include "libc/testlib/testlib.h" #include "libc/x/x.h" @@ -45,3 +46,12 @@ TEST(log1pf, test) { EXPECT_STREQ("-INFINITY", gc(xdtoaf(log1pf(-1)))); EXPECT_STREQ("-NAN", gc(xdtoaf(log1pf(-2)))); } + +BENCH(log1p, bench) { + double _log1p(double) asm("log1p"); + float _log1pf(float) asm("log1pf"); + long double _log1pl(long double) asm("log1pl"); + EZBENCH2("log1p", donothing, _log1p(.7)); /* ~17ns */ + EZBENCH2("log1pf", donothing, _log1pf(.7)); /* ~17ns */ + EZBENCH2("log1pl", donothing, _log1pl(.7)); /* ~28ns */ +} diff --git a/test/libc/tinymath/powl_test.c b/test/libc/tinymath/powl_test.c index db90d1ffd..64a1ab0bb 100644 --- a/test/libc/tinymath/powl_test.c +++ b/test/libc/tinymath/powl_test.c @@ -563,6 +563,6 @@ BENCH(powl, bench) { float _powf(float, float) asm("powf"); long double _powl(long double, long double) asm("powl"); EZBENCH2("pow", donothing, _pow(.7, .2)); /* ~18ns */ - EZBENCH2("powf", donothing, _powf(.7, .2)); /* ~56ns */ + EZBENCH2("powf", donothing, _powf(.7, .2)); /* ~16ns */ EZBENCH2("powl", donothing, _powl(.7, .2)); /* ~56ns */ } diff --git a/test/libc/tinymath/sin_test.c b/test/libc/tinymath/sin_test.c index 5494bb36c..22a39c32e 100644 --- a/test/libc/tinymath/sin_test.c +++ b/test/libc/tinymath/sin_test.c @@ -72,8 +72,11 @@ TEST(sinf, test) { EXPECT_STARTSWITH(".873283", gc(xdtoaf(sinf(555)))); } -BENCH(sinl, bench) { - EZBENCH(donothing, sinl(0.7)); /* ~30ns */ - EZBENCH(donothing, sin(0.7)); /* ~35ns */ - EZBENCH(donothing, sinf(0.7f)); /* ~35ns */ +BENCH(sin, bench) { + double _sin(double) asm("sin"); + float _sinf(float) asm("sinf"); + long double _sinl(long double) asm("sinl"); + EZBENCH2("sin", donothing, _sin(.7)); /* ~5ns */ + EZBENCH2("sinf", donothing, _sinf(.7)); /* ~5ns */ + EZBENCH2("sinl", donothing, _sinl(.7)); /* ~28ns */ } diff --git a/test/libc/tinymath/sincos_test.c b/test/libc/tinymath/sincos_test.c index 746a4daf1..7fdd0ee05 100644 --- a/test/libc/tinymath/sincos_test.c +++ b/test/libc/tinymath/sincos_test.c @@ -29,9 +29,39 @@ TEST(sincos, test) { EXPECT_STREQ("0.995004165278026", gc(xasprintf("%.15g", cosine))); } -BENCH(sincos, bench) { - volatile double x = 31337; - volatile double sine, cosine; - EZBENCH2("sin+cos", donothing, (sin(x), cos(x))); - EZBENCH2("sincos", donothing, sincos(x, &sine, &cosine)); +TEST(sincosf, test) { + float sine, cosine; + sincosf(.1, &sine, &cosine); + EXPECT_STREQ("0.0998334", gc(xasprintf("%.6g", sine))); + EXPECT_STREQ("0.995004", gc(xasprintf("%.6g", cosine))); +} + +TEST(sincosl, test) { + long double sine, cosine; + sincosl(.1, &sine, &cosine); + EXPECT_STREQ("0.0998334166468282", gc(xasprintf("%.15Lg", sine))); + EXPECT_STREQ("0.995004165278026", gc(xasprintf("%.15Lg", cosine))); +} + +#define NUM .123 + +BENCH(sincos, bench) { + double _sin(double) asm("sin"); + float _sinf(float) asm("sinf"); + long double _sinl(long double) asm("sinl"); + double _cos(double) asm("cos"); + float _cosf(float) asm("cosf"); + long double _cosl(long double) asm("cosl"); + double _sincos(double, double*, double*) asm("sincos"); + float _sincosf(float, float*, float*) asm("sincosf"); + long double _sincosl(long double, long double*, long double*) asm("sincosl"); + volatile float sinef, cosinef; + volatile double sine, cosine; + volatile long double sinel, cosinel; + EZBENCH2("sin+cos", donothing, (_sin(NUM), _cos(NUM))); + EZBENCH2("sincos", donothing, _sincos(NUM, &sine, &cosine)); + EZBENCH2("sinf+cosf", donothing, (_sinf(NUM), _cosf(NUM))); + EZBENCH2("sincosf", donothing, _sincosf(NUM, &sinef, &cosinef)); + EZBENCH2("sinl+cosl", donothing, (_sinl(NUM), _cosl(NUM))); + EZBENCH2("sincosl", donothing, _sincosl(NUM, &sinel, &cosinel)); } diff --git a/test/libc/tinymath/tan_test.c b/test/libc/tinymath/tan_test.c index 0131d669e..bfce11b0c 100644 --- a/test/libc/tinymath/tan_test.c +++ b/test/libc/tinymath/tan_test.c @@ -18,6 +18,7 @@ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/math.h" #include "libc/runtime/gc.internal.h" +#include "libc/testlib/ezbench.h" #include "libc/testlib/testlib.h" #include "libc/x/x.h" @@ -38,3 +39,12 @@ TEST(tan, test) { gc(xasprintf("%.15g", tan(__DBL_MIN__)))); EXPECT_STREQ("-0.0049620158744449", gc(xasprintf("%.15g", tan(__DBL_MAX__)))); } + +BENCH(tan, bench) { + double _tan(double) asm("tan"); + float _tanf(float) asm("tanf"); + long double _tanl(long double) asm("tanl"); + EZBENCH2("tan", donothing, _tan(.7)); /* ~19ns */ + EZBENCH2("tanf", donothing, _tanf(.7)); /* ~32ns */ + EZBENCH2("tanl", donothing, _tanl(.7)); /* ~28ns */ +} diff --git a/test/tool/net/ljson_test.lua b/test/tool/net/ljson_test.lua index f270c95dc..6bfe364c2 100644 --- a/test/tool/net/ljson_test.lua +++ b/test/tool/net/ljson_test.lua @@ -26,6 +26,7 @@ assert(EncodeLua(assert(DecodeJson[[ 9.123e6 ]])) == '9123000.') assert(EncodeLua(assert(DecodeJson[[ [{"heh": [1,3,2]}] ]])) == '{{heh={1, 3, 2}}}') assert(EncodeLua(assert(DecodeJson[[ 3.14159 ]])) == '3.14159') assert(EncodeLua(assert(DecodeJson[[ 1e-12 ]])) == '1e-12') +assert(assert(DecodeJson[[ "\u007f" ]]) == '\x7f') assert(EncodeJson(assert(DecodeJson[[ 1e-12 ]])) == '1e-12') assert(EncodeJson(assert(DecodeJson[[ true ]])) == 'true') @@ -162,11 +163,11 @@ function JsonEncodeObject() EncodeJson({["3"]="1", ["4"]="1", ["5"]={["3"]="1", ["4"]="1", ["5"]="9"}}) end -print('JsonParseEmpty', Benchmark(JsonParseEmpty)) -print('JsonParseInteg', Benchmark(JsonParseInteger)) -print('JsonParseDouble', Benchmark(JsonParseDouble)) -print('JsonParseString', Benchmark(JsonParseString)) -print('JsonParseArray', Benchmark(JsonParseArray)) -print('JsonParseObject', Benchmark(JsonParseObject)) -print('JsonEncodeArr', Benchmark(JsonEncodeArray)) -print('JsonEncodeObj', Benchmark(JsonEncodeObject)) +-- print('JsonParseEmpty', Benchmark(JsonParseEmpty)) +-- print('JsonParseInteg', Benchmark(JsonParseInteger)) +-- print('JsonParseDouble', Benchmark(JsonParseDouble)) +-- print('JsonParseString', Benchmark(JsonParseString)) +-- print('JsonParseArray', Benchmark(JsonParseArray)) +-- print('JsonParseObject', Benchmark(JsonParseObject)) +-- print('JsonEncodeArr', Benchmark(JsonEncodeArray)) +-- print('JsonEncodeObj', Benchmark(JsonEncodeObject)) diff --git a/tool/net/ljson.c b/tool/net/ljson.c index 7023b8429..ec358dbc6 100644 --- a/tool/net/ljson.c +++ b/tool/net/ljson.c @@ -306,7 +306,7 @@ static struct DecodeJson Parse(struct lua_State *L, const char *p, goto BadUnicode; } // UTF-8 - if (c < 0x7f) { + if (c <= 0x7f) { w[0] = c; i = 1; } else if (c <= 0x7ff) { @@ -346,7 +346,7 @@ static struct DecodeJson Parse(struct lua_State *L, const char *p, goto StringFailureWithReason; } } - break; + unreachable; StringFailureWithReason: luaL_pushresultsize(&b, 0); lua_pop(L, 1);