diff --git a/libc/math/lib_asin.c b/libc/math/lib_asin.c index 2f45a24a1d9..15bce7913a1 100644 --- a/libc/math/lib_asin.c +++ b/libc/math/lib_asin.c @@ -35,6 +35,8 @@ #include #include +#ifdef CONFIG_HAVE_DOUBLE + /**************************************************************************** * Pre-processor Definitions ****************************************************************************/ @@ -42,16 +44,40 @@ #undef DBL_EPSILON #define DBL_EPSILON 1e-12 +/**************************************************************************** + * Private Functions + ****************************************************************************/ + +/* This lib uses Newton's method to approximate asin(x). Newton's Method + * converges very slowly for x close to 1. We can accelerate convergence + * with the following identy: asin(x)=Sign(x)*(Pi/2-asin(sqrt(1-x^2))) + */ + +static double asin_aux(double x) +{ + long double y; + double y_cos, y_sin; + + y = 0.0; + y_sin = 0.0; + + while (fabs(y_sin - x) > DBL_EPSILON) + { + y_cos = cos(y); + y -= ((long double)y_sin - (long double)x) / (long double)y_cos; + y_sin = sin(y); + } + + return y; +} + /**************************************************************************** * Public Functions ****************************************************************************/ -#ifdef CONFIG_HAVE_DOUBLE double asin(double x) { - long double y; - long double y_sin; - long double y_cos; + double y; /* Verify that the input value is in the domain of the function */ @@ -60,26 +86,19 @@ double asin(double x) return NAN; } - y = 0; + /* if x is > sqrt(2), use identity for faster convergence */ - while (1) + if (fabs(x) > 0.71) { - y_sin = sin(y); - y_cos = cos(y); - - if (y > M_PI_2 || y < -M_PI_2) - { - y = fmod(y, M_PI); - } - - if (y_sin + DBL_EPSILON >= x && y_sin - DBL_EPSILON <= x) - { - break; - } - - y = y - (y_sin - x) / y_cos; + y = M_PI_2 - asin_aux(sqrt(1.0 - x * x)); + y = copysign(y, x); + } + else + { + y = asin_aux(x); } return y; } -#endif + +#endif /* CONFIG_HAVE_DOUBLE */ diff --git a/libc/math/lib_asinf.c b/libc/math/lib_asinf.c index ac17a539441..e60f439f6e7 100644 --- a/libc/math/lib_asinf.c +++ b/libc/math/lib_asinf.c @@ -3,7 +3,7 @@ * * This file is a part of NuttX: * - * Copyright (C) 2012 Gregory Nutt. All rights reserved. + * Copyright (C) 2012, 2016 Gregory Nutt. All rights reserved. * Ported by: Darcy Gong * * It derives from the Rhombs OS math library by Nick Johnson which has @@ -33,33 +33,58 @@ #include /**************************************************************************** - * Public Functions + * Private Functions ****************************************************************************/ -float asinf(float x) +/* This lib uses Newton's method to approximate asin(x). Newton's Method + * converges very slowly for x close to 1. We can accelerate convergence + * with the following identy: asin(x)=Sign(x)*(Pi/2-asin(sqrt(1-x^2))) + */ + +static float asinf_aux(float x) { - long double y, y_sin, y_cos; + double y; + float y_sin, y_cos; - y = 0; + y = 0.0; + y_sin = 0.0F; - while (1) + while (fabsf(y_sin - x) > FLT_EPSILON) { - y_sin = sinf(y); y_cos = cosf(y); - - if (y > M_PI_2_F || y < -M_PI_2_F) - { - y = fmodf(y, M_PI_F); - } - - if (y_sin + FLT_EPSILON >= x && y_sin - FLT_EPSILON <= x) - { - break; - } - - y = y - (y_sin - x) / y_cos; + y -= ((double)y_sin - (double)x) / (double)y_cos; + y_sin = sinf(y); } return y; } +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +float asinf(float x) +{ + float y; + + /* Verify that the input value is in the domain of the function */ + + if (x < -1.0F || x > 1.0F || isnan(x)) + { + return NAN_F; + } + + /* if x is > sqrt(2), use identity for faster convergence */ + + if (fabsf(x) > 0.71F) + { + y = M_PI_2_F - asinf_aux(sqrtf(1.0F - x * x)); + y = copysignf(y, x); + } + else + { + y = asinf_aux(x); + } + + return y; +} diff --git a/libc/math/lib_asinl.c b/libc/math/lib_asinl.c index 466910daf31..a8cd862d950 100644 --- a/libc/math/lib_asinl.c +++ b/libc/math/lib_asinl.c @@ -3,7 +3,7 @@ * * This file is a part of NuttX: * - * Copyright (C) 2012 Gregory Nutt. All rights reserved. + * Copyright (C) 2012, 2016 Gregory Nutt. All rights reserved. * Ported by: Darcy Gong * * It derives from the Rhombs OS math library by Nick Johnson which has @@ -35,35 +35,57 @@ #include #include +#ifdef CONFIG_HAVE_LONG_DOUBLE + /**************************************************************************** * Public Functions ****************************************************************************/ -#ifdef CONFIG_HAVE_LONG_DOUBLE -long double asinl(long double x) +static long double asinl_aux(long double x) { - long double y, y_sin, y_cos; + long double y, y_cos, y_sin; - y = 0; + y = 0.0; + y_sin = 0.0; - while (1) + while (fabsl(y_sin - x) > DBL_EPSILON) { - y_sin = sinl(y); y_cos = cosl(y); - - if (y > M_PI_2 || y < -M_PI_2) - { - y = fmodl(y, M_PI); - } - - if (y_sin + LDBL_EPSILON >= x && y_sin - LDBL_EPSILON <= x) - { - break; - } - - y = y - (y_sin - x) / y_cos; + y -= (y_sin - x) / y_cos; + y_sin = sinl(y); } return y; } -#endif + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +long double asinl(long double x) +{ + long double y; + + /* Verify that the input value is in the domain of the function */ + + if (x < -1.0 || x > 1.0 || isnan(x)) + { + return NAN; + } + + /* if x is > sqrt(2), use identity for faster convergence */ + + if (fabsl(x) > 0.71) + { + y = M_PI_2 - asinl_aux(sqrtl(1.0 - x * x)); + y = copysignl(y, x); + } + else + { + y = asinl_aux(x); + } + + return y; +} + +endif /* CONFIG_HAVE_LONG_DOUBLE */ diff --git a/libc/math/lib_copysignf.c b/libc/math/lib_copysignf.c index 9684f68a746..c438e7d43ae 100644 --- a/libc/math/lib_copysignf.c +++ b/libc/math/lib_copysignf.c @@ -1,10 +1,20 @@ /**************************************************************************** * libc/math/lib_copysignf.c * - * Copyright (C) 2015 Gregory Nutt. All rights reserved. + * Copyright (C) 2015, 2016 Gregory Nutt. All rights reserved. * Authors: Gregory Nutt * Dave Marples * + * Replaced on 2016-07-30 by David Alession with a faster version of + * copysignf() from NetBSD with the following Copyright: + * + * 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. + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: @@ -42,6 +52,44 @@ #include #include +#include + +/**************************************************************************** + * Pre-processor Definitions + ****************************************************************************/ + +/* 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) + +/**************************************************************************** + * Private Types + ****************************************************************************/ + +/* 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; /**************************************************************************** * Public Functions @@ -49,10 +97,12 @@ float copysignf(float x, float y) { - if (y < 0) - { - return -fabsf(x); - } + uint32_t ix; + uint32_t iy; - return fabsf(x); + GET_FLOAT_WORD(ix, x); + GET_FLOAT_WORD(iy, y); + SET_FLOAT_WORD(x, (ix & 0x7fffffff) | (iy & 0x80000000)); + + return x; } diff --git a/libc/math/lib_erf.c b/libc/math/lib_erf.c index 215c1017204..28f266196f0 100644 --- a/libc/math/lib_erf.c +++ b/libc/math/lib_erf.c @@ -1,6 +1,7 @@ /**************************************************************************** * libc/math/lib_erf.c * + * Copyright (C) 2016 Gregory Nutt. All rights reserved. * Copyright (C) 2015 Brennan Ashton. All rights reserved. * Author: Brennan Ashton * @@ -42,32 +43,42 @@ #include +#ifdef CONFIG_HAVE_DOUBLE + +/**************************************************************************** + * Pre-processor Definitions + ****************************************************************************/ + +#define A1 0.254829592 +#define A2 (-0.284496736) +#define A3 1.421413741 +#define A4 (-1.453152027) +#define A5 1.061405429 +#define P 0.3275911 + /**************************************************************************** * Public Functions ****************************************************************************/ -#ifdef CONFIG_HAVE_DOUBLE +/**************************************************************************** + * Name: erf + * + * Description: + * This implementation comes from the Handbook of Mathmatical Functions + * The implementations in this book are not protected by copyright. + * erf comes from formula 7.1.26 + * + ****************************************************************************/ + double erf(double x) { - /* This implementation comes from the Handbook of Mathmatical Functions - * The implementations in this book are not protected by copyright. - * erf comes from formula 7.1.26 - */ - - char sign; double t; - double a1, a2, a3, a4, a5, p; + double z; - a1 = 0.254829592; - a2 = -0.284496736; - a3 = 1.421413741; - a4 = -1.453152027; - a5 = 1.061405429; - p = 0.3275911; - - sign = (x >= 0 ? 1 : -1); - t = 1.0/(1.0 + p*x); - return sign * (1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * - (double)expf(-x * x)); + z = fabs(x); + t = 1.0 / (1.0 + P * z); + t = 1.0 - (((((A5 * t + A4) * t) + A3) * t + A2) * t + A1) * t * exp(-z * z); + return copysign(t, x); } -#endif + +#endif /* CONFIG_HAVE_DOUBLE */ diff --git a/libc/math/lib_erff.c b/libc/math/lib_erff.c index 364b7fe03c0..8ca0fefb0fe 100644 --- a/libc/math/lib_erff.c +++ b/libc/math/lib_erff.c @@ -1,6 +1,7 @@ /**************************************************************************** * libc/math/lib_erff.c * + * Copyright (C) 2016 Gregory Nutt. All rights reserved. * Copyright (C) 2015 Brennan Ashton. All rights reserved. * Author: Brennan Ashton * @@ -42,29 +43,38 @@ #include +/**************************************************************************** + * Pre-processor Definitions + ****************************************************************************/ + +#define A1 0.254829592F +#define A2 (-0.284496736F) +#define A3 1.421413741F +#define A4 (-1.453152027F) +#define A5 1.061405429F +#define P 0.3275911F + /**************************************************************************** * Public Functions ****************************************************************************/ +/**************************************************************************** + * Name: erff + * + * Description: + * This implementation comes from the Handbook of Mathmatical Functions + * The implementations in this book are not protected by copyright. + * erf comes from formula 7.1.26 + * + ****************************************************************************/ + float erff(float x) { - /* This implementation comes from the Handbook of Mathmatical Functions - * The implementations in this book are not protected by copyright. - * erf comes from formula 7.1.26 - */ - - char sign; float t; - float a1, a2, a3, a4, a5, p; + float z; - a1 = 0.254829592F; - a2 = -0.284496736F; - a3 = 1.421413741F; - a4 = -1.453152027F; - a5 = 1.061405429F; - p = 0.3275911F; - - sign = (x >= 0 ? 1 : -1); - t = 1.0F/(1.0F + p*x); - return sign * (1.0F - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * expf(-x * x)); + z = fabsf(x); + t = 1.0F / (1.0F + P * z); + t = 1.0F - (((((A5 * t + A4) * t) + A3) * t + A2) * t + A1) * t * expf(-z * z); + return copysignf(t, x); } diff --git a/libc/math/lib_erfl.c b/libc/math/lib_erfl.c index 44d62389808..2e674bb1497 100644 --- a/libc/math/lib_erfl.c +++ b/libc/math/lib_erfl.c @@ -46,6 +46,13 @@ * Public Functions ****************************************************************************/ +#define A1 0.254829592 +#define A2 (-0.284496736) +#define A3 1.421413741 +#define A4 (-1.453152027) +#define A5 1.061405429 +#define P 0.3275911 + #ifdef CONFIG_HAVE_LONG_DOUBLE long double erfl(long double x) { @@ -54,19 +61,11 @@ long double erfl(long double x) * erf comes from formula 7.1.26 */ - char sign; - long double t; - long double a1, a2, a3, a4, a5, p; + long double t, z; - a1 = 0.254829592; - a2 = -0.284496736; - a3 = 1.421413741; - a4 = -1.453152027; - a5 = 1.061405429; - p = 0.3275911; - - sign = (x >= 0 ? 1 : -1); - t = 1.0/(1.0 + p*x); - return sign * (1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * expf(-x * x)); + z = fabsl(x); + t = 1.0 / (1.0 + P * z); + t = 1.0 - (((((A5 * t + A4) * t) + A3) * t + A2) * t + A1) * t * expl(-z * z); + return copysignl(t, x); } #endif