Skip to content

Commit 20beff9

Browse files
dpgeorgedpgeorge
committed
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 adafruit#851.
1 parent 9530743 commit 20beff9

5 files changed

Lines changed: 148 additions & 13 deletions

File tree

‎lib/libm/asinfacosf.c‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
/*****************************************************************************/
2+
/*****************************************************************************/
3+
// asinf from musl-0.9.15
4+
/*****************************************************************************/
5+
/*****************************************************************************/
6+
7+
/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.c */
8+
/*
9+
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
10+
*/
11+
/*
12+
* ====================================================
13+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
14+
*
15+
* Developed at SunPro, a Sun Microsystems, Inc. business.
16+
* Permission to use, copy, modify, and distribute this
17+
* software is freely granted, provided that this notice
18+
* is preserved.
19+
* ====================================================
20+
*/
21+
22+
#include "libm.h"
23+
24+
// dpgeorge: pio2 was double in original implementation of asinf
25+
static const float
26+
pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
27+
pio2_lo = 7.5497894159e-08; /* 0x33a22168 */
28+
29+
static const float
30+
/* coefficients for R(x^2) */
31+
pS0 = 1.6666586697e-01,
32+
pS1 = -4.2743422091e-02,
33+
pS2 = -8.6563630030e-03,
34+
qS1 = -7.0662963390e-01;
35+
36+
static float R(float z)
37+
{
38+
float_t p, q;
39+
p = z*(pS0+z*(pS1+z*pS2));
40+
q = 1.0f+z*qS1;
41+
return p/q;
42+
}
43+
44+
float asinf(float x)
45+
{
46+
// dpgeorge: s was double in original implementation
47+
float s,z;
48+
uint32_t hx,ix;
49+
50+
GET_FLOAT_WORD(hx, x);
51+
ix = hx & 0x7fffffff;
52+
if (ix >= 0x3f800000) { /* |x| >= 1 */
53+
if (ix == 0x3f800000) /* |x| == 1 */
54+
return x*pio2_hi + 0x1p-120f; /* asin(+-1) = +-pi/2 with inexact */
55+
return 0/(x-x); /* asin(|x|>1) is NaN */
56+
}
57+
if (ix < 0x3f000000) { /* |x| < 0.5 */
58+
/* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */
59+
if (ix < 0x39800000 && ix >= 0x00800000)
60+
return x;
61+
return x + x*R(x*x);
62+
}
63+
/* 1 > |x| >= 0.5 */
64+
z = (1 - fabsf(x))*0.5f;
65+
s = sqrtf(z);
66+
x = pio2_hi - (2*(s+s*R(z)) - pio2_lo); // dpgeorge: use pio2_hi and pio2_lo
67+
if (hx >> 31)
68+
return -x;
69+
return x;
70+
}
71+
72+
/*****************************************************************************/
73+
/*****************************************************************************/
74+
// acosf from musl-0.9.15
75+
/*****************************************************************************/
76+
/*****************************************************************************/
77+
78+
/* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.c */
79+
/*
80+
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
81+
*/
82+
/*
83+
* ====================================================
84+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
85+
*
86+
* Developed at SunPro, a Sun Microsystems, Inc. business.
87+
* Permission to use, copy, modify, and distribute this
88+
* software is freely granted, provided that this notice
89+
* is preserved.
90+
* ====================================================
91+
*/
92+
93+
float acosf(float x)
94+
{
95+
float z,w,s,c,df;
96+
uint32_t hx,ix;
97+
98+
GET_FLOAT_WORD(hx, x);
99+
ix = hx & 0x7fffffff;
100+
/* |x| >= 1 or nan */
101+
if (ix >= 0x3f800000) {
102+
if (ix == 0x3f800000) {
103+
if (hx >> 31)
104+
return 2*pio2_hi + 0x1p-120f;
105+
return 0;
106+
}
107+
return 0/(x-x);
108+
}
109+
/* |x| < 0.5 */
110+
if (ix < 0x3f000000) {
111+
if (ix <= 0x32800000) /* |x| < 2**-26 */
112+
return pio2_hi + 0x1p-120f;
113+
return pio2_hi - (x - (pio2_lo-x*R(x*x)));
114+
}
115+
/* x < -0.5 */
116+
if (hx >> 31) {
117+
z = (1+x)*0.5f;
118+
s = sqrtf(z);
119+
w = R(z)*s-pio2_lo;
120+
return 2*(pio2_hi - (s+w));
121+
}
122+
/* x > 0.5 */
123+
z = (1-x)*0.5f;
124+
s = sqrtf(z);
125+
GET_FLOAT_WORD(hx,s);
126+
SET_FLOAT_WORD(df,hx&0xfffff000);
127+
c = (z-df*df)/(s+df);
128+
w = R(z)*s+c;
129+
return 2*(df+w);
130+
}

