From 20beff9ae3ae50b44022595bcb6d27d195a08963 Mon Sep 17 00:00:00 2001 From: Damien George Date: Thu, 11 Sep 2014 22:24:45 +0100 Subject: [PATCH] py and libm: Add asinf,acosf; print higher precision for float. Also use less stack space when printing single precision float. Addition of asinf and acosf addresses issue #851. --- lib/libm/asinfacosf.c | 130 ++++++++++++++++++++++++++++++++++++++++++ lib/libm/math.c | 2 - py/objcomplex.c | 18 ++++-- py/objfloat.c | 10 ++-- stmhal/Makefile | 1 + 5 files changed, 148 insertions(+), 13 deletions(-) create mode 100644 lib/libm/asinfacosf.c diff --git a/lib/libm/asinfacosf.c b/lib/libm/asinfacosf.c new file mode 100644 index 0000000000..07ecad3f3f --- /dev/null +++ b/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+00, /* 0x3fc90fda */ +pio2_lo = 7.5497894159e-08; /* 0x33a22168 */ + +static const float +/* coefficients for R(x^2) */ +pS0 = 1.6666586697e-01, +pS1 = -4.2743422091e-02, +pS2 = -8.6563630030e-03, +qS1 = -7.0662963390e-01; + +static float R(float z) +{ + float_t p, q; + p = z*(pS0+z*(pS1+z*pS2)); + q = 1.0f+z*qS1; + return p/q; +} + +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/lib/libm/math.c b/lib/libm/math.c index 23218fe56b..07a563e58c 100644 --- a/lib/libm/math.c +++ b/lib/libm/math.c @@ -117,8 +117,6 @@ 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 fmodf(float x, float y) { return 0.0; } float tgammaf(float x) { return 0.0; } float lgammaf(float x) { return 0.0; } diff --git a/py/objcomplex.c b/py/objcomplex.c index 8e0f638c3e..2768a01e65 100644 --- a/py/objcomplex.c +++ b/py/objcomplex.c @@ -25,6 +25,7 @@ */ #include +#include #include #include "mpconfig.h" @@ -55,21 +56,26 @@ mp_obj_t mp_obj_new_complex(mp_float_t real, mp_float_t imag); STATIC void complex_print(void (*print)(void *env, const char *fmt, ...), void *env, mp_obj_t o_in, mp_print_kind_t kind) { mp_obj_complex_t *o = o_in; #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT - char buf[32]; + char buf[16]; if (o->real == 0) { - format_float(o->imag, buf, sizeof(buf), 'g', 6, '\0'); + format_float(o->imag, buf, sizeof(buf), 'g', 7, '\0'); print(env, "%sj", buf); } else { - format_float(o->real, buf, sizeof(buf), 'g', 6, '\0'); + format_float(o->real, buf, sizeof(buf), 'g', 7, '\0'); print(env, "(%s+", buf); - format_float(o->imag, buf, sizeof(buf), 'g', 6, '\0'); + format_float(o->imag, buf, sizeof(buf), 'g', 7, '\0'); print(env, "%sj)", buf); } #else + char buf[32]; if (o->real == 0) { - print(env, "%.8gj", (double) o->imag); + sprintf(buf, "%.16g", (double)o->imag); + print(env, "%sj", buf); } else { - print(env, "(%.8g+%.8gj)", (double) o->real, (double) o->imag); + sprintf(buf, "%.16g", (double)o->real); + print(env, "(%s+", buf); + sprintf(buf, "%.16g", (double)o->imag); + print(env, "%sj)", buf); } #endif } diff --git a/py/objfloat.c b/py/objfloat.c index 52d484135f..c6734ee780 100644 --- a/py/objfloat.c +++ b/py/objfloat.c @@ -48,18 +48,18 @@ STATIC void float_print(void (*print)(void *env, const char *fmt, ...), void *env, mp_obj_t o_in, mp_print_kind_t kind) { mp_obj_float_t *o = o_in; #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT - char buf[32]; - format_float(o->value, buf, sizeof(buf), 'g', 6, '\0'); + char buf[16]; + format_float(o->value, buf, sizeof(buf), 'g', 7, '\0'); print(env, "%s", buf); - if (strchr(buf, '.') == NULL) { + if (strchr(buf, '.') == NULL && strchr(buf, 'e') == NULL) { // Python floats always have decimal point print(env, ".0"); } #else char buf[32]; - sprintf(buf, "%.17g", (double) o->value); + sprintf(buf, "%.16g", (double) o->value); print(env, buf); - if (strchr(buf, '.') == NULL) { + if (strchr(buf, '.') == NULL && strchr(buf, 'e') == NULL) { // Python floats always have decimal point print(env, ".0"); } diff --git a/stmhal/Makefile b/stmhal/Makefile index 3ae1805b17..82d229b873 100644 --- a/stmhal/Makefile +++ b/stmhal/Makefile @@ -64,6 +64,7 @@ endif SRC_LIB = $(addprefix lib/,\ libm/math.c \ libm/mathsincos.c \ + libm/asinfacosf.c \ libm/atanf.c \ libm/atan2f.c \ )