From 8707ea34218a0aae071d4fb8d25e23273c828ea1 Mon Sep 17 00:00:00 2001 From: Damien George Date: Fri, 29 Aug 2014 22:42:26 +0100 Subject: [PATCH] lib: Add lib and libm, moving current files from stmhal. Top-level lib directory is for standard C libraries that we want to provide our own versions of (for efficiency and stand-alone reasons). It currently has libm in it for math functions. Also add atanf and atan2f, which addresses issue #837. --- lib/libm/atan2f.c | 89 ++++++++++++++++++++++++++ lib/libm/atanf.c | 100 ++++++++++++++++++++++++++++++ lib/libm/libm.h | 52 ++++++++++++++++ {stmhal => lib/libm}/math.c | 56 +---------------- {stmhal => lib/libm}/mathsincos.c | 0 stmhal/Makefile | 10 ++- 6 files changed, 250 insertions(+), 57 deletions(-) create mode 100644 lib/libm/atan2f.c create mode 100644 lib/libm/atanf.c create mode 100644 lib/libm/libm.h rename {stmhal => lib/libm}/math.c (91%) rename {stmhal => lib/libm}/mathsincos.c (100%) diff --git a/lib/libm/atan2f.c b/lib/libm/atan2f.c new file mode 100644 index 0000000000..03d000c9e1 --- /dev/null +++ b/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+00, /* 0x40490fdb */ +pi_lo = -8.7422776573e-08; /* 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/lib/libm/atanf.c b/lib/libm/atanf.c new file mode 100644 index 0000000000..053fc1b6d6 --- /dev/null +++ b/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-01, /* atan(0.5)hi 0x3eed6338 */ + 7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */ + 9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */ + 1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */ +}; + +static const float atanlo[] = { + 5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */ + 3.7748947079e-08, /* atan(1.0)lo 0x33222168 */ + 3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */ + 7.5497894159e-08, /* atan(inf)lo 0x33a22168 */ +}; + +static const float aT[] = { + 3.3333328366e-01, + -1.9999158382e-01, + 1.4253635705e-01, + -1.0648017377e-01, + 6.1687607318e-02, +}; + +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/lib/libm/libm.h b/lib/libm/libm.h new file mode 100644 index 0000000000..c87c558219 --- /dev/null +++ b/lib/libm/libm.h @@ -0,0 +1,52 @@ +/*****************************************************************************/ +/*****************************************************************************/ +// 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 +#include + +#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/stmhal/math.c b/lib/libm/math.c similarity index 91% rename from stmhal/math.c rename to lib/libm/math.c index 9f0d312488..8433081b97 100644 --- a/stmhal/math.c +++ b/lib/libm/math.c @@ -24,8 +24,7 @@ * THE SOFTWARE. */ -#include -#include +#include "libm.h" typedef float float_t; typedef union { @@ -117,11 +116,8 @@ float tanhf(float x) { return sinhf(x) / coshf(x); } float acoshf(float x) { return 0.0; } float asinhf(float x) { return 0.0; } float atanhf(float x) { return 0.0; } -float tanf(float x) { return 0.0; } float acosf(float x) { return 0.0; } float asinf(float x) { return 0.0; } -float atanf(float x) { return 0.0; } -float atan2f(float x, float y) { return 0.0; } float fmodf(float x, float y) { return 0.0; } float tgammaf(float x) { return 0.0; } float lgammaf(float x) { return 0.0; } @@ -131,56 +127,6 @@ float modff(float x, float *y) { return 0.0; } float frexpf(float x, int *exp) { return 0.0; } float ldexpf(float x, int exp) { return 0.0; } -/*****************************************************************************/ -/*****************************************************************************/ -// 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. - * ==================================================== - */ - -#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) - /*****************************************************************************/ /*****************************************************************************/ // __fpclassifyf from musl-0.9.15 diff --git a/stmhal/mathsincos.c b/lib/libm/mathsincos.c similarity index 100% rename from stmhal/mathsincos.c rename to lib/libm/mathsincos.c diff --git a/stmhal/Makefile b/stmhal/Makefile index ba965481a3..707b898a52 100644 --- a/stmhal/Makefile +++ b/stmhal/Makefile @@ -60,6 +60,13 @@ endif # uncomment this if you want libgcc #LIBS += $(shell $(CC) -print-libgcc-file-name) +SRC_LIB = $(addprefix lib/,\ + libm/math.c \ + libm/mathsincos.c \ + libm/atanf.c \ + libm/atan2f.c \ + ) + SRC_C = \ main.c \ string0.c \ @@ -84,8 +91,6 @@ SRC_C = \ uart.c \ usb.c \ printf.c \ - math.c \ - mathsincos.c \ gccollect.c \ pybstdio.c \ readline.c \ @@ -193,6 +198,7 @@ SRC_CC3K = $(addprefix $(CC3K_DIR)/,\ OBJ = OBJ += $(PY_O) +OBJ += $(addprefix $(BUILD)/, $(SRC_LIB:.c=.o)) OBJ += $(addprefix $(BUILD)/, $(SRC_C:.c=.o)) OBJ += $(addprefix $(BUILD)/, $(SRC_S:.s=.o)) OBJ += $(addprefix $(BUILD)/, $(SRC_HAL:.c=.o))