Part I of port of Rhombus math library

git-svn-id: svn://svn.code.sf.net/p/nuttx/code/trunk@5269 42af7a65-404d-4744-a932-0658087f49c3
This commit is contained in:
patacongo
2012-10-28 14:31:57 +00:00
parent 1d0fad5f70
commit 918e70fe8e
233 changed files with 1649 additions and 214 deletions
+19 -4
View File
@@ -30,12 +30,27 @@ config LIB_HOMEDIR
---help---
The home directory to use with operations like such as 'cd ~'
config HAVE_LIBM
bool "Architecture-specific libm.a"
config LIBM
bool "Math library"
default n
depends on !ARCH_MATH_H
---help---
Architecture specific logic provides an implementation of libm.a
and a math.h header file that can be found at include/arch/math.h.
By default, no math library will be provided by NuttX. In this this case, it
is assumed that (1) no math library is required, or (2) you will be using the
math.h header file and the libm library provided by your toolchain.
This is may be a very good choice is possible because your toolchain may have
have a highly optimized version of libm.
Another possibility is that you have a custom, architecture-specific math
libary and that the corresponding math.h file resides at arch/<architecture>/include/math.h.
The option is selected via ARCH_MATH_H. If ARCH_MATH_H is selected,then the include/nuttx/math.h
header file will be copied to include/math.h where it can be used by your applications.
If ARCH_MATH_H is not defined, then this option can be selected to build a generic,
math library built into NuttX. This math library comes from the Rhombus OS and
was written by Nick Johnson. The Rhombus OS math library port was contributed by
Darcy Gong.
config NOPRINTF_FIELDWIDTH
bool "Disable sprintf support fieldwidth"
+1 -1
View File
@@ -50,7 +50,7 @@ include pthread/Make.defs
include semaphore/Make.defs
include signal/Make.defs
include mqueue/Make.defs
include math/Make.defs
include fixedmath/Make.defs
include net/Make.defs
include time/Make.defs
include libgen/Make.defs
+2 -1
View File
@@ -26,7 +26,8 @@ in the include/ directory provides the prototype for library functions. So
we have:
libgen - libgen.h
math - math.h and fixedmath.h
fixedmath - fixedmath.h
math - math.h
mqueue - pthread.h
net - Various network-related header files: netinet/ether.h, arpa/inet.h
pthread - pthread.h
+43
View File
@@ -0,0 +1,43 @@
############################################################################
# lib/fixedmath/Make.defs
#
# Copyright (C) 2011-2012 Gregory Nutt. All rights reserved.
# Author: Gregory Nutt <gnutt@nuttx.org>
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
# 3. Neither the name NuttX nor the names of its contributors may be
# used to endorse or promote products derived from this software
# without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
# AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#
############################################################################
# Add the fixed precision math C files to the build
CSRCS += lib_rint.c lib_fixedmath.c lib_b16sin.c lib_b16cos.c lib_b16atan2.c
# Add the fixed precision math directory to the build
DEPPATH += --dep-path fixedmath
VPATH += :fixedmath
@@ -1,5 +1,5 @@
/****************************************************************************
* lib/math/lib_b16atan2.c
* lib/fixedmath/lib_b16atan2.c
*
* Copyright (C) 2011 Gregory Nutt. All rights reserved.
* Author: Gregory Nutt <gnutt@nuttx.org>
@@ -1,5 +1,5 @@
/****************************************************************************
* lib/math/lib_b16cos.c
* lib/fixedmath/lib_b16cos.c
*
* Copyright (C) 2007, 2008, 2011 Gregory Nutt. All rights reserved.
* Author: Gregory Nutt <gnutt@nuttx.org>
@@ -1,5 +1,5 @@
/****************************************************************************
* lib/math/lib_b16sin.c
* lib/fixedmath/lib_b16sin.c
*
* Copyright (C) 2008, 2011 Gregory Nutt. All rights reserved.
* Author: Gregory Nutt <gnutt@nuttx.org>
@@ -1,5 +1,5 @@
/************************************************************
* lib/math/lib_rint.c
* lib/fixedmath/lib_rint.c
*
* Copyright (C) 2007, 2011 Gregory Nutt. All rights reserved.
* Author: Gregory Nutt <gnutt@nuttx.org>
+4
View File
@@ -0,0 +1,4 @@
#
# For a description of the syntax of this configuration file,
# see misc/tools/kconfig-language.txt.
#
+10 -4
View File
@@ -1,7 +1,7 @@
############################################################################
# lib/math/Make.defs
#
# Copyright (C) 2011-2012 Gregory Nutt. All rights reserved.
# Copyright (C) 2012 Gregory Nutt. All rights reserved.
# Author: Gregory Nutt <gnutt@nuttx.org>
#
# Redistribution and use in source and binary forms, with or without
@@ -33,11 +33,17 @@
#
############################################################################
# Add the math C files to the build
if ($(CONFIG_LIBM),y)
CSRCS += lib_rint.c lib_fixedmath.c lib_b16sin.c lib_b16cos.c lib_b16atan2.c
# Add the floating point math C files to the build
# Add the math directory to the build
CSRCS += lib_acos.c lib_asin.c lib_atan.c lib_atan2.c lib_ceil.c lib_cos.c
CSRCS += lib_cosh.c lib_exp.c lib_fabs.c lib_floor.c lib_fmod.c lib_frexp.c
CRRCS += lib_ldexp.c lib_log.c lib_log10.c lib_log2.c lib_modf.c lib_pow.c
CSRCS += lib_sin.c lib_sinh.c lib_sqrt.c lib_tan.c lib_tanh.c
# Add the floating point math directory to the build
DEPPATH += --dep-path math
VPATH += :math
endif
+100
View File
@@ -0,0 +1,100 @@
############################################################################
# nuttx/lib/math/Makefile
#
# Copyright (C) 2011-2012 Gregory Nutt. All rights reserved.
# Author: Gregory Nutt <gnutt@nuttx.org>
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
# 3. Neither the name NuttX nor the names of its contributors may be
# used to endorse or promote products derived from this software
# without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
# AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#
############################################################################
-include $(TOPDIR)/.config
-include $(TOPDIR)/Make.defs
include $(APPDIR)/Make.defs
# NSH Library
ASRCS =
CSRCS =
ifeq ($(CONFIG_MATHEXT_LIBRARY),y)
CSRCS += m_cos.c m_log2.c m_fmod.c m_sin.c m_fabs.c m_exp.c m_asin.c m_tan.c \
m_modf.c m_ldexp.c m_acos.c m_frexp.c m_atan.c m_log.c m_tanh.c m_ceil.c \
m_sinh.c m_sqrt.c m_pow.c m_cosh.c m_floor.c m_atan2.c m_log10.c
endif
AOBJS = $(ASRCS:.S=$(OBJEXT))
COBJS = $(CSRCS:.c=$(OBJEXT))
SRCS = $(ASRCS) $(CSRCS)
OBJS = $(AOBJS) $(COBJS)
ifeq ($(WINTOOL),y)
BIN = "${shell cygpath -w $(APPDIR)/libapps$(LIBEXT)}"
else
BIN = "$(APPDIR)/libapps$(LIBEXT)"
endif
ROOTDEPPATH = --dep-path .
VPATH =
# Build targets
all: .built
.PHONY: context .depend depend clean distclean
$(AOBJS): %$(OBJEXT): %.S
$(call ASSEMBLE, $<, $@)
$(COBJS): %$(OBJEXT): %.c
$(call COMPILE, $<, $@)
.built: $(OBJS)
@( for obj in $(OBJS) ; do \
$(call ARCHIVE, $(BIN), $${obj}); \
done ; )
@touch .built
context:
.depend: Makefile $(SRCS)
@$(MKDEP) $(ROOTDEPPATH) \
$(CC) -- $(CFLAGS) -- $(SRCS) >Make.dep
@touch $@
depend: .depend
clean:
@rm -f *.o *~ .*.swp .built
$(call CLEAN)
distclean: clean
@rm -f Make.dep .depend
-include Make.dep
+30
View File
@@ -0,0 +1,30 @@
/*
* Copyright (C) 2009, 2010 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>
float acosf(float x) {
return (M_PI_2 - asinf(x));
}
double acos(double x) {
return (M_PI_2 - asin(x));
}
long double acosl(long double x) {
return (M_PI_2 - asinl(x));
}
+88
View File
@@ -0,0 +1,88 @@
/*
* 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>
float asinf(float x) {
long double y, y_sin, y_cos;
y = 0;
while (1) {
y_sin = sinf(y);
y_cos = cosf(y);
if (y > M_PI_2 || y < -M_PI_2) {
y = fmodf(y, M_PI);
}
if (y_sin + DBL_EPSILON >= x && y_sin - DBL_EPSILON <= x) {
break;
}
y = y - (y_sin - x) / y_cos;
}
return y;
}
double asin(double x) {
long double y, y_sin, y_cos;
y = 0;
while (1) {
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;
}
return y;
}
long double asinl(long double x) {
long double y, y_sin, y_cos;
y = 0;
while (1) {
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;
}
return y;
}
+32
View File
@@ -0,0 +1,32 @@
/*
* 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 <stddef.h>
#include <stdint.h>
float atanf(float x) {
return asinf(x / sqrtf(x * x + 1));
}
double atan(double x) {
return asin(x / sqrt(x * x + 1));
}
long double atanl(long double x) {
return asinl(x / sqrtl(x * x + 1));
}
+120
View File
@@ -0,0 +1,120 @@
/*
* Copyright (C) 2009, 2010 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>
float atan2f(float y, float x) {
if (y == 0.0) {
if (x >= 0.0) {
return 0.0;
}
else {
return M_PI;
}
}
else if (y > 0.0) {
if (x == 0.0) {
return M_PI_2;
}
else if (x > 0.0) {
return atanf(y / x);
}
else {
return M_PI - atanf(y / x);
}
}
else {
if (x == 0.0) {
return M_PI + M_PI_2;
}
else if (x > 0.0) {
return 2 * M_PI - atanf(y / x);
}
else {
return M_PI + atanf(y / x);
}
}
}
double atan2(double y, double x) {
if (y == 0.0) {
if (x >= 0.0) {
return 0.0;
}
else {
return M_PI;
}
}
else if (y > 0.0) {
if (x == 0.0) {
return M_PI_2;
}
else if (x > 0.0) {
return atan(y / x);
}
else {
return M_PI - atan(y / x);
}
}
else {
if (x == 0.0) {
return M_PI + M_PI_2;
}
else if (x > 0.0) {
return 2 * M_PI - atan(y / x);
}
else {
return M_PI + atan(y / x);
}
}
}
long double atan2l(long double y, long double x) {
if (y == 0.0) {
if (x >= 0.0) {
return 0.0;
}
else {
return M_PI;
}
}
else if (y > 0.0) {
if (x == 0.0) {
return M_PI_2;
}
else if (x > 0.0) {
return atanl(y / x);
}
else {
return M_PI - atanl(y / x);
}
}
else {
if (x == 0.0) {
return M_PI + M_PI_2;
}
else if (x > 0.0) {
return 2 * M_PI - atanl(y / x);
}
else {
return M_PI + atanl(y / x);
}
}
}
+37
View File
@@ -0,0 +1,37 @@
/*
* 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>
float ceilf(float x) {
modff(x, &x);
if (x > 0.0) x += 1.0;
return x;
}
double ceil(double x) {
modf(x, &x);
if (x > 0.0) x += 1.0;
return x;
}
long double ceill(long double x) {
modfl(x, &x);
if (x > 0.0) x += 1.0;
return x;
}
+31
View File
@@ -0,0 +1,31 @@
/*
* Copyright (C) 2009, 2010 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 <stdint.h>
#include <float.h>
#include <apps/math.h>
float cosf(float x) {
return sinf(x + M_PI_2);
}
double cos(double x) {
return sin(x + M_PI_2);
}
long double cosl(long double x) {
return sinl(x + M_PI_2);
}
+33
View File
@@ -0,0 +1,33 @@
/*
* Copyright (C) 2009, 2010 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>
float coshf(float x) {
x = expf(x);
return ((x + (1.0 / x)) / 2.0);
}
double cosh(double x) {
x = exp(x);
return ((x + (1.0 / x)) / 2.0);
}
long double coshl(long double x) {
x = expl(x);
return ((x + (1.0 / x)) / 2.0);
}
+253
View File
@@ -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;
}
}
+30
View File
@@ -0,0 +1,30 @@
/*
* Copyright (C) 2009, 2010 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>
float fabsf(float x) {
return ((x < 0) ? -x : x);
}
double fabs(double x) {
return ((x < 0) ? -x : x);
}
long double fabsl(long double x) {
return ((x < 0) ? -x : x);
}
+37
View File
@@ -0,0 +1,37 @@
/*
* 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>
float floorf(float x) {
modff(x, &x);
if (x < 0.0) x -= 1.0;
return x;
}
double floor(double x) {
modf(x, &x);
if (x < 0.0) x -= 1.0;
return x;
}
long double floorl(long double x) {
modfl(x, &x);
if (x < 0.0) x -= 1.0;
return x;
}
+67
View File
@@ -0,0 +1,67 @@
/*
* 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>
/* If GCC/CLang builtins are available, use them */
#ifdef __GNUC__
float fmodf(float x, float div) {
return __builtin_fmodf(x, div);
}
double fmod(double x, double div) {
return __builtin_fmod(x, div);
}
long double fmodl(long double x, long double div) {
return __builtin_fmodl(x, div);
}
#else
float fmodf(float x, float div) {
float n0;
x /= div;
x = modff(x, &n0);
x *= div;
return x;
}
double fmod(double x, double div) {
double n0;
x /= div;
x = modf(x, &n0);
x *= div;
return x;
}
long double fmodl(long double x, long double div) {
long double n0;
x /= div;
x = modfl(x, &n0);
x *= div;
return x;
}
#endif
+34
View File
@@ -0,0 +1,34 @@
/*
* 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>
float frexpf(float x, int *exp) {
*exp = (int) ceilf(log2f(x));
return x / ldexpf(1.0, *exp);
}
double frexp(double x, int *exp) {
*exp = (int) ceil(log2(x));
return x / ldexp(1.0, *exp);
}
long double frexpl(long double x, int *exp) {
*exp = (int) ceill(log2(x));
return x / ldexpl(1.0, *exp);
}
+31
View File
@@ -0,0 +1,31 @@
/*
* 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>
float ldexpf(float x, int n) {
return (x * powf(2.0, (float) n));
}
double ldexp(double x, int n) {
return (x * pow(2.0, (double) n));
}
long double ldexpl(long double x, int n) {
return (x * powl(2.0, (long double) n));
}
+113
View File
@@ -0,0 +1,113 @@
/*
* Copyright (C) 2009, 2010 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>
float logf(float x) {
float y, y_old, ey, epsilon;
y = 0.0;
y_old = 1.0;
epsilon = FLT_EPSILON;
while (y > y_old + epsilon || y < y_old - epsilon) {
y_old = y;
ey = exp(y);
y -= (ey - x) / ey;
if (y > 700.0) {
y = 700.0;
}
if (y < -700.0) {
y = -700.0;
}
epsilon = (fabs(y) > 1.0) ? fabs(y) * FLT_EPSILON : FLT_EPSILON;
}
if (y == 700.0) {
return INFINITY;
}
if (y == -700.0) {
return INFINITY;
}
return y;
}
double log(double x) {
double y, y_old, ey, epsilon;
y = 0.0;
y_old = 1.0;
epsilon = DBL_EPSILON;
while (y > y_old + epsilon || y < y_old - epsilon) {
y_old = y;
ey = exp(y);
y -= (ey - x) / ey;
if (y > 700.0) {
y = 700.0;
}
if (y < -700.0) {
y = -700.0;
}
epsilon = (fabs(y) > 1.0) ? fabs(y) * DBL_EPSILON : DBL_EPSILON;
}
if (y == 700.0) {
return INFINITY;
}
if (y == -700.0) {
return INFINITY;
}
return y;
}
long double logl(long double x) {
long double y, y_old, ey, epsilon;
y = 0.0;
y_old = 1.0;
epsilon = 1e-6; //fixme
while (y > y_old + epsilon || y < y_old - epsilon) {
y_old = y;
ey = expl(y);
y -= (ey - x) / ey;
if (y > 700.0) {
y = 700.0;
}
if (y < -700.0) {
y = -700.0;
}
}
if (y == 700.0) {
return INFINITY;
}
if (y == -700.0) {
return INFINITY;
}
return y;
}
+30
View File
@@ -0,0 +1,30 @@
/*
* Copyright (C) 2009, 2010 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>
float log10f(float x) {
return (logf(x) / M_LN10);
}
double log10(double x) {
return (log(x) / M_LN10);
}
long double log10l(long double x) {
return (logl(x) / M_LN10);
}
+30
View File
@@ -0,0 +1,30 @@
/*
* Copyright (C) 2009, 2010 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>
float log2f(float x) {
return (logf(x) / M_LN2);
}
double log2(double x) {
return (log(x) / M_LN2);
}
long double log2l(long double x) {
return (logl(x) / M_LN2);
}
+65
View File
@@ -0,0 +1,65 @@
/*
* 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>
float modff(float x, float *iptr) {
if (fabsf(x) >= 8388608.0) {
*iptr = x;
return 0.0;
}
else if (fabs(x) < 1.0) {
*iptr = 0.0;
return x;
}
else {
*iptr = (float) (int) x;
return (x - *iptr);
}
}
double modf(double x, double *iptr) {
if (fabs(x) >= 4503599627370496.0) {
*iptr = x;
return 0.0;
}
else if (fabs(x) < 1.0) {
*iptr = 0.0;
return x;
}
else {
*iptr = (double) (int64_t) x;
return (x - *iptr);
}
}
long double modfl(long double x, long double *iptr) {
if (fabs(x) >= 4503599627370496.0) {
*iptr = x;
return 0.0;
}
else if (fabs(x) < 1.0) {
*iptr = 0.0;
return x;
}
else {
*iptr = (long double) (int64_t) x;
return (x - *iptr);
}
}
+30
View File
@@ -0,0 +1,30 @@
/*
* Copyright (C) 2009, 2010 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>
float powf(float b, float e) {
return expf(e * logf(b));
}
double pow(double b, double e) {
return exp(e * log(b));
}
long double powl(long double b, long double e) {
return expl(e * logl(b));
}
+152
View File
@@ -0,0 +1,152 @@
/*
* 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 <stdint.h>
#include <unistd.h>
#include <float.h>
#include <apps/math.h>
static float _flt_inv_fact[] = {
1.0 / 1.0, // 1 / 1!
1.0 / 6.0, // 1 / 3!
1.0 / 120.0, // 1 / 5!
1.0 / 5040.0, // 1 / 7!
1.0 / 362880.0, // 1 / 9!
1.0 / 39916800.0, // 1 / 11!
};
float sinf(float x) {
float x_squared;
float sin_x;
size_t i;
/* move x to [-pi, pi) */
x = fmodf(x, 2 * M_PI);
if (x >= M_PI) x -= 2 * M_PI;
if (x < -M_PI) x += 2 * M_PI;
/* move x to [-pi/2, pi/2) */
if (x >= M_PI_2) x = M_PI - x;
if (x < -M_PI_2) x = -M_PI - x;
x_squared = x * x;
sin_x = 0.0;
/* perform Taylor series approximation for sin(x) with six terms */
for (i = 0; i < 6; i++) {
if (i % 2 == 0) {
sin_x += x * _flt_inv_fact[i];
}
else {
sin_x -= x * _flt_inv_fact[i];
}
x *= x_squared;
}
return sin_x;
}
static double _dbl_inv_fact[] = {
1.0 / 1.0, // 1 / 1!
1.0 / 6.0, // 1 / 3!
1.0 / 120.0, // 1 / 5!
1.0 / 5040.0, // 1 / 7!
1.0 / 362880.0, // 1 / 9!
1.0 / 39916800.0, // 1 / 11!
1.0 / 6227020800.0, // 1 / 13!
1.0 / 1307674368000.0, // 1 / 15!
1.0 / 355687428096000.0, // 1 / 17!
1.0 / 121645100408832000.0, // 1 / 19!
};
double sin(double x) {
double x_squared;
double sin_x;
size_t i;
/* move x to [-pi, pi) */
x = fmod(x, 2 * M_PI);
if (x >= M_PI) x -= 2 * M_PI;
if (x < -M_PI) x += 2 * M_PI;
/* move x to [-pi/2, pi/2) */
if (x >= M_PI_2) x = M_PI - x;
if (x < -M_PI_2) x = -M_PI - x;
x_squared = x * x;
sin_x = 0.0;
/* perform Taylor series approximation for sin(x) with ten terms */
for (i = 0; i < 10; i++) {
if (i % 2 == 0) {
sin_x += x * _dbl_inv_fact[i];
}
else {
sin_x -= x * _dbl_inv_fact[i];
}
x *= x_squared;
}
return sin_x;
}
static long double _ldbl_inv_fact[] = {
1.0 / 1.0, // 1 / 1!
1.0 / 6.0, // 1 / 3!
1.0 / 120.0, // 1 / 5!
1.0 / 5040.0, // 1 / 7!
1.0 / 362880.0, // 1 / 9!
1.0 / 39916800.0, // 1 / 11!
1.0 / 6227020800.0, // 1 / 13!
1.0 / 1307674368000.0, // 1 / 15!
1.0 / 355687428096000.0, // 1 / 17!
1.0 / 121645100408832000.0, // 1 / 19!
};
long double sinl(long double x) {
long double x_squared;
long double sin_x;
size_t i;
/* move x to [-pi, pi) */
x = fmodl(x, 2 * M_PI);
if (x >= M_PI) x -= 2 * M_PI;
if (x < -M_PI) x += 2 * M_PI;
/* move x to [-pi/2, pi/2) */
if (x >= M_PI_2) x = M_PI - x;
if (x < -M_PI_2) x = -M_PI - x;
x_squared = x * x;
sin_x = 0.0;
/* perform Taylor series approximation for sin(x) with ten terms */
for (i = 0; i < 10; i++) {
if (i % 2 == 0) {
sin_x += x * _ldbl_inv_fact[i];
}
else {
sin_x -= x * _ldbl_inv_fact[i];
}
x *= x_squared;
}
return sin_x;
}
+33
View File
@@ -0,0 +1,33 @@
/*
* Copyright (C) 2009, 2010 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>
float sinhf(float x) {
x = expf(x);
return ((x - (1.0 / x)) / 2.0);
}
double sinh(double x) {
x = exp(x);
return ((x - (1.0 / x)) / 2.0);
}
long double sinhl(long double x) {
x = expl(x);
return ((x - (1.0 / x)) / 2.0);
}
+116
View File
@@ -0,0 +1,116 @@
/*
* 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 <stdint.h>
#include <float.h>
#include <errno.h>
#include <apps/math.h>
static float __sqrt_approx(float x) {
int32_t i;
// floats + bit manipulation = +inf fun!
i = *((int32_t*) &x);
i = 0x1FC00000 + (i >> 1);
x = *((float*) &i);
return x;
}
float sqrtf(float x) {
float y;
// filter out invalid/trivial inputs
if (x < 0.0) { errno = EDOM; return NAN; }
if (isnan(x)) return NAN;
if (isinf(x)) return INFINITY;
if (x == 0.0) return 0.0;
// guess square root (using bit manipulation)
y = __sqrt_approx(x);
// perform three iterations of approximation
// this number (3) is definitely optimal
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
return y;
}
double sqrt(double x) {
long double y, y1;
// filter out invalid/trivial inputs
if (x < 0.0) { errno = EDOM; return NAN; }
if (isnan(x)) return NAN;
if (isinf(x)) return INFINITY;
if (x == 0.0) return 0.0;
// guess square root (using bit manipulation)
y = __sqrt_approx(x);
// perform four iterations of approximation
// this number (4) is definitely optimal
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
// if guess was terribe (out of range of float)
// repeat approximation until convergence
if (y * y < x - 1.0 || y * y > x + 1.0) {
y1 = -1.0;
while (y != y1) {
y1 = y;
y = 0.5 * (y + x / y);
}
}
return y;
}
long double sqrtl(long double x) {
long double y, y1;
// filter out invalid/trivial inputs
if (x < 0.0) { errno = EDOM; return NAN; }
if (isnan(x)) return NAN;
if (isinf(x)) return INFINITY;
if (x == 0.0) return 0.0;
// guess square root (using bit manipulation)
y = __sqrt_approx(x);
// perform four iterations of approximation
// this number (4) is definitely optimal
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
y = 0.5 * (y + x / y);
// if guess was terribe (out of range of float)
// repeat approximation until convergence
if (y * y < x - 1.0 || y * y > x + 1.0) {
y1 = -1.0;
while (y != y1) {
y1 = y;
y = 0.5 * (y + x / y);
}
}
return y;
}
+31
View File
@@ -0,0 +1,31 @@
/*
* Copyright (C) 2009, 2010 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 <stdint.h>
#include <float.h>
#include <apps/math.h>
float tanf(float x) {
return (sinf(x) / cosf(x));
}
double tan(double x) {
return (sin(x) / cos(x));
}
long double tanl(long double x) {
return (sinl(x) / cosl(x));
}
+39
View File
@@ -0,0 +1,39 @@
/*
* Copyright (C) 2009, 2010 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>
float tanhf(float x) {
float x0 = expf(x);
float x1 = 1.0 / x0;
return ((x0 + x1) / (x0 - x1));
}
double tanh(double x) {
double x0 = exp(x);
double x1 = 1.0 / x0;
return ((x0 + x1) / (x0 - x1));
}
long double tanhl(long double x) {
long double x0 = exp(x);
long double x1 = 1.0 / x0;
return ((x0 + x1) / (x0 - x1));
}