diff options
Diffstat (limited to 'nuttx/lib/math/lib_exp.c')
| -rw-r--r-- | nuttx/lib/math/lib_exp.c | 253 |
1 files changed, 253 insertions, 0 deletions
diff --git a/nuttx/lib/math/lib_exp.c b/nuttx/lib/math/lib_exp.c new file mode 100644 index 000000000..a70cb10e7 --- /dev/null +++ b/nuttx/lib/math/lib_exp.c @@ -0,0 +1,253 @@ +/* + * Copyright (C) 2009-2011 Nick Johnson <nickbjohnson4224 at gmail.com> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +#include <apps/math.h> +#include <float.h> +#include <stdint.h> +#include <stdbool.h> +#include <unistd.h> + +#define M_E2 (M_E * M_E) +#define M_E4 (M_E2 * M_E2) +#define M_E8 (M_E4 * M_E4) +#define M_E16 (M_E8 * M_E8) +#define M_E32 (M_E16 * M_E16) +#define M_E64 (M_E32 * M_E32) +#define M_E128 (M_E64 * M_E64) +#define M_E256 (M_E128 * M_E128) +#define M_E512 (M_E256 * M_E256) +#define M_E1024 (M_E512 * M_E512) + +static double _expi_square_tbl[11] = { + M_E, // e^1 + M_E2, // e^2 + M_E4, // e^4 + M_E8, // e^8 + M_E16, // e^16 + M_E32, // e^32 + M_E64, // e^64 + M_E128, // e^128 + M_E256, // e^256 + M_E512, // e^512 + M_E1024, // e^1024 +}; + +static double _expi(size_t n) { + size_t i; + double val; + + if (n > 1024) { + return INFINITY; + } + + val = 1.0; + + for (i = 0; n; i++) { + if (n & (1 << i)) { + n &= ~(1 << i); + val *= _expi_square_tbl[i]; + } + } + + return val; +} + +static float _flt_inv_fact[] = { + 1.0 / 1.0, // 1/0! + 1.0 / 1.0, // 1/1! + 1.0 / 2.0, // 1/2! + 1.0 / 6.0, // 1/3! + 1.0 / 24.0, // 1/4! + 1.0 / 120.0, // 1/5! + 1.0 / 720.0, // 1/6! + 1.0 / 5040.0, // 1/7! + 1.0 / 40320.0, // 1/8! + 1.0 / 362880.0, // 1/9! + 1.0 / 3628800.0, // 1/10! +}; + +float expf(float x) { + size_t int_part; + bool invert; + float value; + float x0; + size_t i; + + if (x == 0) { + return 1; + } + else if (x < 0) { + invert = true; + x = -x; + } + else { + invert = false; + } + + /* extract integer component */ + int_part = (size_t) x; + + /* set x to fractional component */ + x -= (float) int_part; + + /* perform Taylor series approximation with eleven terms */ + value = 0.0; + x0 = 1.0; + for (i = 0; i < 10; i++) { + value += x0 * _flt_inv_fact[i]; + x0 *= x; + } + + /* multiply by exp of the integer component */ + value *= _expi(int_part); + + if (invert) { + return (1.0 / value); + } + else { + return value; + } +} + +static double _dbl_inv_fact[] = { + 1.0 / 1.0, // 1 / 0! + 1.0 / 1.0, // 1 / 1! + 1.0 / 2.0, // 1 / 2! + 1.0 / 6.0, // 1 / 3! + 1.0 / 24.0, // 1 / 4! + 1.0 / 120.0, // 1 / 5! + 1.0 / 720.0, // 1 / 6! + 1.0 / 5040.0, // 1 / 7! + 1.0 / 40320.0, // 1 / 8! + 1.0 / 362880.0, // 1 / 9! + 1.0 / 3628800.0, // 1 / 10! + 1.0 / 39916800.0, // 1 / 11! + 1.0 / 479001600.0, // 1 / 12! + 1.0 / 6227020800.0, // 1 / 13! + 1.0 / 87178291200.0, // 1 / 14! + 1.0 / 1307674368000.0, // 1 / 15! + 1.0 / 20922789888000.0, // 1 / 16! + 1.0 / 355687428096000.0, // 1 / 17! + 1.0 / 6402373705728000.0, // 1 / 18! +}; + +double exp(double x) { + size_t int_part; + bool invert; + double value; + double x0; + size_t i; + + if (x == 0) { + return 1; + } + else if (x < 0) { + invert = true; + x = -x; + } + else { + invert = false; + } + + /* extract integer component */ + int_part = (size_t) x; + + /* set x to fractional component */ + x -= (double) int_part; + + /* perform Taylor series approximation with nineteen terms */ + value = 0.0; + x0 = 1.0; + for (i = 0; i < 19; i++) { + value += x0 * _dbl_inv_fact[i]; + x0 *= x; + } + + /* multiply by exp of the integer component */ + value *= _expi(int_part); + + if (invert) { + return (1.0 / value); + } + else { + return value; + } +} + +static long double _ldbl_inv_fact[] = { + 1.0 / 1.0, // 1 / 0! + 1.0 / 1.0, // 1 / 1! + 1.0 / 2.0, // 1 / 2! + 1.0 / 6.0, // 1 / 3! + 1.0 / 24.0, // 1 / 4! + 1.0 / 120.0, // 1 / 5! + 1.0 / 720.0, // 1 / 6! + 1.0 / 5040.0, // 1 / 7! + 1.0 / 40320.0, // 1 / 8! + 1.0 / 362880.0, // 1 / 9! + 1.0 / 3628800.0, // 1 / 10! + 1.0 / 39916800.0, // 1 / 11! + 1.0 / 479001600.0, // 1 / 12! + 1.0 / 6227020800.0, // 1 / 13! + 1.0 / 87178291200.0, // 1 / 14! + 1.0 / 1307674368000.0, // 1 / 15! + 1.0 / 20922789888000.0, // 1 / 16! + 1.0 / 355687428096000.0, // 1 / 17! + 1.0 / 6402373705728000.0, // 1 / 18! +}; + +long double expl(long double x) { + size_t int_part; + bool invert; + long double value; + long double x0; + size_t i; + + if (x == 0) { + return 1; + } + else if (x < 0) { + invert = true; + x = -x; + } + else { + invert = false; + } + + /* extract integer component */ + int_part = (size_t) x; + + /* set x to fractional component */ + x -= (long double) int_part; + + /* perform Taylor series approximation with nineteen terms */ + value = 0.0; + x0 = 1.0; + for (i = 0; i < 19; i++) { + value += x0 * _ldbl_inv_fact[i]; + x0 *= x; + } + + /* multiply by exp of the integer component */ + value *= _expi(int_part); + + if (invert) { + return (1.0 / value); + } + else { + return value; + } +} |