summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorManuel Novoa III <mjn3@codepoet.org>2001-05-22 15:03:49 +0000
committerManuel Novoa III <mjn3@codepoet.org>2001-05-22 15:03:49 +0000
commitca3c705bbb793559bc5adc8fd4e0bb7e4b350f2e (patch)
tree458362b6d010adda88da063a7962486b6630aabd
parent54265bebea80fcb705f88209d30ae5d52c9bc945 (diff)
Added file for non-Cephes double routines; currently only fmod and modf.
-rw-r--r--libm/double/Makefile2
-rw-r--r--libm/double/noncephes.c127
2 files changed, 128 insertions, 1 deletions
diff --git a/libm/double/Makefile b/libm/double/Makefile
index b64116df8..e594b0f85 100644
--- a/libm/double/Makefile
+++ b/libm/double/Makefile
@@ -35,7 +35,7 @@ CSRC=acosh.c airy.c asin.c asinh.c atan.c atanh.c bdtr.c beta.c \
polevl.c polmisc.c polylog.c polyn.c pow.c powi.c psi.c rgamma.c round.c \
shichi.c sici.c sin.c sindg.c sinh.c spence.c stdtr.c struve.c \
tan.c tandg.c tanh.c unity.c yn.c zeta.c zetac.c \
- sqrt.c floor.c setprec.c mtherr.c
+ sqrt.c floor.c setprec.c mtherr.c noncephes.c
COBJS=$(patsubst %.c,%.o, $(CSRC))
diff --git a/libm/double/noncephes.c b/libm/double/noncephes.c
new file mode 100644
index 000000000..72f129d55
--- /dev/null
+++ b/libm/double/noncephes.c
@@ -0,0 +1,127 @@
+/*
+ * This file contains math functions missing from the Cephes library.
+ *
+ * May 22, 2001 Manuel Novoa III
+ *
+ * Added modf and fmod.
+ *
+ * TODO:
+ * Break out functions into seperate object files as is done
+ * by (for example) stdio. Also do this with cephes files.
+ */
+
+#include <math.h>
+#include <errno.h>
+
+#undef UNK
+
+/* Set this to nonzero to enable a couple of shortcut tests in fmod. */
+#define SPEED_OVER_SIZE 0
+
+/**********************************************************************/
+
+double modf(double x, double *iptr)
+{
+ double y;
+
+#ifdef UNK
+ mtherr( "modf", DOMAIN );
+ *iptr = NAN;
+ return NAN;
+#endif
+
+#ifdef NANS
+ if( isnan(x) ) {
+ *iptr = x;
+ return x;
+ }
+#endif
+
+#ifdef INFINITIES
+ if(!isfinite(x)) {
+ *iptr = x; /* Matches glibc, but returning NAN */
+ return 0; /* makes more sense to me... */
+ }
+#endif
+
+ if (x < 0) { /* Round towards 0. */
+ y = ceil(x);
+ } else {
+ y = floor(x);
+ }
+
+ *iptr = y;
+ return x - y;
+}
+
+/**********************************************************************/
+
+extern double NAN;
+
+double fmod(double x, double y)
+{
+ double z;
+ int negative, ex, ey;
+
+#ifdef UNK
+ mtherr( "fmod", DOMAIN );
+ return NAN;
+#endif
+
+#ifdef NANS
+ if( isnan(x) || isnan(y) ) {
+ errno = EDOM;
+ return NAN;
+ }
+#endif
+
+ if (y == 0) {
+ errno = EDOM;
+ return NAN;
+ }
+
+#ifdef INFINITIES
+ if(!isfinite(x)) {
+ errno = EDOM;
+ return NAN;
+ }
+
+#if SPEED_OVER_SIZE
+ if(!isfinite(y)) {
+ return x;
+ }
+#endif
+#endif
+
+#if SPEED_OVER_SIZE
+ if (x == 0) {
+ return 0;
+ }
+#endif
+
+ negative = 0;
+ if (x < 0) {
+ negative = 1;
+ x = -x;
+ }
+
+ if (y < 0) {
+ y = -y;
+ }
+
+ frexp(y,&ey);
+ while (x >= y) {
+ frexp(x,&ex);
+ z = ldexp(y,ex-ey);
+ if (z > x) {
+ z /= 2;
+ }
+ x -= z;
+ }
+
+ if (negative) {
+ return -x;
+ } else {
+ return x;
+ }
+}