diff --git a/libc/tinymath/complex.h b/libc/tinymath/complex.h new file mode 100644 index 000000000..908098eed --- /dev/null +++ b/libc/tinymath/complex.h @@ -0,0 +1,160 @@ +/*-*- 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#ifndef _COMPLEX_H +#define _COMPLEX_H + +#ifdef __cplusplus +extern "C" { +#endif + +#define complex _Complex +#ifdef __GNUC__ +#define _Complex_I (__extension__ (0.0f+1.0fi)) +#else +#define _Complex_I (0.0f+1.0fi) +#endif +#define I _Complex_I + +double complex cacos(double complex); +float complex cacosf(float complex); +long double complex cacosl(long double complex); + +double complex casin(double complex); +float complex casinf(float complex); +long double complex casinl(long double complex); + +double complex catan(double complex); +float complex catanf(float complex); +long double complex catanl(long double complex); + +double complex ccos(double complex); +float complex ccosf(float complex); +long double complex ccosl(long double complex); + +double complex csin(double complex); +float complex csinf(float complex); +long double complex csinl(long double complex); + +double complex ctan(double complex); +float complex ctanf(float complex); +long double complex ctanl(long double complex); + +double complex cacosh(double complex); +float complex cacoshf(float complex); +long double complex cacoshl(long double complex); + +double complex casinh(double complex); +float complex casinhf(float complex); +long double complex casinhl(long double complex); + +double complex catanh(double complex); +float complex catanhf(float complex); +long double complex catanhl(long double complex); + +double complex ccosh(double complex); +float complex ccoshf(float complex); +long double complex ccoshl(long double complex); + +double complex csinh(double complex); +float complex csinhf(float complex); +long double complex csinhl(long double complex); + +double complex ctanh(double complex); +float complex ctanhf(float complex); +long double complex ctanhl(long double complex); + +double complex cexp(double complex); +float complex cexpf(float complex); +long double complex cexpl(long double complex); + +double complex clog(double complex); +float complex clogf(float complex); +long double complex clogl(long double complex); + +double cabs(double complex); +float cabsf(float complex); +long double cabsl(long double complex); + +double complex cpow(double complex, double complex); +float complex cpowf(float complex, float complex); +long double complex cpowl(long double complex, long double complex); + +double complex csqrt(double complex); +float complex csqrtf(float complex); +long double complex csqrtl(long double complex); + +double carg(double complex); +float cargf(float complex); +long double cargl(long double complex); + +double cimag(double complex); +float cimagf(float complex); +long double cimagl(long double complex); + +double complex conj(double complex); +float complex conjf(float complex); +long double complex conjl(long double complex); + +double complex cproj(double complex); +float complex cprojf(float complex); +long double complex cprojl(long double complex); + +double creal(double complex); +float crealf(float complex); +long double creall(long double complex); + +#ifndef __cplusplus +#define __CIMAG(x, t) \ + (+(union { _Complex t __z; t __xy[2]; }){(_Complex t)(x)}.__xy[1]) + +#define creal(x) ((double)(x)) +#define crealf(x) ((float)(x)) +#define creall(x) ((long double)(x)) + +#define cimag(x) __CIMAG(x, double) +#define cimagf(x) __CIMAG(x, float) +#define cimagl(x) __CIMAG(x, long double) +#endif + +#if __STDC_VERSION__ >= 201112L +#if defined(_Imaginary_I) +#define __CMPLX(x, y, t) ((t)(x) + _Imaginary_I*(t)(y)) +#elif defined(__clang__) +#define __CMPLX(x, y, t) (+(_Complex t){ (t)(x), (t)(y) }) +#else +#define __CMPLX(x, y, t) (__builtin_complex((t)(x), (t)(y))) +#endif +#define CMPLX(x, y) __CMPLX(x, y, double) +#define CMPLXF(x, y) __CMPLX(x, y, float) +#define CMPLXL(x, y) __CMPLX(x, y, long double) +#endif + +#ifdef __cplusplus +} +#endif +#endif diff --git a/libc/tinymath/complex_impl.h b/libc/tinymath/complex_impl.h new file mode 100644 index 000000000..c30ad9fd0 --- /dev/null +++ b/libc/tinymath/complex_impl.h @@ -0,0 +1,49 @@ +/*-*- 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#ifndef _COMPLEX_IMPL_H +#define _COMPLEX_IMPL_H + +#include +#include "libm.h" + +#undef __CMPLX +#undef CMPLX +#undef CMPLXF +#undef CMPLXL + +#define __CMPLX(x, y, t) \ + ((union { _Complex t __z; t __xy[2]; }){.__xy = {(x),(y)}}.__z) + +#define CMPLX(x, y) __CMPLX(x, y, double) +#define CMPLXF(x, y) __CMPLX(x, y, float) +#define CMPLXL(x, y) __CMPLX(x, y, long double) + +hidden double complex __ldexp_cexp(double complex,int); +hidden float complex __ldexp_cexpf(float complex,int); + +#endif diff --git a/libc/tinymath/csqrt.c b/libc/tinymath/csqrt.c index 61c042fb9..3e542eafd 100644 --- a/libc/tinymath/csqrt.c +++ b/libc/tinymath/csqrt.c @@ -52,7 +52,7 @@ * SUCH DAMAGE. */ -#include "complex_impl.h" +#include "libc/tinymath/complex_impl.h" /* * gcc doesn't implement complex multiplication or division correctly, diff --git a/libc/tinymath/libm.h b/libc/tinymath/libm.h new file mode 100644 index 000000000..56393d0c9 --- /dev/null +++ b/libc/tinymath/libm.h @@ -0,0 +1,301 @@ +/*-*- 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. │ +│ │ +╚─────────────────────────────────────────────────────────────────────────────*/ +#ifndef _LIBM_H +#define _LIBM_H + +#include +#include +#include +#include +// #include "fp_arch.h" + +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t m; + uint16_t se; + } i; +}; +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN +/* This is the m68k variant of 80-bit long double, and this definition only works + * on archs where the alignment requirement of uint64_t is <= 4. */ +union ldshape { + long double f; + struct { + uint16_t se; + uint16_t pad; + uint64_t m; + } i; +}; +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t lo; + uint32_t mid; + uint16_t top; + uint16_t se; + } i; + struct { + uint64_t lo; + uint64_t hi; + } i2; +}; +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN +union ldshape { + long double f; + struct { + uint16_t se; + uint16_t top; + uint32_t mid; + uint64_t lo; + } i; + struct { + uint64_t hi; + uint64_t lo; + } i2; +}; +#else +#error Unsupported long double representation +#endif + +/* Support non-nearest rounding mode. */ +#define WANT_ROUNDING 1 +/* Support signaling NaNs. */ +#define WANT_SNAN 0 + +#if WANT_SNAN +#error SNaN is unsupported +#else +#define issignalingf_inline(x) 0 +#define issignaling_inline(x) 0 +#endif + +#ifndef TOINT_INTRINSICS +#define TOINT_INTRINSICS 0 +#endif + +#if TOINT_INTRINSICS +/* Round x to nearest int in all rounding modes, ties have to be rounded + consistently with converttoint so the results match. If the result + would be outside of [-2^31, 2^31-1] then the semantics is unspecified. */ +static double_t roundtoint(double_t); + +/* Convert x to nearest int in all rounding modes, ties have to be rounded + consistently with roundtoint. If the result is not representible in an + int32_t then the semantics is unspecified. */ +static int32_t converttoint(double_t); +#endif + +/* Helps static branch prediction so hot path can be better optimized. */ +#ifdef __GNUC__ +#define predict_true(x) __builtin_expect(!!(x), 1) +#define predict_false(x) __builtin_expect(x, 0) +#else +#define predict_true(x) (x) +#define predict_false(x) (x) +#endif + +/* Evaluate an expression as the specified type. With standard excess + precision handling a type cast or assignment is enough (with + -ffloat-store an assignment is required, in old compilers argument + passing and return statement may not drop excess precision). */ + +static inline float eval_as_float(float x) +{ + float y = x; + return y; +} + +static inline double eval_as_double(double x) +{ + double y = x; + return y; +} + +/* fp_barrier returns its input, but limits code transformations + as if it had a side-effect (e.g. observable io) and returned + an arbitrary value. */ + +#ifndef fp_barrierf +#define fp_barrierf fp_barrierf +static inline float fp_barrierf(float x) +{ + volatile float y = x; + return y; +} +#endif + +#ifndef fp_barrier +#define fp_barrier fp_barrier +static inline double fp_barrier(double x) +{ + volatile double y = x; + return y; +} +#endif + +#ifndef fp_barrierl +#define fp_barrierl fp_barrierl +static inline long double fp_barrierl(long double x) +{ + volatile long double y = x; + return y; +} +#endif + +/* fp_force_eval ensures that the input value is computed when that's + otherwise unused. To prevent the constant folding of the input + expression, an additional fp_barrier may be needed or a compilation + mode that does so (e.g. -frounding-math in gcc). Then it can be + used to evaluate an expression for its fenv side-effects only. */ + +#ifndef fp_force_evalf +#define fp_force_evalf fp_force_evalf +static inline void fp_force_evalf(float x) +{ + volatile float y; + y = x; +} +#endif + +#ifndef fp_force_eval +#define fp_force_eval fp_force_eval +static inline void fp_force_eval(double x) +{ + volatile double y; + y = x; +} +#endif + +#ifndef fp_force_evall +#define fp_force_evall fp_force_evall +static inline void fp_force_evall(long double x) +{ + volatile long double y; + y = x; +} +#endif + +#define FORCE_EVAL(x) do { \ + if (sizeof(x) == sizeof(float)) { \ + fp_force_evalf(x); \ + } else if (sizeof(x) == sizeof(double)) { \ + fp_force_eval(x); \ + } else { \ + fp_force_evall(x); \ + } \ +} while(0) + +#define asuint(f) ((union{float _f; uint32_t _i;}){f})._i +#define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f +#define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i +#define asdouble(i) ((union{uint64_t _i; double _f;}){i})._f + +#define EXTRACT_WORDS(hi,lo,d) \ +do { \ + uint64_t __u = asuint64(d); \ + (hi) = __u >> 32; \ + (lo) = (uint32_t)__u; \ +} while (0) + +#define GET_HIGH_WORD(hi,d) \ +do { \ + (hi) = asuint64(d) >> 32; \ +} while (0) + +#define GET_LOW_WORD(lo,d) \ +do { \ + (lo) = (uint32_t)asuint64(d); \ +} while (0) + +#define INSERT_WORDS(d,hi,lo) \ +do { \ + (d) = asdouble(((uint64_t)(hi)<<32) | (uint32_t)(lo)); \ +} while (0) + +#define SET_HIGH_WORD(d,hi) \ + INSERT_WORDS(d, hi, (uint32_t)asuint64(d)) + +#define SET_LOW_WORD(d,lo) \ + INSERT_WORDS(d, asuint64(d)>>32, lo) + +#define GET_FLOAT_WORD(w,d) \ +do { \ + (w) = asuint(d); \ +} while (0) + +#define SET_FLOAT_WORD(d,w) \ +do { \ + (d) = asfloat(w); \ +} while (0) + +hidden int __rem_pio2_large(double*,double*,int,int,int); + +hidden int __rem_pio2(double,double*); +hidden double __sin(double,double,int); +hidden double __cos(double,double); +hidden double __tan(double,double,int); +hidden double __expo2(double,double); + +hidden int __rem_pio2f(float,double*); +hidden float __sindf(double); +hidden float __cosdf(double); +hidden float __tandf(double,int); +hidden float __expo2f(float,float); + +hidden int __rem_pio2l(long double, long double *); +hidden long double __sinl(long double, long double, int); +hidden long double __cosl(long double, long double); +hidden long double __tanl(long double, long double, int); + +hidden long double __polevll(long double, const long double *, int); +hidden long double __p1evll(long double, const long double *, int); + +extern int __signgam; +hidden double __lgamma_r(double, int *); +hidden float __lgammaf_r(float, int *); + +/* error handling functions */ +hidden float __math_xflowf(uint32_t, float); +hidden float __math_uflowf(uint32_t); +hidden float __math_oflowf(uint32_t); +hidden float __math_divzerof(uint32_t); +hidden float __math_invalidf(float); +hidden double __math_xflow(uint32_t, double); +hidden double __math_uflow(uint32_t); +hidden double __math_oflow(uint32_t); +hidden double __math_divzero(uint32_t); +hidden double __math_invalid(double); +#if LDBL_MANT_DIG != DBL_MANT_DIG +hidden long double __math_invalidl(long double); +#endif + +#endif