From cca63462fcc1f9c9f0b1dcdb0ea194458718f8d6 Mon Sep 17 00:00:00 2001
From: bumbread <thebumboni@gmail.com>
Date: Thu, 16 Jun 2022 21:51:55 +1100
Subject: [PATCH] math generic exp function

---
 bake.cmd               |  2 +-
 code/math/exp/exp.c    |  8 --------
 code/math/generic.c    | 25 +++++++++++++++++++++++++
 code/math/generic.h    | 42 ++++++++++++++++++++++++++++++++++++++++++
 code/stdio/fmt_print.h |  1 +
 inc/math.h             |  2 +-
 test/test_math.c       | 14 ++++++++++++++
 7 files changed, 84 insertions(+), 10 deletions(-)
 delete mode 100644 code/math/exp/exp.c
 create mode 100644 code/math/generic.c
 create mode 100644 code/math/generic.h

diff --git a/bake.cmd b/bake.cmd
index 7631633..a243b17 100644
--- a/bake.cmd
+++ b/bake.cmd
@@ -51,5 +51,5 @@ del build\*.obj
 
 :skip_crt_compilation
 echo Compiling test..
-clang -fno-builtin test\test_wctype.c ciabatta.lib -std=c11 -lkernel32 -luser32 -lshell32 -nostdlib %CIABATTA_OPTIONS%
+clang -fno-builtin test\test_math.c ciabatta.lib -std=c11 -lkernel32 -luser32 -lshell32 -nostdlib %CIABATTA_OPTIONS%
 ::cl test\test_math.c /Iinc -D_CRT_SECURE_NO_WARNINGS /Z7 /link ciabatta.lib kernel32.lib user32.lib shell32.lib -nostdlib -nodefaultlibs
diff --git a/code/math/exp/exp.c b/code/math/exp/exp.c
deleted file mode 100644
index e0a9a37..0000000
--- a/code/math/exp/exp.c
+++ /dev/null
@@ -1,8 +0,0 @@
-
-#include <math.h>
-#include <stdint.h>
-#include <float.h>
-#include <fenv.h>
-
-#define rln2 1.4426950408889634073599246810018921374266459541529859341354494069
-
diff --git a/code/math/generic.c b/code/math/generic.c
new file mode 100644
index 0000000..b61b39b
--- /dev/null
+++ b/code/math/generic.c
@@ -0,0 +1,25 @@
+
+#include <math.h>
+#include <fenv.h>
+#include <errno.h>
+#include <stdint.h>
+
+#define ln2 0.69314718055994530941723212145817
+
+#define ftype float
+#define suffix(name) name ## f
+#include "generic.h"
+#undef ftype
+#undef suffix
+
+#define ftype double
+#define suffix(name) name
+#include "generic.h"
+#undef ftype
+#undef suffix
+
+#define ftype long double
+#define suffix(name) name ## l
+#include "generic.h"
+#undef ftype
+#undef suffix
diff --git a/code/math/generic.h b/code/math/generic.h
new file mode 100644
index 0000000..872f5c0
--- /dev/null
+++ b/code/math/generic.h
@@ -0,0 +1,42 @@
+
+ftype suffix(exp)(ftype x) {
+    if(isnan(x)) return NAN;
+    if(x > 0 && isinf(x)) return +INFINITY;
+    if(x < 0 && isinf(x)) return +0.0;
+    if(x == 0) return 1.0;
+    if(x > 709.8) {
+        #if math_errhandling & MATH_ERREXCEPT
+            feraiseexcept(FE_OVERFLOW);
+        #endif
+        #if math_errhandling & MATH_ERRNO
+            errno = ERANGE;
+        #endif
+        return +INFINITY;
+    }
+    if(x < -708.4) {
+        #if math_errhandling & MATH_ERREXCEPT
+            feraiseexcept(FE_OVERFLOW);
+        #endif
+        #if math_errhandling & MATH_ERRNO
+            errno = ERANGE;
+        #endif
+        return 0;
+    }
+    ftype e = 1.0;
+    ftype xp = 1.0;
+    ftype f = 1;
+    for(uint64_t i = 1; i != 10; ++i) {
+        f *= i;
+        xp *= x;
+        e += xp / f;
+    }
+    return e;
+}
+
+ftype suffix(exp2)(ftype x) {
+    return suffix(exp)(x * ln2);
+}
+
+ftype suffix(expm1)(ftype x) {
+    return suffix(exp)(x) - 1.0;
+}
diff --git a/code/stdio/fmt_print.h b/code/stdio/fmt_print.h
index 9e846d0..5efd24f 100644
--- a/code/stdio/fmt_print.h
+++ b/code/stdio/fmt_print.h
@@ -4,6 +4,7 @@
 #include <stdint.h>
 #include <limits.h>
 #include <string.h>
+#include <math.h>
 
 // This stuff is kinda related to what's going on in this file, so I left it
 // in a rather ugly manner here.
diff --git a/inc/math.h b/inc/math.h
index ad53ca8..d93c220 100644
--- a/inc/math.h
+++ b/inc/math.h
@@ -20,7 +20,7 @@ typedef double double_t;
 #define MATH_ERRNO     1
 #define MATH_ERREXCEPT 2
 
-#define math_errhandling MATH_ERRNO
+#define math_errhandling MATH_ERRNO | MATH_ERREXCEPT
 
 // Classification
 #define FP_INFINITE  0
diff --git a/test/test_math.c b/test/test_math.c
index 1c3d49c..ab73bb0 100644
--- a/test/test_math.c
+++ b/test/test_math.c
@@ -4,6 +4,7 @@
 #include <float.h>
 #include <fenv.h>
 #include <inttypes.h>
+#include <errno.h>
 
 const char *show_classification(double x) {
     switch(fpclassify(x)) {
@@ -103,5 +104,18 @@ int main() {
     printf("lround(LONG_MAX+1.5) = %ld\n", lround(LONG_MAX+1.5));
     if(fetestexcept(FE_INVALID)) printf("    FE_INVALID was raised\n");
 
+    printf("\n\n=== exp === \n");
+    printf("exp(1) = %f\n", exp(1));
+    printf("FV of $100, continuously compounded at 3%% for 1 year = %f\n",
+            100*exp(0.03));
+    // special values
+    printf("exp(-0) = %f\n", exp(-0.0));
+    printf("exp(-Inf) = %f\n", exp(-INFINITY));
+    //error handling
+    errno = 0; feclearexcept(FE_ALL_EXCEPT);
+    printf("exp(710) = %f\n", exp(710));
+    if(errno == ERANGE) printf("    errno == ERANGE\n");
+    if(fetestexcept(FE_OVERFLOW)) printf("    FE_OVERFLOW raised\n");
+
     return 0;
 }