micropython: add micropython component

components/language/micropython/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;
+}

components/language/micropython/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);
+}

components/language/micropython/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;
+}

components/language/micropython/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(-,-) */
+	}
+}

components/language/micropython/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;
+}

components/language/micropython/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;
+}

components/language/micropython/lib/libm/ef_rem_pio2.c

@@ -0,0 +1,202 @@
+/*
+ * 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(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&&(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;
+	}
+    /* 
+     * 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;
+	while(tx[nx-1]==zero) nx--;	/* skip zero term */
+	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;
+}

components/language/micropython/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;
+}

components/language/micropython/lib/libm/erf_lgamma.c

@@ -0,0 +1,255 @@
+/*
+ * 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"
+
+#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;
+}

components/language/micropython/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)

components/language/micropython/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;
+}

components/language/micropython/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));
+	}
+}

components/language/micropython/lib/libm/kf_rem_pio2.c

@@ -0,0 +1,218 @@
+/*
+ * 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 */
+	if(z==zero) {
+	    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 */
+	if(z==(float)0.0) {
+	    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;
+}

components/language/micropython/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);
+}

components/language/micropython/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);
+	}
+}

components/language/micropython/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)

components/language/micropython/lib/libm/log1pf.c

@@ -0,0 +1,83 @@
+/*****************************************************************************/
+/*****************************************************************************/
+// 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"
+
+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;
+}

components/language/micropython/lib/libm/math.c

@@ -0,0 +1,826 @@
+/*
+ * This file is part of the MicroPython project, http://micropython.org/
+ *
+ * The MIT License (MIT)
+ *
+ * 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; }
+
+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 == 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
+/*****************************************************************************/
+/*****************************************************************************/
+
+/* 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;
+
+/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
+float __expo2f(float x)
+{
+	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;
+}

components/language/micropython/lib/libm/nearbyintf.c

@@ -0,0 +1,21 @@
+// 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;
+	if (y == 0)
+		return s ? -0.0f : 0.0f;
+	return y;
+}

components/language/micropython/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;
+}

components/language/micropython/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) */

components/language/micropython/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) */

components/language/micropython/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) */

components/language/micropython/lib/libm/sf_ldexp.c

@@ -0,0 +1,52 @@
+/*
+ * 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
+{
+	if(!isfinite(value)||value==(float)0.0) return value;
+	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) */

components/language/micropython/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) */

components/language/micropython/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) */

components/language/micropython/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) */

components/language/micropython/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;
+}

components/language/micropython/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) */

components/language/micropython/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) */