‎lib/libm/math.c‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -117,8 +117,6 @@ float acoshf(float x) { return 0.0; }
117117
float asinhf(float x) { return 0.0; }
118118
float atanhf(float x) { return 0.0; }
119119
float tanf(float x) { return 0.0; }
120-
float acosf(float x) { return 0.0; }
121-
float asinf(float x) { return 0.0; }
122120
float fmodf(float x, float y) { return 0.0; }
123121
float tgammaf(float x) { return 0.0; }
124122
float lgammaf(float x) { return 0.0; }

‎py/objcomplex.c‎

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
*/
2626

2727
#include <stdlib.h>
28+
#include <stdio.h>
2829
#include <assert.h>
2930

3031
#include "mpconfig.h"
@@ -55,21 +56,26 @@ mp_obj_t mp_obj_new_complex(mp_float_t real, mp_float_t imag);
5556
STATIC void complex_print(void (*print)(void *env, const char *fmt, ...), void *env, mp_obj_t o_in, mp_print_kind_t kind) {
5657
mp_obj_complex_t *o = o_in;
5758
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
58-
char buf[32];
59+
char buf[16];
5960
if (o->real == 0) {
60-
format_float(o->imag, buf, sizeof(buf), 'g', 6, '\0');
61+
format_float(o->imag, buf, sizeof(buf), 'g', 7, '\0');
6162
print(env, "%sj", buf);
6263
} else {
63-
format_float(o->real, buf, sizeof(buf), 'g', 6, '\0');
64+
format_float(o->real, buf, sizeof(buf), 'g', 7, '\0');
6465
print(env, "(%s+", buf);
65-
format_float(o->imag, buf, sizeof(buf), 'g', 6, '\0');
66+
format_float(o->imag, buf, sizeof(buf), 'g', 7, '\0');
6667
print(env, "%sj)", buf);
6768
}
6869
#else
70+
char buf[32];
6971
if (o->real == 0) {
70-
print(env, "%.8gj", (double) o->imag);
72+
sprintf(buf, "%.16g", (double)o->imag);
73+
print(env, "%sj", buf);
7174
} else {
72-
print(env, "(%.8g+%.8gj)", (double) o->real, (double) o->imag);
75+
sprintf(buf, "%.16g", (double)o->real);
76+
print(env, "(%s+", buf);
77+
sprintf(buf, "%.16g", (double)o->imag);
78+
print(env, "%sj)", buf);
7379
}
7480
#endif
7581
}

‎py/objfloat.c‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -48,18 +48,18 @@
4848
STATIC void float_print(void (*print)(void *env, const char *fmt, ...), void *env, mp_obj_t o_in, mp_print_kind_t kind) {
4949
mp_obj_float_t *o = o_in;
5050
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
51-
char buf[32];
52-
format_float(o->value, buf, sizeof(buf), 'g', 6, '\0');
51+
char buf[16];
52+
format_float(o->value, buf, sizeof(buf), 'g', 7, '\0');
5353
print(env, "%s", buf);
54-
if (strchr(buf, '.') == NULL) {
54+
if (strchr(buf, '.') == NULL && strchr(buf, 'e') == NULL) {
5555
// Python floats always have decimal point
5656
print(env, ".0");
5757
}
5858
#else
5959
char buf[32];
60-
sprintf(buf, "%.17g", (double) o->value);
60+
sprintf(buf, "%.16g", (double) o->value);
6161
print(env, buf);
62-
if (strchr(buf, '.') == NULL) {
62+
if (strchr(buf, '.') == NULL && strchr(buf, 'e') == NULL) {
6363
// Python floats always have decimal point
6464
print(env, ".0");
6565
}

‎stmhal/Makefile‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,7 @@ endif
6464
SRC_LIB = $(addprefix lib/,\
6565
libm/math.c \
6666
libm/mathsincos.c \
67+
libm/asinfacosf.c \
6768
libm/atanf.c \
6869
libm/atan2f.c \
6970
)

0 commit comments

Comments
 (0)