diff options
Diffstat (limited to 'circuitpython/lib/libm')
30 files changed, 3617 insertions, 0 deletions
diff --git a/circuitpython/lib/libm/acoshf.c b/circuitpython/lib/libm/acoshf.c new file mode 100644 index 0000000..8a8409f --- /dev/null +++ b/circuitpython/lib/libm/acoshf.c @@ -0,0 +1,32 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// acoshf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +#include "libm.h" + +#if FLT_EVAL_METHOD==2 +#undef sqrtf +#define sqrtf sqrtl +#elif FLT_EVAL_METHOD==1 +#undef sqrtf +#define sqrtf sqrt +#endif + +/* acosh(x) = log(x + sqrt(x*x-1)) */ +float acoshf(float x) +{ + union {float f; uint32_t i;} u = {x}; + uint32_t a = u.i & 0x7fffffff; + + if (a < 0x3f800000+(1<<23)) + /* |x| < 2, invalid if x < 1 or nan */ + /* up to 2ulp error in [1,1.125] */ + return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1))); + if (a < 0x3f800000+(12<<23)) + /* |x| < 0x1p12 */ + return logf(2*x - 1/(x+sqrtf(x*x-1))); + /* x >= 0x1p12 */ + return logf(x) + 0.693147180559945309417232121458176568f; +} diff --git a/circuitpython/lib/libm/asinfacosf.c b/circuitpython/lib/libm/asinfacosf.c new file mode 100644 index 0000000..22dbb3d --- /dev/null +++ b/circuitpython/lib/libm/asinfacosf.c @@ -0,0 +1,130 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// asinf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* 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. + * ==================================================== + */ + +#include "libm.h" + +// dpgeorge: pio2 was double in original implementation of asinf +static const float +pio2_hi = 1.5707962513e+00f, /* 0x3fc90fda */ +pio2_lo = 7.5497894159e-08f; /* 0x33a22168 */ + +static const float +/* coefficients for R(x^2) */ +pS0 = 1.6666586697e-01f, +pS1 = -4.2743422091e-02f, +pS2 = -8.6563630030e-03f, +qS1 = -7.0662963390e-01f; + +static float R(float z) +{ + float_t p, q; + p = z*(pS0+z*(pS1+z*pS2)); + q = 1.0f+z*qS1; + return p/q; +} + +float asinf(float x) +{ + // dpgeorge: s was double in original implementation + float s,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_hi + 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 = sqrtf(z); + x = pio2_hi - (2*(s+s*R(z)) - pio2_lo); // dpgeorge: use pio2_hi and pio2_lo + if (hx >> 31) + return -x; + return x; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// acosf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* 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. + * ==================================================== + */ + +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); +} diff --git a/circuitpython/lib/libm/asinhf.c b/circuitpython/lib/libm/asinhf.c new file mode 100644 index 0000000..4bcb3f9 --- /dev/null +++ b/circuitpython/lib/libm/asinhf.c @@ -0,0 +1,34 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// asinhf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +#include "libm.h" + +/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */ +float asinhf(float x) +{ + union {float f; uint32_t i;} u = {.f = x}; + uint32_t i = u.i & 0x7fffffff; + unsigned s = u.i >> 31; + + /* |x| */ + u.i = i; + x = u.f; + + if (i >= 0x3f800000 + (12<<23)) { + /* |x| >= 0x1p12 or inf or nan */ + x = logf(x) + 0.693147180559945309417232121458176568f; + } else if (i >= 0x3f800000 + (1<<23)) { + /* |x| >= 2 */ + x = logf(2*x + 1/(sqrtf(x*x+1)+x)); + } else if (i >= 0x3f800000 - (12<<23)) { + /* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */ + x = log1pf(x + x*x/(sqrtf(x*x+1)+1)); + } else { + /* |x| < 0x1p-12, raise inexact if x!=0 */ + FORCE_EVAL(x + 0x1p120f); + } + return s ? -x : x; +} diff --git a/circuitpython/lib/libm/atan2f.c b/circuitpython/lib/libm/atan2f.c new file mode 100644 index 0000000..dcbaf82 --- /dev/null +++ b/circuitpython/lib/libm/atan2f.c @@ -0,0 +1,89 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// atan2f from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2f.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. + * ==================================================== + */ + +#include "libm.h" + +static const float +pi = 3.1415927410e+00f, /* 0x40490fdb */ +pi_lo = -8.7422776573e-08f; /* 0xb3bbbd2e */ + +float atan2f(float y, float x) +{ + float z; + uint32_t m,ix,iy; + + if (isnan(x) || isnan(y)) + return x+y; + GET_FLOAT_WORD(ix, x); + GET_FLOAT_WORD(iy, y); + if (ix == 0x3f800000) /* x=1.0 */ + return atanf(y); + m = ((iy>>31)&1) | ((ix>>30)&2); /* 2*sign(x)+sign(y) */ + ix &= 0x7fffffff; + iy &= 0x7fffffff; + + /* when y = 0 */ + if (iy == 0) { + switch (m) { + case 0: + case 1: return y; /* atan(+-0,+anything)=+-0 */ + case 2: return pi; /* atan(+0,-anything) = pi */ + case 3: return -pi; /* atan(-0,-anything) =-pi */ + } + } + /* when x = 0 */ + if (ix == 0) + return m&1 ? -pi/2 : pi/2; + /* when x is INF */ + if (ix == 0x7f800000) { + if (iy == 0x7f800000) { + switch (m) { + case 0: return pi/4; /* atan(+INF,+INF) */ + case 1: return -pi/4; /* atan(-INF,+INF) */ + case 2: return 3*pi/4; /*atan(+INF,-INF)*/ + case 3: return -3*pi/4; /*atan(-INF,-INF)*/ + } + } else { + switch (m) { + case 0: return 0.0f; /* atan(+...,+INF) */ + case 1: return -0.0f; /* atan(-...,+INF) */ + case 2: return pi; /* atan(+...,-INF) */ + case 3: return -pi; /* atan(-...,-INF) */ + } + } + } + /* |y/x| > 0x1p26 */ + if (ix+(26<<23) < iy || iy == 0x7f800000) + return m&1 ? -pi/2 : pi/2; + + /* z = atan(|y/x|) with correct underflow */ + if ((m&2) && iy+(26<<23) < ix) /*|y/x| < 0x1p-26, x < 0 */ + z = 0.0; + else + z = atanf(fabsf(y/x)); + switch (m) { + case 0: return z; /* atan(+,+) */ + case 1: return -z; /* atan(-,+) */ + case 2: return pi - (z-pi_lo); /* atan(+,-) */ + default: /* case 3 */ + return (z-pi_lo) - pi; /* atan(-,-) */ + } +} diff --git a/circuitpython/lib/libm/atanf.c b/circuitpython/lib/libm/atanf.c new file mode 100644 index 0000000..40bc363 --- /dev/null +++ b/circuitpython/lib/libm/atanf.c @@ -0,0 +1,100 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// atanf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* 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. + * ==================================================== + */ + + +#include "libm.h" + +static const float atanhi[] = { + 4.6364760399e-01f, /* atan(0.5)hi 0x3eed6338 */ + 7.8539812565e-01f, /* atan(1.0)hi 0x3f490fda */ + 9.8279368877e-01f, /* atan(1.5)hi 0x3f7b985e */ + 1.5707962513e+00f, /* atan(inf)hi 0x3fc90fda */ +}; + +static const float atanlo[] = { + 5.0121582440e-09f, /* atan(0.5)lo 0x31ac3769 */ + 3.7748947079e-08f, /* atan(1.0)lo 0x33222168 */ + 3.4473217170e-08f, /* atan(1.5)lo 0x33140fb4 */ + 7.5497894159e-08f, /* atan(inf)lo 0x33a22168 */ +}; + +static const float aT[] = { + 3.3333328366e-01f, + -1.9999158382e-01f, + 1.4253635705e-01f, + -1.0648017377e-01f, + 6.1687607318e-02f, +}; + +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; +} diff --git a/circuitpython/lib/libm/atanhf.c b/circuitpython/lib/libm/atanhf.c new file mode 100644 index 0000000..6f95f49 --- /dev/null +++ b/circuitpython/lib/libm/atanhf.c @@ -0,0 +1,34 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// atanhf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +#include "libm.h" + +/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ +float atanhf(float x) +{ + union {float f; uint32_t i;} u = {.f = x}; + unsigned s = u.i >> 31; + float_t y; + + /* |x| */ + u.i &= 0x7fffffff; + y = u.f; + + if (u.i < 0x3f800000 - (1<<23)) { + if (u.i < 0x3f800000 - (32<<23)) { + /* handle underflow */ + if (u.i < (1<<23)) + FORCE_EVAL((float)(y*y)); + } else { + /* |x| < 0.5, up to 1.7ulp error */ + y = 0.5f*log1pf(2*y + 2*y*y/(1-y)); + } + } else { + /* avoid overflow */ + y = 0.5f*log1pf(2*(y/(1-y))); + } + return s ? -y : y; +} diff --git a/circuitpython/lib/libm/ef_rem_pio2.c b/circuitpython/lib/libm/ef_rem_pio2.c new file mode 100644 index 0000000..386c6fd --- /dev/null +++ b/circuitpython/lib/libm/ef_rem_pio2.c @@ -0,0 +1,211 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* ef_rem_pio2.c -- float version of e_rem_pio2.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. + * ==================================================== + * + */ + +/* __ieee754_rem_pio2f(x,y) + * + * return the remainder of x rem pi/2 in y[0]+y[1] + * use __kernel_rem_pio2f() + */ + +#include "fdlibm.h" + +/* + * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi + */ +#ifdef __STDC__ +static const __uint8_t two_over_pi[] = { +#else +static __uint8_t two_over_pi[] = { +#endif +0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC, +0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62, +0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63, +0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A, +0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09, +0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29, +0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44, +0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41, +0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C, +0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8, +0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11, +0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF, +0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E, +0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5, +0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92, +0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08, +0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0, +0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3, +0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85, +0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80, +0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA, +0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B, +}; + +/* This array is like the one in e_rem_pio2.c, but the numbers are + single precision and the last 8 bits are forced to 0. */ +#ifdef __STDC__ +static const __int32_t npio2_hw[] = { +#else +static __int32_t npio2_hw[] = { +#endif +0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00, +0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00, +0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100, +0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00, +0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00, +0x4242c700, 0x42490f00 +}; + +/* + * invpio2: 24 bits of 2/pi + * pio2_1: first 17 bit of pi/2 + * pio2_1t: pi/2 - pio2_1 + * pio2_2: second 17 bit of pi/2 + * pio2_2t: pi/2 - (pio2_1+pio2_2) + * pio2_3: third 17 bit of pi/2 + * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) + */ + +#ifdef __STDC__ +static const float +#else +static float +#endif +zero = 0.0000000000e+00f, /* 0x00000000 */ +half = 5.0000000000e-01f, /* 0x3f000000 */ +two8 = 2.5600000000e+02f, /* 0x43800000 */ +invpio2 = 6.3661980629e-01f, /* 0x3f22f984 */ +pio2_1 = 1.5707855225e+00f, /* 0x3fc90f80 */ +pio2_1t = 1.0804334124e-05f, /* 0x37354443 */ +pio2_2 = 1.0804273188e-05f, /* 0x37354400 */ +pio2_2t = 6.0770999344e-11f, /* 0x2e85a308 */ +pio2_3 = 6.0770943833e-11f, /* 0x2e85a300 */ +pio2_3t = 6.1232342629e-17f; /* 0x248d3132 */ + +#ifdef __STDC__ + __int32_t __ieee754_rem_pio2f(float x, float *y) +#else + __int32_t __ieee754_rem_pio2f(x,y) + float x,y[]; +#endif +{ + float z,w,t,r,fn; + float tx[3]; + __int32_t i,j,n,ix,hx; + int e0,nx; + + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; + if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */ + {y[0] = x; y[1] = 0; return 0;} + if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */ + if(hx>0) { + z = x - pio2_1; + if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ + y[0] = z - pio2_1t; + y[1] = (z-y[0])-pio2_1t; + } else { /* near pi/2, use 24+24+24 bit pi */ + z -= pio2_2; + y[0] = z - pio2_2t; + y[1] = (z-y[0])-pio2_2t; + } + return 1; + } else { /* negative x */ + z = x + pio2_1; + if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ + y[0] = z + pio2_1t; + y[1] = (z-y[0])+pio2_1t; + } else { /* near pi/2, use 24+24+24 bit pi */ + z += pio2_2; + y[0] = z + pio2_2t; + y[1] = (z-y[0])+pio2_2t; + } + return -1; + } + } +#if CIRCUITPY_FULL_BUILD + if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */ + t = fabsf(x); + n = (__int32_t) (t*invpio2+half); + fn = (float)n; + r = t-fn*pio2_1; + w = fn*pio2_1t; /* 1st round good to 40 bit */ + if(n<32&&(__int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) { + y[0] = r-w; /* quick check no cancellation */ + } else { + __uint32_t high; + j = ix>>23; + y[0] = r-w; + GET_FLOAT_WORD(high,y[0]); + i = j-((high>>23)&0xff); + if(i>8) { /* 2nd iteration needed, good to 57 */ + t = r; + w = fn*pio2_2; + r = t-w; + w = fn*pio2_2t-((t-r)-w); + y[0] = r-w; + GET_FLOAT_WORD(high,y[0]); + i = j-((high>>23)&0xff); + if(i>25) { /* 3rd iteration need, 74 bits acc */ + t = r; /* will cover all possible cases */ + w = fn*pio2_3; + r = t-w; + w = fn*pio2_3t-((t-r)-w); + y[0] = r-w; + } + } + } + y[1] = (r-y[0])-w; + if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} + else return n; + } +#else + // Suppress "defined but not used" diagnostics + (void) j; (void) fn; (void) r; (void) t; (void) w; (void) pio2_3t; + (void) pio2_3; (void) invpio2; (void)half; (void)npio2_hw; +#endif + /* + * all other (large) arguments + */ + if(!FLT_UWORD_IS_FINITE(ix)) { + y[0]=y[1]=x-x; return 0; + } + /* set z = scalbn(|x|,ilogb(x)-7) */ + e0 = (int)((ix>>23)-134); /* e0 = ilogb(z)-7; */ + SET_FLOAT_WORD(z, ix - ((__int32_t)e0<<23)); + for(i=0;i<2;i++) { + tx[i] = (float)((__int32_t)(z)); + z = (z-tx[i])*two8; + } + tx[2] = z; + nx = 3; +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wfloat-equal" + while(tx[nx-1]==zero) nx--; /* skip zero term */ +#pragma GCC diagnostic pop + n = __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi); + if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} + return n; +} diff --git a/circuitpython/lib/libm/ef_sqrt.c b/circuitpython/lib/libm/ef_sqrt.c new file mode 100644 index 0000000..2c1ae6f --- /dev/null +++ b/circuitpython/lib/libm/ef_sqrt.c @@ -0,0 +1,102 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* ef_sqrtf.c -- float version of e_sqrt.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float one = 1.0f, tiny=1.0e-30f; +#else +static float one = 1.0f, tiny=1.0e-30f; +#endif + +// sqrtf is exactly __ieee754_sqrtf when _IEEE_LIBM defined +float sqrtf(float x) +/* +#ifdef __STDC__ + float __ieee754_sqrtf(float x) +#else + float __ieee754_sqrtf(x) + float x; +#endif +*/ +{ + float z; + __uint32_t r,hx; + __int32_t ix,s,q,m,t,i; + + GET_FLOAT_WORD(ix,x); + hx = ix&0x7fffffff; + + /* take care of Inf and NaN */ + if(!FLT_UWORD_IS_FINITE(hx)) + return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf + sqrt(-inf)=sNaN */ + /* take care of zero and -ves */ + if(FLT_UWORD_IS_ZERO(hx)) return x;/* sqrt(+-0) = +-0 */ + if(ix<0) return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ + + /* normalize x */ + m = (ix>>23); + if(FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */ + for(i=0;(ix&0x00800000L)==0;i++) ix<<=1; + m -= i-1; + } + m -= 127; /* unbias exponent */ + ix = (ix&0x007fffffL)|0x00800000L; + if(m&1) /* odd m, double x to make it even */ + ix += ix; + m >>= 1; /* m = [m/2] */ + + /* generate sqrt(x) bit by bit */ + ix += ix; + q = s = 0; /* q = sqrt(x) */ + r = 0x01000000L; /* r = moving bit from right to left */ + + while(r!=0) { + t = s+r; + if(t<=ix) { + s = t+r; + ix -= t; + q += r; + } + ix += ix; + r>>=1; + } + + /* use floating add to find out rounding direction */ + if(ix!=0) { + z = one-tiny; /* trigger inexact flag */ + if (z>=one) { + z = one+tiny; + if (z>one) + q += 2; + else + q += (q&1); + } + } + ix = (q>>1)+0x3f000000L; + ix += (m <<23); + SET_FLOAT_WORD(z,ix); + return z; +} diff --git a/circuitpython/lib/libm/erf_lgamma.c b/circuitpython/lib/libm/erf_lgamma.c new file mode 100644 index 0000000..6d675f1 --- /dev/null +++ b/circuitpython/lib/libm/erf_lgamma.c @@ -0,0 +1,256 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* erf_lgamma.c -- float version of er_lgamma.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. + * ==================================================== + * + */ + +#include "fdlibm.h" +#pragma GCC diagnostic ignored "-Wfloat-equal" + +#define __ieee754_logf logf + +#ifdef __STDC__ +static const float +#else +static float +#endif +two23= 8.3886080000e+06f, /* 0x4b000000 */ +half= 5.0000000000e-01f, /* 0x3f000000 */ +one = 1.0000000000e+00f, /* 0x3f800000 */ +pi = 3.1415927410e+00f, /* 0x40490fdb */ +a0 = 7.7215664089e-02f, /* 0x3d9e233f */ +a1 = 3.2246702909e-01f, /* 0x3ea51a66 */ +a2 = 6.7352302372e-02f, /* 0x3d89f001 */ +a3 = 2.0580807701e-02f, /* 0x3ca89915 */ +a4 = 7.3855509982e-03f, /* 0x3bf2027e */ +a5 = 2.8905137442e-03f, /* 0x3b3d6ec6 */ +a6 = 1.1927076848e-03f, /* 0x3a9c54a1 */ +a7 = 5.1006977446e-04f, /* 0x3a05b634 */ +a8 = 2.2086278477e-04f, /* 0x39679767 */ +a9 = 1.0801156895e-04f, /* 0x38e28445 */ +a10 = 2.5214456400e-05f, /* 0x37d383a2 */ +a11 = 4.4864096708e-05f, /* 0x383c2c75 */ +tc = 1.4616321325e+00f, /* 0x3fbb16c3 */ +tf = -1.2148628384e-01f, /* 0xbdf8cdcd */ +/* tt = -(tail of tf) */ +tt = 6.6971006518e-09f, /* 0x31e61c52 */ +t0 = 4.8383611441e-01f, /* 0x3ef7b95e */ +t1 = -1.4758771658e-01f, /* 0xbe17213c */ +t2 = 6.4624942839e-02f, /* 0x3d845a15 */ +t3 = -3.2788541168e-02f, /* 0xbd064d47 */ +t4 = 1.7970675603e-02f, /* 0x3c93373d */ +t5 = -1.0314224288e-02f, /* 0xbc28fcfe */ +t6 = 6.1005386524e-03f, /* 0x3bc7e707 */ +t7 = -3.6845202558e-03f, /* 0xbb7177fe */ +t8 = 2.2596477065e-03f, /* 0x3b141699 */ +t9 = -1.4034647029e-03f, /* 0xbab7f476 */ +t10 = 8.8108185446e-04f, /* 0x3a66f867 */ +t11 = -5.3859531181e-04f, /* 0xba0d3085 */ +t12 = 3.1563205994e-04f, /* 0x39a57b6b */ +t13 = -3.1275415677e-04f, /* 0xb9a3f927 */ +t14 = 3.3552918467e-04f, /* 0x39afe9f7 */ +u0 = -7.7215664089e-02f, /* 0xbd9e233f */ +u1 = 6.3282704353e-01f, /* 0x3f2200f4 */ +u2 = 1.4549225569e+00f, /* 0x3fba3ae7 */ +u3 = 9.7771751881e-01f, /* 0x3f7a4bb2 */ +u4 = 2.2896373272e-01f, /* 0x3e6a7578 */ +u5 = 1.3381091878e-02f, /* 0x3c5b3c5e */ +v1 = 2.4559779167e+00f, /* 0x401d2ebe */ +v2 = 2.1284897327e+00f, /* 0x4008392d */ +v3 = 7.6928514242e-01f, /* 0x3f44efdf */ +v4 = 1.0422264785e-01f, /* 0x3dd572af */ +v5 = 3.2170924824e-03f, /* 0x3b52d5db */ +s0 = -7.7215664089e-02f, /* 0xbd9e233f */ +s1 = 2.1498242021e-01f, /* 0x3e5c245a */ +s2 = 3.2577878237e-01f, /* 0x3ea6cc7a */ +s3 = 1.4635047317e-01f, /* 0x3e15dce6 */ +s4 = 2.6642270386e-02f, /* 0x3cda40e4 */ +s5 = 1.8402845599e-03f, /* 0x3af135b4 */ +s6 = 3.1947532989e-05f, /* 0x3805ff67 */ +r1 = 1.3920053244e+00f, /* 0x3fb22d3b */ +r2 = 7.2193557024e-01f, /* 0x3f38d0c5 */ +r3 = 1.7193385959e-01f, /* 0x3e300f6e */ +r4 = 1.8645919859e-02f, /* 0x3c98bf54 */ +r5 = 7.7794247773e-04f, /* 0x3a4beed6 */ +r6 = 7.3266842264e-06f, /* 0x36f5d7bd */ +w0 = 4.1893854737e-01f, /* 0x3ed67f1d */ +w1 = 8.3333335817e-02f, /* 0x3daaaaab */ +w2 = -2.7777778450e-03f, /* 0xbb360b61 */ +w3 = 7.9365057172e-04f, /* 0x3a500cfd */ +w4 = -5.9518753551e-04f, /* 0xba1c065c */ +w5 = 8.3633989561e-04f, /* 0x3a5b3dd2 */ +w6 = -1.6309292987e-03f; /* 0xbad5c4e8 */ + +#ifdef __STDC__ +static const float zero= 0.0000000000e+00f; +#else +static float zero= 0.0000000000e+00f; +#endif + +#ifdef __STDC__ + static float sin_pif(float x) +#else + static float sin_pif(x) + float x; +#endif +{ + float y,z; + __int32_t n,ix; + + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; + + if(ix<0x3e800000) return __kernel_sinf(pi*x,zero,0); + y = -x; /* x is assume negative */ + + /* + * argument reduction, make sure inexact flag not raised if input + * is an integer + */ + z = floorf(y); + if(z!=y) { /* inexact anyway */ + y *= (float)0.5; + y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */ + n = (__int32_t) (y*(float)4.0); + } else { + if(ix>=0x4b800000) { + y = zero; n = 0; /* y must be even */ + } else { + if(ix<0x4b000000) z = y+two23; /* exact */ + GET_FLOAT_WORD(n,z); + n &= 1; + y = n; + n<<= 2; + } + } + switch (n) { + case 0: y = __kernel_sinf(pi*y,zero,0); break; + case 1: + case 2: y = __kernel_cosf(pi*((float)0.5-y),zero); break; + case 3: + case 4: y = __kernel_sinf(pi*(one-y),zero,0); break; + case 5: + case 6: y = -__kernel_cosf(pi*(y-(float)1.5),zero); break; + default: y = __kernel_sinf(pi*(y-(float)2.0),zero,0); break; + } + return -y; +} + + +#ifdef __STDC__ + float __ieee754_lgammaf_r(float x, int *signgamp) +#else + float __ieee754_lgammaf_r(x,signgamp) + float x; int *signgamp; +#endif +{ + float t,y,z,nadj = 0.0,p,p1,p2,p3,q,r,w; + __int32_t i,hx,ix; + + GET_FLOAT_WORD(hx,x); + + /* purge off +-inf, NaN, +-0, and negative arguments */ + *signgamp = 1; + ix = hx&0x7fffffff; + if(ix>=0x7f800000) return x*x; + if(ix==0) return one/zero; + if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */ + if(hx<0) { + *signgamp = -1; + return -__ieee754_logf(-x); + } else return -__ieee754_logf(x); + } + if(hx<0) { + if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ + return one/zero; + t = sin_pif(x); + if(t==zero) return one/zero; /* -integer */ + nadj = __ieee754_logf(pi/fabsf(t*x)); + if(t<zero) *signgamp = -1; + x = -x; + } + + /* purge off 1 and 2 */ + if (ix==0x3f800000||ix==0x40000000) r = 0; + /* for x < 2.0 */ + else if(ix<0x40000000) { + if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ + r = -__ieee754_logf(x); + if(ix>=0x3f3b4a20) {y = one-x; i= 0;} + else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;} + else {y = x; i=2;} + } else { + r = zero; + if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */ + else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */ + else {y=x-one;i=2;} + } + switch(i) { + case 0: + z = y*y; + p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); + p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); + p = y*p1+p2; + r += (p-(float)0.5*y); break; + case 1: + z = y*y; + w = z*y; + p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ + p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); + p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); + p = z*p1-(tt-w*(p2+y*p3)); + r += (tf + p); break; + case 2: + p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); + p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); + r += (-(float)0.5*y + p1/p2); + } + } + else if(ix<0x41000000) { /* x < 8.0 */ + i = (__int32_t)x; + t = zero; + y = x-(float)i; + p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); + q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); + r = half*y+p/q; + z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ + switch(i) { + case 7: z *= (y+(float)6.0); /* FALLTHRU */ + case 6: z *= (y+(float)5.0); /* FALLTHRU */ + case 5: z *= (y+(float)4.0); /* FALLTHRU */ + case 4: z *= (y+(float)3.0); /* FALLTHRU */ + case 3: z *= (y+(float)2.0); /* FALLTHRU */ + r += __ieee754_logf(z); break; + } + /* 8.0 <= x < 2**58 */ + } else if (ix < 0x5c800000) { + t = __ieee754_logf(x); + z = one/x; + y = z*z; + w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); + r = (x-half)*(t-one)+w; + } else + /* 2**58 <= x <= inf */ + r = x*(__ieee754_logf(x)-one); + if(hx<0) r = nadj - r; + return r; +} diff --git a/circuitpython/lib/libm/fdlibm.h b/circuitpython/lib/libm/fdlibm.h new file mode 100644 index 0000000..91698a6 --- /dev/null +++ b/circuitpython/lib/libm/fdlibm.h @@ -0,0 +1,227 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * This file is adapted from from newlib-nano-2, the newlib/libm/common/fdlib.h, + * available from https://github.com/32bitmicro/newlib-nano-2. The main change + * is removal of anything to do with double precision. + * + * Appropriate copyright headers are reproduced below. + */ + +/* @(#)fdlibm.h 5.1 93/09/24 */ +/* + * ==================================================== + * 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. + * ==================================================== + */ + +#include <math.h> + +/* Default to XOPEN_MODE. */ +#define _XOPEN_MODE + +/* Most routines need to check whether a float is finite, infinite, or not a + number, and many need to know whether the result of an operation will + overflow. These conditions depend on whether the largest exponent is + used for NaNs & infinities, or whether it's used for finite numbers. The + macros below wrap up that kind of information: + + FLT_UWORD_IS_FINITE(X) + True if a positive float with bitmask X is finite. + + FLT_UWORD_IS_NAN(X) + True if a positive float with bitmask X is not a number. + + FLT_UWORD_IS_INFINITE(X) + True if a positive float with bitmask X is +infinity. + + FLT_UWORD_MAX + The bitmask of FLT_MAX. + + FLT_UWORD_HALF_MAX + The bitmask of FLT_MAX/2. + + FLT_UWORD_EXP_MAX + The bitmask of the largest finite exponent (129 if the largest + exponent is used for finite numbers, 128 otherwise). + + FLT_UWORD_LOG_MAX + The bitmask of log(FLT_MAX), rounded down. This value is the largest + input that can be passed to exp() without producing overflow. + + FLT_UWORD_LOG_2MAX + The bitmask of log(2*FLT_MAX), rounded down. This value is the + largest input than can be passed to cosh() without producing + overflow. + + FLT_LARGEST_EXP + The largest biased exponent that can be used for finite numbers + (255 if the largest exponent is used for finite numbers, 254 + otherwise) */ + +#ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL +#define FLT_UWORD_IS_FINITE(x) 1 +#define FLT_UWORD_IS_NAN(x) 0 +#define FLT_UWORD_IS_INFINITE(x) 0 +#define FLT_UWORD_MAX 0x7fffffff +#define FLT_UWORD_EXP_MAX 0x43010000 +#define FLT_UWORD_LOG_MAX 0x42b2d4fc +#define FLT_UWORD_LOG_2MAX 0x42b437e0 +#define HUGE ((float)0X1.FFFFFEP128) +#else +#define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L) +#define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L) +#define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L) +#define FLT_UWORD_MAX 0x7f7fffffL +#define FLT_UWORD_EXP_MAX 0x43000000 +#define FLT_UWORD_LOG_MAX 0x42b17217 +#define FLT_UWORD_LOG_2MAX 0x42b2d4fc +#define HUGE ((float)3.40282346638528860e+38) +#endif +#define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23)) +#define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23) + +/* Many routines check for zero and subnormal numbers. Such things depend + on whether the target supports denormals or not: + + FLT_UWORD_IS_ZERO(X) + True if a positive float with bitmask X is +0. Without denormals, + any float with a zero exponent is a +0 representation. With + denormals, the only +0 representation is a 0 bitmask. + + FLT_UWORD_IS_SUBNORMAL(X) + True if a non-zero positive float with bitmask X is subnormal. + (Routines should check for zeros first.) + + FLT_UWORD_MIN + The bitmask of the smallest float above +0. Call this number + REAL_FLT_MIN... + + FLT_UWORD_EXP_MIN + The bitmask of the float representation of REAL_FLT_MIN's exponent. + + FLT_UWORD_LOG_MIN + The bitmask of |log(REAL_FLT_MIN)|, rounding down. + + FLT_SMALLEST_EXP + REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported, + -22 if they are). +*/ + +#ifdef _FLT_NO_DENORMALS +#define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L) +#define FLT_UWORD_IS_SUBNORMAL(x) 0 +#define FLT_UWORD_MIN 0x00800000 +#define FLT_UWORD_EXP_MIN 0x42fc0000 +#define FLT_UWORD_LOG_MIN 0x42aeac50 +#define FLT_SMALLEST_EXP 1 +#else +#define FLT_UWORD_IS_ZERO(x) ((x)==0) +#define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L) +#define FLT_UWORD_MIN 0x00000001 +#define FLT_UWORD_EXP_MIN 0x43160000 +#define FLT_UWORD_LOG_MIN 0x42cff1b5 +#define FLT_SMALLEST_EXP -22 +#endif + +#ifdef __STDC__ +#undef __P +#define __P(p) p +#else +#define __P(p) () +#endif + +/* + * set X_TLOSS = pi*2**52, which is possibly defined in <values.h> + * (one may replace the following line by "#include <values.h>") + */ + +#define X_TLOSS 1.41484755040568800000e+16 + +/* Functions that are not documented, and are not in <math.h>. */ + +/* Undocumented float functions. */ +#ifdef _SCALB_INT +extern float scalbf __P((float, int)); +#else +extern float scalbf __P((float, float)); +#endif +extern float significandf __P((float)); + +/* ieee style elementary float functions */ +extern float __ieee754_sqrtf __P((float)); +extern float __ieee754_acosf __P((float)); +extern float __ieee754_acoshf __P((float)); +extern float __ieee754_logf __P((float)); +extern float __ieee754_atanhf __P((float)); +extern float __ieee754_asinf __P((float)); +extern float __ieee754_atan2f __P((float,float)); +extern float __ieee754_expf __P((float)); +extern float __ieee754_coshf __P((float)); +extern float __ieee754_fmodf __P((float,float)); +extern float __ieee754_powf __P((float,float)); +extern float __ieee754_lgammaf_r __P((float,int *)); +extern float __ieee754_gammaf_r __P((float,int *)); +extern float __ieee754_log10f __P((float)); +extern float __ieee754_sinhf __P((float)); +extern float __ieee754_hypotf __P((float,float)); +extern float __ieee754_j0f __P((float)); +extern float __ieee754_j1f __P((float)); +extern float __ieee754_y0f __P((float)); +extern float __ieee754_y1f __P((float)); +extern float __ieee754_jnf __P((int,float)); +extern float __ieee754_ynf __P((int,float)); +extern float __ieee754_remainderf __P((float,float)); +extern __int32_t __ieee754_rem_pio2f __P((float,float*)); +#ifdef _SCALB_INT +extern float __ieee754_scalbf __P((float,int)); +#else +extern float __ieee754_scalbf __P((float,float)); +#endif + +/* float versions of fdlibm kernel functions */ +extern float __kernel_sinf __P((float,float,int)); +extern float __kernel_cosf __P((float,float)); +extern float __kernel_tanf __P((float,float,int)); +extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const __uint8_t*)); + +/* A union which permits us to convert between a float and a 32 bit + int. */ + +typedef union +{ + float value; + __uint32_t word; +} ieee_float_shape_type; + +/* Get a 32 bit int from a float. */ + +#define GET_FLOAT_WORD(i,d) \ +do { \ + ieee_float_shape_type gf_u; \ + gf_u.value = (d); \ + (i) = gf_u.word; \ +} while (0) + +/* Set a float from a 32 bit int. */ + +#define SET_FLOAT_WORD(d,i) \ +do { \ + ieee_float_shape_type sf_u; \ + sf_u.word = (i); \ + (d) = sf_u.value; \ +} while (0) + +/* Macros to avoid undefined behaviour that can arise if the amount + of a shift is exactly equal to the size of the shifted operand. */ + +#define SAFE_LEFT_SHIFT(op,amt) \ + (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0) + +#define SAFE_RIGHT_SHIFT(op,amt) \ + (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0) diff --git a/circuitpython/lib/libm/fmodf.c b/circuitpython/lib/libm/fmodf.c new file mode 100644 index 0000000..69a9ad9 --- /dev/null +++ b/circuitpython/lib/libm/fmodf.c @@ -0,0 +1,70 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// fmodf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +#include "libm.h" + +float fmodf(float x, float y) +{ + union {float f; uint32_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>23 & 0xff; + int ey = uy.i>>23 & 0xff; + uint32_t sx = ux.i & 0x80000000; + uint32_t i; + uint32_t uxi = ux.i; + + if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) + return (x*y)/(x*y); + if (uxi<<1 <= uy.i<<1) { + if (uxi<<1 == uy.i<<1) + return 0*x; + return x; + } + + /* normalize x and y */ + if (!ex) { + for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1U >> 9; + uxi |= 1U << 23; + } + if (!ey) { + for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1U >> 9; + uy.i |= 1U << 23; + } + + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + uxi <<= 1; + } + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + for (; uxi>>23 == 0; uxi <<= 1, ex--); + + /* scale result up */ + if (ex > 0) { + uxi -= 1U << 23; + uxi |= (uint32_t)ex << 23; + } else { + uxi >>= -ex + 1; + } + uxi |= sx; + ux.i = uxi; + return ux.f; +} diff --git a/circuitpython/lib/libm/kf_cos.c b/circuitpython/lib/libm/kf_cos.c new file mode 100644 index 0000000..73b4859 --- /dev/null +++ b/circuitpython/lib/libm/kf_cos.c @@ -0,0 +1,68 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* kf_cos.c -- float version of k_cos.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +one = 1.0000000000e+00f, /* 0x3f800000 */ +C1 = 4.1666667908e-02f, /* 0x3d2aaaab */ +C2 = -1.3888889225e-03f, /* 0xbab60b61 */ +C3 = 2.4801587642e-05f, /* 0x37d00d01 */ +C4 = -2.7557314297e-07f, /* 0xb493f27c */ +C5 = 2.0875723372e-09f, /* 0x310f74f6 */ +C6 = -1.1359647598e-11f; /* 0xad47d74e */ + +#ifdef __STDC__ + float __kernel_cosf(float x, float y) +#else + float __kernel_cosf(x, y) + float x,y; +#endif +{ + float a,hz,z,r,qx; + __int32_t ix; + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; /* ix = |x|'s high word*/ + if(ix<0x32000000) { /* if x < 2**27 */ + if(((int)x)==0) return one; /* generate inexact */ + } + z = x*x; + r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); + if(ix < 0x3e99999a) /* if |x| < 0.3 */ + return one - ((float)0.5*z - (z*r - x*y)); + else { + if(ix > 0x3f480000) { /* x > 0.78125 */ + qx = (float)0.28125; + } else { + SET_FLOAT_WORD(qx,ix-0x01000000); /* x/4 */ + } + hz = (float)0.5*z-qx; + a = one-qx; + return a - (hz - (z*r-x*y)); + } +} diff --git a/circuitpython/lib/libm/kf_rem_pio2.c b/circuitpython/lib/libm/kf_rem_pio2.c new file mode 100644 index 0000000..cce9eb1 --- /dev/null +++ b/circuitpython/lib/libm/kf_rem_pio2.c @@ -0,0 +1,224 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* kf_rem_pio2.c -- float version of k_rem_pio2.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. + * ==================================================== + */ + +#include "fdlibm.h" + +/* In the float version, the input parameter x contains 8 bit + integers, not 24 bit integers. 113 bit precision is not supported. */ + +#ifdef __STDC__ +static const int init_jk[] = {4,7,9}; /* initial value for jk */ +#else +static int init_jk[] = {4,7,9}; +#endif + +#ifdef __STDC__ +static const float PIo2[] = { +#else +static float PIo2[] = { +#endif + 1.5703125000e+00f, /* 0x3fc90000 */ + 4.5776367188e-04f, /* 0x39f00000 */ + 2.5987625122e-05f, /* 0x37da0000 */ + 7.5437128544e-08f, /* 0x33a20000 */ + 6.0026650317e-11f, /* 0x2e840000 */ + 7.3896444519e-13f, /* 0x2b500000 */ + 5.3845816694e-15f, /* 0x27c20000 */ + 5.6378512969e-18f, /* 0x22d00000 */ + 8.3009228831e-20f, /* 0x1fc40000 */ + 3.2756352257e-22f, /* 0x1bc60000 */ + 6.3331015649e-25f, /* 0x17440000 */ +}; + +#ifdef __STDC__ +static const float +#else +static float +#endif +zero = 0.0f, +one = 1.0f, +two8 = 2.5600000000e+02f, /* 0x43800000 */ +twon8 = 3.9062500000e-03f; /* 0x3b800000 */ + +#ifdef __STDC__ + int __kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec, const __uint8_t *ipio2) +#else + int __kernel_rem_pio2f(x,y,e0,nx,prec,ipio2) + float x[], y[]; int e0,nx,prec; __uint8_t ipio2[]; +#endif +{ + __int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; + float z,fw,f[20],fq[20],q[20]; + + /* initialize jk*/ + jk = init_jk[prec]; + jp = jk; + + /* determine jx,jv,q0, note that 3>q0 */ + jx = nx-1; + jv = (e0-3)/8; if(jv<0) jv=0; + q0 = e0-8*(jv+1); + + /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ + j = jv-jx; m = jx+jk; + for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (float) ipio2[j]; + + /* compute q[0],q[1],...q[jk] */ + for (i=0;i<=jk;i++) { + for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; + q[i] = fw; + } + + jz = jk; +recompute: + /* distill q[] into iq[] reversingly */ + for(i=0,j=jz,z=q[jz];j>0;i++,j--) { + fw = (float)((__int32_t)(twon8* z)); + iq[i] = (__int32_t)(z-two8*fw); + z = q[j-1]+fw; + } + + /* compute n */ + z = scalbnf(z,(int)q0); /* actual value of z */ + z -= (float)8.0*floorf(z*(float)0.125); /* trim off integer >= 8 */ + n = (__int32_t) z; + z -= (float)n; + ih = 0; + if(q0>0) { /* need iq[jz-1] to determine n */ + i = (iq[jz-1]>>(8-q0)); n += i; + iq[jz-1] -= i<<(8-q0); + ih = iq[jz-1]>>(7-q0); + } + else if(q0==0) ih = iq[jz-1]>>8; + else if(z>=(float)0.5) ih=2; + + if(ih>0) { /* q > 0.5 */ + n += 1; carry = 0; + for(i=0;i<jz ;i++) { /* compute 1-q */ + j = iq[i]; + if(carry==0) { + if(j!=0) { + carry = 1; iq[i] = 0x100- j; + } + } else iq[i] = 0xff - j; + } + if(q0>0) { /* rare case: chance is 1 in 12 */ + switch(q0) { + case 1: + iq[jz-1] &= 0x7f; break; + case 2: + iq[jz-1] &= 0x3f; break; + } + } + if(ih==2) { + z = one - z; + if(carry!=0) z -= scalbnf(one,(int)q0); + } + } + + /* check if recomputation is needed */ +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wfloat-equal" + if(z==zero) { +#pragma GCC diagnostic pop + j = 0; + for (i=jz-1;i>=jk;i--) j |= iq[i]; + if(j==0) { /* need recomputation */ + for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ + + for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ + f[jx+i] = (float) ipio2[jv+i]; + for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; + q[i] = fw; + } + jz += k; + goto recompute; + } + } + + /* chop off zero terms */ +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wfloat-equal" + if(z==(float)0.0) { +#pragma GCC diagnostic pop + jz -= 1; q0 -= 8; + while(iq[jz]==0) { jz--; q0-=8;} + } else { /* break z into 8-bit if necessary */ + z = scalbnf(z,-(int)q0); + if(z>=two8) { + fw = (float)((__int32_t)(twon8*z)); + iq[jz] = (__int32_t)(z-two8*fw); + jz += 1; q0 += 8; + iq[jz] = (__int32_t) fw; + } else iq[jz] = (__int32_t) z ; + } + + /* convert integer "bit" chunk to floating-point value */ + fw = scalbnf(one,(int)q0); + for(i=jz;i>=0;i--) { + q[i] = fw*(float)iq[i]; fw*=twon8; + } + + /* compute PIo2[0,...,jp]*q[jz,...,0] */ + for(i=jz;i>=0;i--) { + for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; + fq[jz-i] = fw; + } + + /* compress fq[] into y[] */ + switch(prec) { + case 0: + fw = 0.0; + for (i=jz;i>=0;i--) fw += fq[i]; + y[0] = (ih==0)? fw: -fw; + break; + case 1: + case 2: + fw = 0.0; + for (i=jz;i>=0;i--) fw += fq[i]; + y[0] = (ih==0)? fw: -fw; + fw = fq[0]-fw; + for (i=1;i<=jz;i++) fw += fq[i]; + y[1] = (ih==0)? fw: -fw; + break; + case 3: /* painful */ + for (i=jz;i>0;i--) { + fw = fq[i-1]+fq[i]; + fq[i] += fq[i-1]-fw; + fq[i-1] = fw; + } + for (i=jz;i>1;i--) { + fw = fq[i-1]+fq[i]; + fq[i] += fq[i-1]-fw; + fq[i-1] = fw; + } + for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; + if(ih==0) { + y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; + } else { + y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; + } + } + return n&7; +} diff --git a/circuitpython/lib/libm/kf_sin.c b/circuitpython/lib/libm/kf_sin.c new file mode 100644 index 0000000..9f22906 --- /dev/null +++ b/circuitpython/lib/libm/kf_sin.c @@ -0,0 +1,58 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* kf_sin.c -- float version of k_sin.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +half = 5.0000000000e-01f,/* 0x3f000000 */ +S1 = -1.6666667163e-01f, /* 0xbe2aaaab */ +S2 = 8.3333337680e-03f, /* 0x3c088889 */ +S3 = -1.9841270114e-04f, /* 0xb9500d01 */ +S4 = 2.7557314297e-06f, /* 0x3638ef1b */ +S5 = -2.5050759689e-08f, /* 0xb2d72f34 */ +S6 = 1.5896910177e-10f; /* 0x2f2ec9d3 */ + +#ifdef __STDC__ + float __kernel_sinf(float x, float y, int iy) +#else + float __kernel_sinf(x, y, iy) + float x,y; int iy; /* iy=0 if y is zero */ +#endif +{ + float z,r,v; + __int32_t ix; + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; /* high word of x */ + if(ix<0x32000000) /* |x| < 2**-27 */ + {if((int)x==0) return x;} /* generate inexact */ + z = x*x; + v = z*x; + r = S2+z*(S3+z*(S4+z*(S5+z*S6))); + if(iy==0) return x+v*(S1+z*r); + else return x-((z*(half*y-v*r)-y)-v*S1); +} diff --git a/circuitpython/lib/libm/kf_tan.c b/circuitpython/lib/libm/kf_tan.c new file mode 100644 index 0000000..305645b --- /dev/null +++ b/circuitpython/lib/libm/kf_tan.c @@ -0,0 +1,105 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* kf_tan.c -- float version of k_tan.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. + * ==================================================== + */ + +#include "libm.h" +#ifdef __STDC__ +static const float +#else +static float +#endif +one = 1.0000000000e+00f, /* 0x3f800000 */ +pio4 = 7.8539812565e-01f, /* 0x3f490fda */ +pio4lo= 3.7748947079e-08f, /* 0x33222168 */ +T[] = { + 3.3333334327e-01f, /* 0x3eaaaaab */ + 1.3333334029e-01f, /* 0x3e088889 */ + 5.3968254477e-02f, /* 0x3d5d0dd1 */ + 2.1869488060e-02f, /* 0x3cb327a4 */ + 8.8632395491e-03f, /* 0x3c11371f */ + 3.5920790397e-03f, /* 0x3b6b6916 */ + 1.4562094584e-03f, /* 0x3abede48 */ + 5.8804126456e-04f, /* 0x3a1a26c8 */ + 2.4646313977e-04f, /* 0x398137b9 */ + 7.8179444245e-05f, /* 0x38a3f445 */ + 7.1407252108e-05f, /* 0x3895c07a */ + -1.8558637748e-05f, /* 0xb79bae5f */ + 2.5907305826e-05f, /* 0x37d95384 */ +}; + +#ifdef __STDC__ + float __kernel_tanf(float x, float y, int iy) +#else + float __kernel_tanf(x, y, iy) + float x,y; int iy; +#endif +{ + float z,r,v,w,s; + __int32_t ix,hx; + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; /* high word of |x| */ + if(ix<0x31800000) /* x < 2**-28 */ + {if((int)x==0) { /* generate inexact */ + if((ix|(iy+1))==0) return one/fabsf(x); + else return (iy==1)? x: -one/x; + } + } + if(ix>=0x3f2ca140) { /* |x|>=0.6744 */ + if(hx<0) {x = -x; y = -y;} + z = pio4-x; + w = pio4lo-y; + x = z+w; y = 0.0; + } + z = x*x; + w = z*z; + /* Break x^5*(T[1]+x^2*T[2]+...) into + * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + + * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12])) + */ + r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11])))); + v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12]))))); + s = z*x; + r = y + z*(s*(r+v)+y); + r += T[0]*s; + w = x+r; + if(ix>=0x3f2ca140) { + v = (float)iy; + return (float)(1-((hx>>30)&2))*(v-(float)2.0*(x-(w*w/(w+v)-r))); + } + if(iy==1) return w; + else { /* if allow error up to 2 ulp, + simply return -1.0/(x+r) here */ + /* compute -1.0/(x+r) accurately */ + float a,t; + __int32_t i; + z = w; + GET_FLOAT_WORD(i,z); + SET_FLOAT_WORD(z,i&0xfffff000); + v = r-(z - x); /* z+v = r+x */ + t = a = -(float)1.0/w; /* a = -1.0/w */ + GET_FLOAT_WORD(i,t); + SET_FLOAT_WORD(t,i&0xfffff000); + s = (float)1.0+t*z; + return t+a*(s+t*v); + } +} diff --git a/circuitpython/lib/libm/libm.h b/circuitpython/lib/libm/libm.h new file mode 100644 index 0000000..f782249 --- /dev/null +++ b/circuitpython/lib/libm/libm.h @@ -0,0 +1,54 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// portions extracted from musl-0.9.15 libm.h +/*****************************************************************************/ +/*****************************************************************************/ + +/* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */ +/* + * ==================================================== + * 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. + * ==================================================== + */ + +#include <stdint.h> +#include <math.h> + +#define FLT_EVAL_METHOD 0 + +#define FORCE_EVAL(x) do { \ + if (sizeof(x) == sizeof(float)) { \ + volatile float __x; \ + __x = (x); \ + (void)__x; \ + } else if (sizeof(x) == sizeof(double)) { \ + volatile double __x; \ + __x = (x); \ + (void)__x; \ + } else { \ + volatile long double __x; \ + __x = (x); \ + (void)__x; \ + } \ +} while(0) + +/* Get a 32 bit int from a float. */ +#define GET_FLOAT_WORD(w,d) \ +do { \ + union {float f; uint32_t i;} __u; \ + __u.f = (d); \ + (w) = __u.i; \ +} while (0) + +/* Set a float from a 32 bit int. */ +#define SET_FLOAT_WORD(d,w) \ +do { \ + union {float f; uint32_t i;} __u; \ + __u.i = (w); \ + (d) = __u.f; \ +} while (0) diff --git a/circuitpython/lib/libm/log1pf.c b/circuitpython/lib/libm/log1pf.c new file mode 100644 index 0000000..306b1cd --- /dev/null +++ b/circuitpython/lib/libm/log1pf.c @@ -0,0 +1,85 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// log1pf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.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. + * ==================================================== + */ + +#include "libm.h" + +#pragma GCC diagnostic ignored "-Wfloat-equal" + +static const float +ln2_hi = 6.9313812256e-01f, /* 0x3f317180 */ +ln2_lo = 9.0580006145e-06f, /* 0x3717f7d1 */ +/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ +Lg1 = 0xaaaaaa.0p-24f, /* 0.66666662693 */ +Lg2 = 0xccce13.0p-25f, /* 0.40000972152 */ +Lg3 = 0x91e9ee.0p-25f, /* 0.28498786688 */ +Lg4 = 0xf89e26.0p-26f; /* 0.24279078841 */ + +float log1pf(float x) +{ + union {float f; uint32_t i;} u = {x}; + float_t hfsq,f,c,s,z,R,w,t1,t2,dk; + uint32_t ix,iu; + int k; + + ix = u.i; + k = 1; + if (ix < 0x3ed413d0 || ix>>31) { /* 1+x < sqrt(2)+ */ + if (ix >= 0xbf800000) { /* x <= -1.0 */ + if (x == -1) + return x/0.0f; /* log1p(-1)=+inf */ + return (x-x)/0.0f; /* log1p(x<-1)=NaN */ + } + if (ix<<1 < 0x33800000<<1) { /* |x| < 2**-24 */ + /* underflow if subnormal */ + if ((ix&0x7f800000) == 0) + FORCE_EVAL(x*x); + return x; + } + if (ix <= 0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ + k = 0; + c = 0; + f = x; + } + } else if (ix >= 0x7f800000) + return x; + if (k) { + u.f = 1 + x; + iu = u.i; + iu += 0x3f800000 - 0x3f3504f3; + k = (int)(iu>>23) - 0x7f; + /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ + if (k < 25) { + c = k >= 2 ? 1-(u.f-x) : x-(u.f-1); + c /= u.f; + } else + c = 0; + /* reduce u into [sqrt(2)/2, sqrt(2)] */ + iu = (iu&0x007fffff) + 0x3f3504f3; + u.i = iu; + f = u.f - 1; + } + s = f/(2.0f + f); + z = s*s; + w = z*z; + t1= w*(Lg2+w*Lg4); + t2= z*(Lg1+w*Lg3); + R = t2 + t1; + hfsq = 0.5f*f*f; + dk = k; + return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi; +} diff --git a/circuitpython/lib/libm/math.c b/circuitpython/lib/libm/math.c new file mode 100644 index 0000000..09f4fc2 --- /dev/null +++ b/circuitpython/lib/libm/math.c @@ -0,0 +1,829 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * The MIT License (MIT) + * + * SPDX-FileCopyrightText: Copyright (c) 2013, 2014 Damien P. George + * + * 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 "libm.h" + +typedef float float_t; +typedef union { + float f; + struct { + uint32_t m : 23; + uint32_t e : 8; + uint32_t s : 1; + }; +} float_s_t; + +int __signbitf(float f) { + float_s_t u = {.f = f}; + return u.s; +} + +#ifndef NDEBUG +float copysignf(float x, float y) { + float_s_t fx={.f = x}; + float_s_t fy={.f = y}; + + // copy sign bit; + fx.s = fy.s; + + return fx.f; +} +#endif + +static const float _M_LN10 = 2.30258509299404f; // 0x40135d8e +float log10f(float x) { return logf(x) / (float)_M_LN10; } +#undef _M_LN2 +static const float _M_LN2 = 0.6931472; +float log2f(float x) { return logf(x) / (float)_M_LN2; } + +float tanhf(float x) { + int sign = 0; + if (x < 0) { + sign = 1; + x = -x; + } + x = expm1f(-2 * x); + x = x / (x + 2); + return sign ? x : -x; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// __fpclassifyf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +int __fpclassifyf(float x) +{ + union {float f; uint32_t i;} u = {x}; + int e = u.i>>23 & 0xff; + if (!e) return u.i<<1 ? FP_SUBNORMAL : FP_ZERO; + if (e==0xff) return u.i<<9 ? FP_NAN : FP_INFINITE; + return FP_NORMAL; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// scalbnf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +float scalbnf(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; + n += 126; + if (n < -126) { + y *= 0x1p-126f; + n += 126; + if (n < -126) + n = -126; + } + } + u.i = (uint32_t)(0x7f+n)<<23; + x = y * u.f; + return x; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// powf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* origin: FreeBSD /usr/src/lib/msun/src/e_powf.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 +bp[] = {1.0f, 1.5f,}, +dp_h[] = { 0.0f, 5.84960938e-01f,}, /* 0x3f15c000 */ +dp_l[] = { 0.0f, 1.56322085e-06f,}, /* 0x35d1cfdc */ +two24 = 16777216.0f, /* 0x4b800000 */ +huge = 1.0e30f, +tiny = 1.0e-30f, +/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ +L1 = 6.0000002384e-01f, /* 0x3f19999a */ +L2 = 4.2857143283e-01f, /* 0x3edb6db7 */ +L3 = 3.3333334327e-01f, /* 0x3eaaaaab */ +L4 = 2.7272811532e-01f, /* 0x3e8ba305 */ +L5 = 2.3066075146e-01f, /* 0x3e6c3255 */ +L6 = 2.0697501302e-01f, /* 0x3e53f142 */ +P1 = 1.6666667163e-01f, /* 0x3e2aaaab */ +P2 = -2.7777778450e-03f, /* 0xbb360b61 */ +P3 = 6.6137559770e-05f, /* 0x388ab355 */ +P4 = -1.6533901999e-06f, /* 0xb5ddea0e */ +P5 = 4.1381369442e-08f, /* 0x3331bb4c */ +lg2 = 6.9314718246e-01f, /* 0x3f317218 */ +lg2_h = 6.93145752e-01f, /* 0x3f317200 */ +lg2_l = 1.42860654e-06f, /* 0x35bfbe8c */ +ovt = 4.2995665694e-08f, /* -(128-log2(ovfl+.5ulp)) */ +cp = 9.6179670095e-01f, /* 0x3f76384f =2/(3ln2) */ +cp_h = 9.6191406250e-01f, /* 0x3f764000 =12b cp */ +cp_l = -1.1736857402e-04f, /* 0xb8f623c6 =tail of cp_h */ +ivln2 = 1.4426950216e+00f, /* 0x3fb8aa3b =1/ln2 */ +ivln2_h = 1.4426879883e+00f, /* 0x3fb8aa00 =16b 1/ln2*/ +ivln2_l = 7.0526075433e-06f; /* 0x36eca570 =1/ln2 tail*/ + +float powf(float x, float y) +{ + float z,ax,z_h,z_l,p_h,p_l; + float y1,t1,t2,r,s,sn,t,u,v,w; + int32_t i,j,k,yisint,n; + int32_t hx,hy,ix,iy,is; + + GET_FLOAT_WORD(hx, x); + GET_FLOAT_WORD(hy, y); + ix = hx & 0x7fffffff; + iy = hy & 0x7fffffff; + + /* x**0 = 1, even if x is NaN */ + if (iy == 0) + return 1.0f; + /* 1**y = 1, even if y is NaN */ + if (hx == 0x3f800000) + return 1.0f; + /* NaN if either arg is NaN */ + if (ix > 0x7f800000 || iy > 0x7f800000) + return x + y; + + /* determine if y is an odd int when x < 0 + * yisint = 0 ... y is not an integer + * yisint = 1 ... y is an odd int + * yisint = 2 ... y is an even int + */ + yisint = 0; + if (hx < 0) { + if (iy >= 0x4b800000) + yisint = 2; /* even integer y */ + else if (iy >= 0x3f800000) { + k = (iy>>23) - 0x7f; /* exponent */ + j = iy>>(23-k); + if ((j<<(23-k)) == iy) + yisint = 2 - (j & 1); + } + } + + /* special value of y */ + if (iy == 0x7f800000) { /* y is +-inf */ + if (ix == 0x3f800000) /* (-1)**+-inf is 1 */ + return 1.0f; + else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */ + return hy >= 0 ? y : 0.0f; + else if (ix != 0) /* (|x|<1)**+-inf = 0,inf if x!=0 */ + return hy >= 0 ? 0.0f: -y; + } + if (iy == 0x3f800000) /* y is +-1 */ + return hy >= 0 ? x : 1.0f/x; + if (hy == 0x40000000) /* y is 2 */ + return x*x; + if (hy == 0x3f000000) { /* y is 0.5 */ + if (hx >= 0) /* x >= +0 */ + return sqrtf(x); + } + + ax = fabsf(x); + /* special value of x */ + if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */ + z = ax; + if (hy < 0) /* z = (1/|x|) */ + z = 1.0f/z; + if (hx < 0) { + if (((ix-0x3f800000)|yisint) == 0) { + z = (z-z)/(z-z); /* (-1)**non-int is NaN */ + } else if (yisint == 1) + z = -z; /* (x<0)**odd = -(|x|**odd) */ + } + return z; + } + + sn = 1.0f; /* sign of result */ + if (hx < 0) { + if (yisint == 0) /* (x<0)**(non-int) is NaN */ + return (x-x)/(x-x); + if (yisint == 1) /* (x<0)**(odd int) */ + sn = -1.0f; + } + + /* |y| is huge */ + if (iy > 0x4d000000) { /* if |y| > 2**27 */ + /* over/underflow if x is not close to one */ + if (ix < 0x3f7ffff8) + return hy < 0 ? sn*huge*huge : sn*tiny*tiny; + if (ix > 0x3f800007) + return hy > 0 ? sn*huge*huge : sn*tiny*tiny; + /* now |1-x| is tiny <= 2**-20, suffice to compute + log(x) by x-x^2/2+x^3/3-x^4/4 */ + t = ax - 1; /* t has 20 trailing zeros */ + w = (t*t)*(0.5f - t*(0.333333333333f - t*0.25f)); + u = ivln2_h*t; /* ivln2_h has 16 sig. bits */ + v = t*ivln2_l - w*ivln2; + t1 = u + v; + GET_FLOAT_WORD(is, t1); + SET_FLOAT_WORD(t1, is & 0xfffff000); + t2 = v - (t1-u); + } else { + float s2,s_h,s_l,t_h,t_l; + n = 0; + /* take care subnormal number */ + if (ix < 0x00800000) { + ax *= two24; + n -= 24; + GET_FLOAT_WORD(ix, ax); + } + n += ((ix)>>23) - 0x7f; + j = ix & 0x007fffff; + /* determine interval */ + ix = j | 0x3f800000; /* normalize ix */ + if (j <= 0x1cc471) /* |x|<sqrt(3/2) */ + k = 0; + else if (j < 0x5db3d7) /* |x|<sqrt(3) */ + k = 1; + else { + k = 0; + n += 1; + ix -= 0x00800000; + } + SET_FLOAT_WORD(ax, ix); + + /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ + u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ + v = 1.0f/(ax+bp[k]); + s = u*v; + s_h = s; + GET_FLOAT_WORD(is, s_h); + SET_FLOAT_WORD(s_h, is & 0xfffff000); + /* t_h=ax+bp[k] High */ + is = ((ix>>1) & 0xfffff000) | 0x20000000; + SET_FLOAT_WORD(t_h, is + 0x00400000 + (k<<21)); + t_l = ax - (t_h - bp[k]); + s_l = v*((u - s_h*t_h) - s_h*t_l); + /* compute log(ax) */ + s2 = s*s; + r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); + r += s_l*(s_h+s); + s2 = s_h*s_h; + t_h = 3.0f + s2 + r; + GET_FLOAT_WORD(is, t_h); + SET_FLOAT_WORD(t_h, is & 0xfffff000); + t_l = r - ((t_h - 3.0f) - s2); + /* u+v = s*(1+...) */ + u = s_h*t_h; + v = s_l*t_h + t_l*s; + /* 2/(3log2)*(s+...) */ + p_h = u + v; + GET_FLOAT_WORD(is, p_h); + SET_FLOAT_WORD(p_h, is & 0xfffff000); + p_l = v - (p_h - u); + z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ + z_l = cp_l*p_h + p_l*cp+dp_l[k]; + /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ + t = (float)n; + t1 = (((z_h + z_l) + dp_h[k]) + t); + GET_FLOAT_WORD(is, t1); + SET_FLOAT_WORD(t1, is & 0xfffff000); + t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); + } + + /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ + GET_FLOAT_WORD(is, y); + SET_FLOAT_WORD(y1, is & 0xfffff000); + p_l = (y-y1)*t1 + y*t2; + p_h = y1*t1; + z = p_l + p_h; + GET_FLOAT_WORD(j, z); + if (j > 0x43000000) /* if z > 128 */ + return sn*huge*huge; /* overflow */ + else if (j == 0x43000000) { /* if z == 128 */ + if (p_l + ovt > z - p_h) + return sn*huge*huge; /* overflow */ + } else if ((j&0x7fffffff) > 0x43160000) /* z < -150 */ // FIXME: check should be (uint32_t)j > 0xc3160000 + return sn*tiny*tiny; /* underflow */ + else if (j == (int32_t) 0xc3160000) { /* z == -150 */ + if (p_l <= z-p_h) + return sn*tiny*tiny; /* underflow */ + } + /* + * compute 2**(p_h+p_l) + */ + i = j & 0x7fffffff; + k = (i>>23) - 0x7f; + n = 0; + if (i > 0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */ + n = j + (0x00800000>>(k+1)); + k = ((n&0x7fffffff)>>23) - 0x7f; /* new k for n */ + SET_FLOAT_WORD(t, n & ~(0x007fffff>>k)); + n = ((n&0x007fffff)|0x00800000)>>(23-k); + if (j < 0) + n = -n; + p_h -= t; + } + t = p_l + p_h; + GET_FLOAT_WORD(is, t); + SET_FLOAT_WORD(t, is & 0xffff8000); + u = t*lg2_h; + v = (p_l-(t-p_h))*lg2 + t*lg2_l; + z = u + v; + w = v - (z - u); + t = z*z; + t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); + r = (z*t1)/(t1-2.0f) - (w+z*w); + z = 1.0f - (r - z); + GET_FLOAT_WORD(j, z); + j += n<<23; + if ((j>>23) <= 0) /* subnormal output */ + z = scalbnf(z, n); + else + SET_FLOAT_WORD(z, j); + return sn*z; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// expf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* origin: FreeBSD /usr/src/lib/msun/src/e_expf.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 +half[2] = {0.5f,-0.5f}, +ln2hi = 6.9314575195e-1f, /* 0x3f317200 */ +ln2lo = 1.4286067653e-6f, /* 0x35bfbe8e */ +invln2 = 1.4426950216e+0f, /* 0x3fb8aa3b */ +/* + * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]: + * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74 + */ +expf_P1 = 1.6666625440e-1f, /* 0xaaaa8f.0p-26 */ +expf_P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */ + +float expf(float x) +{ + float_t hi, lo, c, xx, y; + int k, sign; + uint32_t hx; + + GET_FLOAT_WORD(hx, x); + sign = hx >> 31; /* sign bit of x */ + hx &= 0x7fffffff; /* high word of |x| */ + + /* special cases */ + if (hx >= 0x42aeac50) { /* if |x| >= -87.33655f or NaN */ + if (hx >= 0x42b17218 && !sign) { /* x >= 88.722839f */ + /* overflow */ + x *= 0x1p127f; + return x; + } + if (sign) { + /* underflow */ + FORCE_EVAL(-0x1p-149f/x); + if (hx >= 0x42cff1b5) /* x <= -103.972084f */ + return 0; + } + } + + /* argument reduction */ + if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ + if (hx > 0x3f851592) /* if |x| > 1.5 ln2 */ + k = (int)(invln2*x + half[sign]); + else + k = 1 - sign - sign; + hi = x - k*ln2hi; /* k*ln2hi is exact here */ + lo = k*ln2lo; + x = hi - lo; + } else if (hx > 0x39000000) { /* |x| > 2**-14 */ + k = 0; + hi = x; + lo = 0; + } else { + /* raise inexact */ + FORCE_EVAL(0x1p127f + x); + return 1 + x; + } + + /* x is now in primary range */ + xx = x*x; + c = x - xx*(expf_P1+xx*expf_P2); + y = 1 + (x*c/(2-c) - lo + hi); + if (k == 0) + return y; + return scalbnf(y, k); +} + +/*****************************************************************************/ +/*****************************************************************************/ +// expm1f from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.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 +o_threshold = 8.8721679688e+01f, /* 0x42b17180 */ +ln2_hi = 6.9313812256e-01f, /* 0x3f317180 */ +ln2_lo = 9.0580006145e-06f, /* 0x3717f7d1 */ +//invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */ +/* + * Domain [-0.34568, 0.34568], range ~[-6.694e-10, 6.696e-10]: + * |6 / x * (1 + 2 * (1 / (exp(x) - 1) - 1 / x)) - q(x)| < 2**-30.04 + * Scaled coefficients: Qn_here = 2**n * Qn_for_q (see s_expm1.c): + */ +Q1 = -3.3333212137e-2f, /* -0x888868.0p-28 */ +Q2 = 1.5807170421e-3f; /* 0xcf3010.0p-33 */ + +float expm1f(float x) +{ + float_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk; + union {float f; uint32_t i;} u = {x}; + uint32_t hx = u.i & 0x7fffffff; + int k, sign = u.i >> 31; + + /* filter out huge and non-finite argument */ + if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */ + if (hx > 0x7f800000) /* NaN */ + return x; + if (sign) + return -1; + if (x > o_threshold) { + x *= 0x1p127f; + return x; + } + } + + /* argument reduction */ + if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ + if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ + if (!sign) { + hi = x - ln2_hi; + lo = ln2_lo; + k = 1; + } else { + hi = x + ln2_hi; + lo = -ln2_lo; + k = -1; + } + } else { + k = (int)(invln2*x + (sign ? -0.5f : 0.5f)); + t = k; + hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ + lo = t*ln2_lo; + } + x = hi-lo; + c = (hi-x)-lo; + } else if (hx < 0x33000000) { /* when |x|<2**-25, return x */ + if (hx < 0x00800000) + FORCE_EVAL(x*x); + return x; + } else + k = 0; + + /* x is now in primary range */ + hfx = 0.5f*x; + hxs = x*hfx; + r1 = 1.0f+hxs*(Q1+hxs*Q2); + t = 3.0f - r1*hfx; + e = hxs*((r1-t)/(6.0f - x*t)); + if (k == 0) /* c is 0 */ + return x - (x*e-hxs); + e = x*(e-c) - c; + e -= hxs; + /* exp(x) ~ 2^k (x_reduced - e + 1) */ + if (k == -1) + return 0.5f*(x-e) - 0.5f; + if (k == 1) { + if (x < -0.25f) + return -2.0f*(e-(x+0.5f)); + return 1.0f + 2.0f*(x-e); + } + u.i = (0x7f+k)<<23; /* 2^k */ + twopk = u.f; + if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ + y = x - e + 1.0f; + if (k == 128) + y = y*2.0f*0x1p127f; + else + y = y*twopk; + return y - 1.0f; + } + u.i = (0x7f-k)<<23; /* 2^-k */ + if (k < 23) + y = (x-e+(1-u.f))*twopk; + else + y = (x-(e+u.f)+1)*twopk; + return y; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// __expo2f from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + + +/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */ +float __expo2f(float x) +{ + /* k is such that k*ln2 has minimal relative error and x - kln2 > log(FLT_MIN) */ + static const int k = 235; + static const float kln2 = 0x1.45c778p+7f; + float scale; + + /* note that k is odd and scale*scale overflows */ + SET_FLOAT_WORD(scale, (uint32_t)(0x7f + k/2) << 23); + /* exp(x - k ln2) * 2**(k-1) */ + return expf(x - kln2) * scale * scale; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// logf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +/* origin: FreeBSD /usr/src/lib/msun/src/e_logf.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 +/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ +Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ +Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ +Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ +Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ + +float logf(float x) +{ + union {float f; uint32_t i;} u = {x}; + float_t hfsq,f,s,z,R,w,t1,t2,dk; + uint32_t ix; + int k; + + ix = u.i; + k = 0; + if (ix < 0x00800000 || ix>>31) { /* x < 2**-126 */ + if (ix<<1 == 0) + return -1/(x*x); /* log(+-0)=-inf */ + if (ix>>31) + return (x-x)/0.0f; /* log(-#) = NaN */ + /* subnormal number, scale up x */ + k -= 25; + x *= 0x1p25f; + u.f = x; + ix = u.i; + } else if (ix >= 0x7f800000) { + return x; + } else if (ix == 0x3f800000) + return 0; + + /* reduce x into [sqrt(2)/2, sqrt(2)] */ + ix += 0x3f800000 - 0x3f3504f3; + k += (int)(ix>>23) - 0x7f; + ix = (ix&0x007fffff) + 0x3f3504f3; + u.i = ix; + x = u.f; + + f = x - 1.0f; + s = f/(2.0f + f); + z = s*s; + w = z*z; + t1= w*(Lg2+w*Lg4); + t2= z*(Lg1+w*Lg3); + R = t2 + t1; + hfsq = 0.5f*f*f; + dk = k; + return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// coshf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +float coshf(float x) +{ + union {float f; uint32_t i;} u = {.f = x}; + uint32_t w; + float t; + + /* |x| */ + u.i &= 0x7fffffff; + x = u.f; + w = u.i; + + /* |x| < log(2) */ + if (w < 0x3f317217) { + if (w < 0x3f800000 - (12<<23)) { + FORCE_EVAL(x + 0x1p120f); + return 1; + } + t = expm1f(x); + return 1 + t*t/(2*(1+t)); + } + + /* |x| < log(FLT_MAX) */ + if (w < 0x42b17217) { + t = expf(x); + return 0.5f*(t + 1/t); + } + + /* |x| > log(FLT_MAX) or nan */ + t = __expo2f(x); + return t; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// sinhf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +float sinhf(float x) +{ + union {float f; uint32_t i;} u = {.f = x}; + uint32_t w; + float t, h, absx; + + h = 0.5; + if (u.i >> 31) + h = -h; + /* |x| */ + u.i &= 0x7fffffff; + absx = u.f; + w = u.i; + + /* |x| < log(FLT_MAX) */ + if (w < 0x42b17217) { + t = expm1f(absx); + if (w < 0x3f800000) { + if (w < 0x3f800000 - (12<<23)) + return x; + return h*(2*t - t*t/(t+1)); + } + return h*(t + t/(t+1)); + } + + /* |x| > logf(FLT_MAX) or nan */ + t = 2*h*__expo2f(absx); + return t; +} + +/*****************************************************************************/ +/*****************************************************************************/ +// ceilf, floorf and truncf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +float ceilf(float x) +{ + union {float f; uint32_t i;} u = {x}; + int e = (int)(u.i >> 23 & 0xff) - 0x7f; + uint32_t m; + + if (e >= 23) + return x; + if (e >= 0) { + m = 0x007fffff >> e; + if ((u.i & m) == 0) + return x; + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31 == 0) + u.i += m; + u.i &= ~m; + } else { + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31) + u.f = -0.0; + else if (u.i << 1) + u.f = 1.0; + } + return u.f; +} + +float floorf(float x) +{ + union {float f; uint32_t i;} u = {x}; + int e = (int)(u.i >> 23 & 0xff) - 0x7f; + uint32_t m; + + if (e >= 23) + return x; + if (e >= 0) { + m = 0x007fffff >> e; + if ((u.i & m) == 0) + return x; + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31) + u.i += m; + u.i &= ~m; + } else { + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31 == 0) + u.i = 0; + else if (u.i << 1) + u.f = -1.0; + } + return u.f; +} + +float truncf(float x) +{ + union {float f; uint32_t i;} u = {x}; + int e = (int)(u.i >> 23 & 0xff) - 0x7f + 9; + uint32_t m; + + if (e >= 23 + 9) + return x; + if (e < 9) + e = 1; + m = -1U >> e; + if ((u.i & m) == 0) + return x; + FORCE_EVAL(x + 0x1p120f); + u.i &= ~m; + return u.f; +} diff --git a/circuitpython/lib/libm/nearbyintf.c b/circuitpython/lib/libm/nearbyintf.c new file mode 100644 index 0000000..5ac3178 --- /dev/null +++ b/circuitpython/lib/libm/nearbyintf.c @@ -0,0 +1,24 @@ +// adapted from the rintf() function from musl-1.1.16 + +#include "libm.h" + +float nearbyintf(float x) +{ + union {float f; uint32_t i;} u = {x}; + int e = u.i>>23 & 0xff; + int s = u.i>>31; + float_t y; + + if (e >= 0x7f+23) + return x; + if (s) + y = x - 0x1p23f + 0x1p23f; + else + y = x + 0x1p23f - 0x1p23f; +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wfloat-equal" + if (y == 0) + return s ? -0.0f : 0.0f; +#pragma GCC diagnostic pop + return y; +} diff --git a/circuitpython/lib/libm/roundf.c b/circuitpython/lib/libm/roundf.c new file mode 100644 index 0000000..3da1f05 --- /dev/null +++ b/circuitpython/lib/libm/roundf.c @@ -0,0 +1,33 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// roundf from musl-0.9.15 +/*****************************************************************************/ +/*****************************************************************************/ + +#include "libm.h" + +float roundf(float x) +{ + union {float f; uint32_t i;} u = {x}; + int e = u.i >> 23 & 0xff; + float_t y; + + if (e >= 0x7f+23) + return x; + if (u.i >> 31) + x = -x; + if (e < 0x7f-1) { + FORCE_EVAL(x + 0x1p23f); + return 0*u.f; + } + y = (float)(x + 0x1p23f) - 0x1p23f - x; + if (y > 0.5f) + y = y + x - 1; + else if (y <= -0.5f) + y = y + x + 1; + else + y = y + x; + if (u.i >> 31) + y = -y; + return y; +} diff --git a/circuitpython/lib/libm/sf_cos.c b/circuitpython/lib/libm/sf_cos.c new file mode 100644 index 0000000..0fe7c8c --- /dev/null +++ b/circuitpython/lib/libm/sf_cos.c @@ -0,0 +1,71 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* sf_cos.c -- float version of s_cos.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float cosf(float x) +#else + float cosf(x) + float x; +#endif +{ + float y[2],z=0.0; + __int32_t n,ix; + + GET_FLOAT_WORD(ix,x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if(ix <= 0x3f490fd8) return __kernel_cosf(x,z); + + /* cos(Inf or NaN) is NaN */ + else if (!FLT_UWORD_IS_FINITE(ix)) return x-x; + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2f(x,y); + switch(n&3) { + case 0: return __kernel_cosf(y[0],y[1]); + case 1: return -__kernel_sinf(y[0],y[1],1); + case 2: return -__kernel_cosf(y[0],y[1]); + default: + return __kernel_sinf(y[0],y[1],1); + } + } +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double cos(double x) +#else + double cos(x) + double x; +#endif +{ + return (double) cosf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/sf_erf.c b/circuitpython/lib/libm/sf_erf.c new file mode 100644 index 0000000..e1b0a78 --- /dev/null +++ b/circuitpython/lib/libm/sf_erf.c @@ -0,0 +1,257 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* sf_erf.c -- float version of s_erf.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#define __ieee754_expf expf + +#ifdef __v810__ +#define const +#endif + +#ifdef __STDC__ +static const float +#else +static float +#endif +tiny = 1e-30f, +half= 5.0000000000e-01f, /* 0x3F000000 */ +one = 1.0000000000e+00f, /* 0x3F800000 */ +two = 2.0000000000e+00f, /* 0x40000000 */ + /* c = (subfloat)0.84506291151 */ +erx = 8.4506291151e-01f, /* 0x3f58560b */ +/* + * Coefficients for approximation to erf on [0,0.84375] + */ +efx = 1.2837916613e-01f, /* 0x3e0375d4 */ +efx8= 1.0270333290e+00f, /* 0x3f8375d4 */ +pp0 = 1.2837916613e-01f, /* 0x3e0375d4 */ +pp1 = -3.2504209876e-01f, /* 0xbea66beb */ +pp2 = -2.8481749818e-02f, /* 0xbce9528f */ +pp3 = -5.7702702470e-03f, /* 0xbbbd1489 */ +pp4 = -2.3763017452e-05f, /* 0xb7c756b1 */ +qq1 = 3.9791721106e-01f, /* 0x3ecbbbce */ +qq2 = 6.5022252500e-02f, /* 0x3d852a63 */ +qq3 = 5.0813062117e-03f, /* 0x3ba68116 */ +qq4 = 1.3249473704e-04f, /* 0x390aee49 */ +qq5 = -3.9602282413e-06f, /* 0xb684e21a */ +/* + * Coefficients for approximation to erf in [0.84375,1.25] + */ +pa0 = -2.3621185683e-03f, /* 0xbb1acdc6 */ +pa1 = 4.1485610604e-01f, /* 0x3ed46805 */ +pa2 = -3.7220788002e-01f, /* 0xbebe9208 */ +pa3 = 3.1834661961e-01f, /* 0x3ea2fe54 */ +pa4 = -1.1089469492e-01f, /* 0xbde31cc2 */ +pa5 = 3.5478305072e-02f, /* 0x3d1151b3 */ +pa6 = -2.1663755178e-03f, /* 0xbb0df9c0 */ +qa1 = 1.0642088205e-01f, /* 0x3dd9f331 */ +qa2 = 5.4039794207e-01f, /* 0x3f0a5785 */ +qa3 = 7.1828655899e-02f, /* 0x3d931ae7 */ +qa4 = 1.2617121637e-01f, /* 0x3e013307 */ +qa5 = 1.3637083583e-02f, /* 0x3c5f6e13 */ +qa6 = 1.1984500103e-02f, /* 0x3c445aa3 */ +/* + * Coefficients for approximation to erfc in [1.25,1/0.35] + */ +ra0 = -9.8649440333e-03f, /* 0xbc21a093 */ +ra1 = -6.9385856390e-01f, /* 0xbf31a0b7 */ +ra2 = -1.0558626175e+01f, /* 0xc128f022 */ +ra3 = -6.2375331879e+01f, /* 0xc2798057 */ +ra4 = -1.6239666748e+02f, /* 0xc322658c */ +ra5 = -1.8460508728e+02f, /* 0xc3389ae7 */ +ra6 = -8.1287437439e+01f, /* 0xc2a2932b */ +ra7 = -9.8143291473e+00f, /* 0xc11d077e */ +sa1 = 1.9651271820e+01f, /* 0x419d35ce */ +sa2 = 1.3765776062e+02f, /* 0x4309a863 */ +sa3 = 4.3456588745e+02f, /* 0x43d9486f */ +sa4 = 6.4538726807e+02f, /* 0x442158c9 */ +sa5 = 4.2900814819e+02f, /* 0x43d6810b */ +sa6 = 1.0863500214e+02f, /* 0x42d9451f */ +sa7 = 6.5702495575e+00f, /* 0x40d23f7c */ +sa8 = -6.0424413532e-02f, /* 0xbd777f97 */ +/* + * Coefficients for approximation to erfc in [1/.35,28] + */ +rb0 = -9.8649431020e-03f, /* 0xbc21a092 */ +rb1 = -7.9928326607e-01f, /* 0xbf4c9dd4 */ +rb2 = -1.7757955551e+01f, /* 0xc18e104b */ +rb3 = -1.6063638306e+02f, /* 0xc320a2ea */ +rb4 = -6.3756646729e+02f, /* 0xc41f6441 */ +rb5 = -1.0250950928e+03f, /* 0xc480230b */ +rb6 = -4.8351919556e+02f, /* 0xc3f1c275 */ +sb1 = 3.0338060379e+01f, /* 0x41f2b459 */ +sb2 = 3.2579251099e+02f, /* 0x43a2e571 */ +sb3 = 1.5367296143e+03f, /* 0x44c01759 */ +sb4 = 3.1998581543e+03f, /* 0x4547fdbb */ +sb5 = 2.5530502930e+03f, /* 0x451f90ce */ +sb6 = 4.7452853394e+02f, /* 0x43ed43a7 */ +sb7 = -2.2440952301e+01f; /* 0xc1b38712 */ + +#ifdef __STDC__ + float erff(float x) +#else + float erff(x) + float x; +#endif +{ + __int32_t hx,ix,i; + float R,S,P,Q,s,y,z,r; + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; + if(!FLT_UWORD_IS_FINITE(ix)) { /* erf(nan)=nan */ + i = ((__uint32_t)hx>>31)<<1; + return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */ + } + + if(ix < 0x3f580000) { /* |x|<0.84375 */ + if(ix < 0x31800000) { /* |x|<2**-28 */ + if (ix < 0x04000000) + /*avoid underflow */ + return (float)0.125*((float)8.0*x+efx8*x); + return x + efx*x; + } + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + return x + x*y; + } + if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */ + s = fabsf(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) return erx + P/Q; else return -erx - P/Q; + } + if (ix >= 0x40c00000) { /* inf>|x|>=6 */ + if(hx>=0) return one-tiny; else return tiny-one; + } + x = fabsf(x); + s = one/(x*x); + if(ix< 0x4036DB6E) { /* |x| < 1/0.35 */ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/0.35 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + GET_FLOAT_WORD(ix,x); + SET_FLOAT_WORD(z,ix&0xfffff000); + r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S); + if(hx>=0) return one-r/x; else return r/x-one; +} + +#ifdef __STDC__ + float erfcf(float x) +#else + float erfcf(x) + float x; +#endif +{ + __int32_t hx,ix; + float R,S,P,Q,s,y,z,r; + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; + if(!FLT_UWORD_IS_FINITE(ix)) { /* erfc(nan)=nan */ + /* erfc(+-inf)=0,2 */ + return (float)(((__uint32_t)hx>>31)<<1)+one/x; + } + + if(ix < 0x3f580000) { /* |x|<0.84375 */ + if(ix < 0x23800000) /* |x|<2**-56 */ + return one-x; + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + if(hx < 0x3e800000) { /* x<1/4 */ + return one-(x+x*y); + } else { + r = x*y; + r += (x-half); + return half - r ; + } + } + if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */ + s = fabsf(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) { + z = one-erx; return z - P/Q; + } else { + z = erx+P/Q; return one+z; + } + } + if (ix < 0x41e00000) { /* |x|<28 */ + x = fabsf(x); + s = one/(x*x); + if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/.35 ~ 2.857143 */ + if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + GET_FLOAT_WORD(ix,x); + SET_FLOAT_WORD(z,ix&0xfffff000); + r = __ieee754_expf(-z*z-(float)0.5625)* + __ieee754_expf((z-x)*(z+x)+R/S); + if(hx>0) return r/x; else return two-r/x; + } else { + if(hx>0) return tiny*tiny; else return two-tiny; + } +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double erf(double x) +#else + double erf(x) + double x; +#endif +{ + return (double) erff((float) x); +} + +#ifdef __STDC__ + double erfc(double x) +#else + double erfc(x) + double x; +#endif +{ + return (double) erfcf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/sf_frexp.c b/circuitpython/lib/libm/sf_frexp.c new file mode 100644 index 0000000..097185c --- /dev/null +++ b/circuitpython/lib/libm/sf_frexp.c @@ -0,0 +1,70 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* sf_frexp.c -- float version of s_frexp.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +two25 = 3.3554432000e+07f; /* 0x4c000000 */ + +#ifdef __STDC__ + float frexpf(float x, int *eptr) +#else + float frexpf(x, eptr) + float x; int *eptr; +#endif +{ + __int32_t hx, ix; + GET_FLOAT_WORD(hx,x); + ix = 0x7fffffff&hx; + *eptr = 0; + if(!FLT_UWORD_IS_FINITE(ix)||FLT_UWORD_IS_ZERO(ix)) return x; /* 0,inf,nan */ + if (FLT_UWORD_IS_SUBNORMAL(ix)) { /* subnormal */ + x *= two25; + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; + *eptr = -25; + } + *eptr += (ix>>23)-126; + hx = (hx&0x807fffff)|0x3f000000; + SET_FLOAT_WORD(x,hx); + return x; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double frexp(double x, int *eptr) +#else + double frexp(x, eptr) + double x; int *eptr; +#endif +{ + return (double) frexpf((float) x, eptr); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/sf_ldexp.c b/circuitpython/lib/libm/sf_ldexp.c new file mode 100644 index 0000000..ebddc1f --- /dev/null +++ b/circuitpython/lib/libm/sf_ldexp.c @@ -0,0 +1,55 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* sf_ldexp.c -- float version of s_ldexp.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float ldexpf(float value, int exp) +#else + float ldexpf(value, exp) + float value; int exp; +#endif +{ +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wfloat-equal" + if(!isfinite(value)||value==(float)0.0) return value; +#pragma GCC diagnostic pop + value = scalbnf(value,exp); + //if(!finitef(value)||value==(float)0.0) errno = ERANGE; + return value; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double ldexp(double value, int exp) +#else + double ldexp(value, exp) + double value; int exp; +#endif +{ + return (double) ldexpf((float) value, exp); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/sf_modf.c b/circuitpython/lib/libm/sf_modf.c new file mode 100644 index 0000000..8141aa8 --- /dev/null +++ b/circuitpython/lib/libm/sf_modf.c @@ -0,0 +1,82 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/common + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* sf_modf.c -- float version of s_modf.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float one = 1.0f; +#else +static float one = 1.0f; +#endif + +#ifdef __STDC__ + float modff(float x, float *iptr) +#else + float modff(x, iptr) + float x,*iptr; +#endif +{ + __int32_t i0,j0; + __uint32_t i; + GET_FLOAT_WORD(i0,x); + j0 = ((i0>>23)&0xff)-0x7f; /* exponent of x */ + if(j0<23) { /* integer part in x */ + if(j0<0) { /* |x|<1 */ + SET_FLOAT_WORD(*iptr,i0&0x80000000); /* *iptr = +-0 */ + return x; + } else { + i = (0x007fffff)>>j0; + if((i0&i)==0) { /* x is integral */ + __uint32_t ix; + *iptr = x; + GET_FLOAT_WORD(ix,x); + SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */ + return x; + } else { + SET_FLOAT_WORD(*iptr,i0&(~i)); + return x - *iptr; + } + } + } else { /* no fraction part */ + __uint32_t ix; + *iptr = x*one; + GET_FLOAT_WORD(ix,x); + SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */ + return x; + } +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double modf(double x, double *iptr) +#else + double modf(x, iptr) + double x,*iptr; +#endif +{ + return (double) modff((float) x, (float *) iptr); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/sf_sin.c b/circuitpython/lib/libm/sf_sin.c new file mode 100644 index 0000000..b9d4049 --- /dev/null +++ b/circuitpython/lib/libm/sf_sin.c @@ -0,0 +1,71 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* sf_sin.c -- float version of s_sin.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float sinf(float x) +#else + float sinf(x) + float x; +#endif +{ + float y[2],z=0.0; + __int32_t n,ix; + + GET_FLOAT_WORD(ix,x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if(ix <= 0x3f490fd8) return __kernel_sinf(x,z,0); + + /* sin(Inf or NaN) is NaN */ + else if (!FLT_UWORD_IS_FINITE(ix)) return x-x; + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2f(x,y); + switch(n&3) { + case 0: return __kernel_sinf(y[0],y[1],1); + case 1: return __kernel_cosf(y[0],y[1]); + case 2: return -__kernel_sinf(y[0],y[1],1); + default: + return -__kernel_cosf(y[0],y[1]); + } + } +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double sin(double x) +#else + double sin(x) + double x; +#endif +{ + return (double) sinf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/sf_tan.c b/circuitpython/lib/libm/sf_tan.c new file mode 100644 index 0000000..42d2df9 --- /dev/null +++ b/circuitpython/lib/libm/sf_tan.c @@ -0,0 +1,66 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* sf_tan.c -- float version of s_tan.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. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float tanf(float x) +#else + float tanf(x) + float x; +#endif +{ + float y[2],z=0.0; + __int32_t n,ix; + + GET_FLOAT_WORD(ix,x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if(ix <= 0x3f490fda) return __kernel_tanf(x,z,1); + + /* tan(Inf or NaN) is NaN */ + else if (!FLT_UWORD_IS_FINITE(ix)) return x-x; /* NaN */ + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2f(x,y); + return __kernel_tanf(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even + -1 -- n odd */ + } +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double tan(double x) +#else + double tan(x) + double x; +#endif +{ + return (double) tanf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/thumb_vfp_sqrtf.c b/circuitpython/lib/libm/thumb_vfp_sqrtf.c new file mode 100644 index 0000000..12ffebf --- /dev/null +++ b/circuitpython/lib/libm/thumb_vfp_sqrtf.c @@ -0,0 +1,11 @@ +// an implementation of sqrtf for Thumb using hardware VFP instructions + +#include <math.h> + +float sqrtf(float x) { + asm volatile ( + "vsqrt.f32 %[r], %[x]\n" + : [r] "=t" (x) + : [x] "t" (x)); + return x; +} diff --git a/circuitpython/lib/libm/wf_lgamma.c b/circuitpython/lib/libm/wf_lgamma.c new file mode 100644 index 0000000..bcf3705 --- /dev/null +++ b/circuitpython/lib/libm/wf_lgamma.c @@ -0,0 +1,96 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* wf_lgamma.c -- float version of w_lgamma.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. + * ==================================================== + * + */ + +#include "fdlibm.h" +#define _IEEE_LIBM 1 + +#ifdef __STDC__ + float lgammaf(float x) +#else + float lgammaf(x) + float x; +#endif +{ +#ifdef _IEEE_LIBM + int sign; + return __ieee754_lgammaf_r(x,&sign); +#else + float y; + struct exception exc; + y = __ieee754_lgammaf_r(x,&(_REENT_SIGNGAM(_REENT))); + if(_LIB_VERSION == _IEEE_) return y; + if(!finitef(y)&&finitef(x)) { +#ifndef HUGE_VAL +#define HUGE_VAL inf + double inf = 0.0; + + SET_HIGH_WORD(inf,0x7ff00000); /* set inf to infinite */ +#endif + exc.name = "lgammaf"; + exc.err = 0; + exc.arg1 = exc.arg2 = (double)x; + if (_LIB_VERSION == _SVID_) + exc.retval = HUGE; + else + exc.retval = HUGE_VAL; + if(floorf(x)==x&&x<=(float)0.0) { + /* lgammaf(-integer) */ + exc.type = SING; + if (_LIB_VERSION == _POSIX_) + errno = EDOM; + else if (!matherr(&exc)) { + errno = EDOM; + } + + } else { + /* lgammaf(finite) overflow */ + exc.type = OVERFLOW; + if (_LIB_VERSION == _POSIX_) + errno = ERANGE; + else if (!matherr(&exc)) { + errno = ERANGE; + } + } + if (exc.err != 0) + errno = exc.err; + return (float)exc.retval; + } else + return y; +#endif +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double lgamma(double x) +#else + double lgamma(x) + double x; +#endif +{ + return (double) lgammaf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/circuitpython/lib/libm/wf_tgamma.c b/circuitpython/lib/libm/wf_tgamma.c new file mode 100644 index 0000000..3ff05f3 --- /dev/null +++ b/circuitpython/lib/libm/wf_tgamma.c @@ -0,0 +1,73 @@ +/* + * This file is part of the MicroPython project, http://micropython.org/ + * + * These math functions are taken from newlib-nano-2, the newlib/libm/math + * directory, available from https://github.com/32bitmicro/newlib-nano-2. + * + * Appropriate copyright headers are reproduced below. + */ + +/* w_gammaf.c -- float version of w_gamma.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. + * ==================================================== + */ + +#include "math.h" +#include "fdlibm.h" +#define _IEEE_LIBM 1 + +#ifdef __STDC__ + float tgammaf(float x) +#else + float tgammaf(x) + float x; +#endif +{ + float y; + int local_signgam; + if (!isfinite(x)) { + /* special cases: tgammaf(nan)=nan, tgammaf(inf)=inf, tgammaf(-inf)=nan */ + return x + INFINITY; + } + y = expf(__ieee754_lgammaf_r(x,&local_signgam)); + if (local_signgam < 0) y = -y; +#ifdef _IEEE_LIBM + return y; +#else + if(_LIB_VERSION == _IEEE_) return y; + + if(!finitef(y)&&finitef(x)) { + if(floorf(x)==x&&x<=(float)0.0) + /* tgammaf pole */ + return (float)__kernel_standard((double)x,(double)x,141); + else + /* tgammaf overflow */ + return (float)__kernel_standard((double)x,(double)x,140); + } + return y; +#endif +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double tgamma(double x) +#else + double tgamma(x) + double x; +#endif +{ + return (double) tgammaf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ |