fenv.h

bake.cmd

@@ -51,4 +51,4 @@ del build\*.obj
 
 :skip_crt_compilation
 echo Compiling test..
-clang test\test5.c ciabatta.lib -std=c11 -lkernel32 -luser32 -lshell32 -nostdlib %CIABATTA_OPTIONS%
+clang test\test6.c ciabatta.lib -std=c11 -lkernel32 -luser32 -lshell32 -nostdlib %CIABATTA_OPTIONS%

code/fenv.c

@@ -0,0 +1,104 @@
+
+#include <fenv.h>
+#include <intrin.h>
+
+#define fe_masks(excepts) (((fexcept_t)(excepts)) << 7)
+#define fe_flags(excepts) ((fexcept_t)(excepts))
+
+fenv_t _fe_dfl_env = 0x1f80; // Based (on my machine)
+
+int feclearexcept(int excepts)
+{
+    if((excepts & FE_ALL_EXCEPT) != excepts) {
+        return 1;
+    }
+    if(excepts == 0) {
+        return 0;
+    }
+    fexcept_t csr = _mm_getcsr();
+    csr |= fe_masks(excepts);
+    _mm_setcsr(csr);
+    return 0;
+}
+
+int fegetexceptflag(fexcept_t *flagp, int excepts) {
+    if((excepts & FE_ALL_EXCEPT) != excepts) {
+        return 1;
+    }
+    *flagp = fe_flags(excepts);
+    return 0;
+}
+
+int feraiseexcept(int excepts) {
+    if((excepts & FE_ALL_EXCEPT) != excepts) {
+        return 1;
+    }
+    if(excepts == 0) {
+        return 0;
+    }
+    fexcept_t csr = _mm_getcsr();
+    csr |= fe_flags(excepts);
+    _mm_setcsr(csr);
+    return 0;
+}
+
+int fesetexceptflag(const fexcept_t *flagp, int excepts) {
+    if((excepts & FE_ALL_EXCEPT) != excepts) {
+        return 1;
+    }
+    if(excepts == 0) {
+        return 0;
+    }
+    fexcept_t flags = *flagp;
+    fexcept_t csr = _mm_getcsr();
+    csr |= flags;
+    _mm_setcsr(csr);
+    return 0;
+}
+
+int fetestexcept(int excepts) {
+    fexcept_t csr = _mm_getcsr();
+    fexcept_t flags = fe_flags(excepts);
+    return (int)(csr & flags);
+}
+
+int fegetround(void) {
+    fexcept_t csr = _mm_getcsr();
+    int round = (csr >> 13) & 0x3;
+    return round;
+}
+
+int fesetround(int round) {
+    if(!(0 <= round && round < 4)) {
+        return 1;
+    }
+    fexcept_t csr = _mm_getcsr();
+    csr &= ~(0x3 << 13);
+    csr |= round << 13;
+    _mm_setcsr(csr);
+    return 0;
+}
+
+int fegetenv(fenv_t *env) {
+    fenv_t csr = _mm_getcsr();
+    *env = csr;
+    return 0;
+}
+
+int fesetenv(fenv_t *env) {
+    _mm_setcsr(*env);
+    return 1;
+}
+
+int feholdexcept(fenv_t *envp) {
+    fegetenv(envp);
+    feclearexcept(FE_ALL_EXCEPT);
+    return 0;
+}
+
+int feupdateenv(fenv_t const *envp) {
+    int excepts = _mm_getcsr() & FE_ALL_EXCEPT;
+    _mm_setcsr(*envp);
+    feraiseexcept(excepts);
+    return 0;
+}

code/math/pow.h

@@ -93,3 +93,11 @@ ftype suffix(sqrt)(ftype x) {
     bits = b_cons(0, exp, man);
     return suffix(f_frombits)(bits);
 }
+
+ftype suffix(hypot)(ftype x, ftype y)
+{
+    if(isinf(x) || isinf(y)) {
+        return INFINITY;
+    }
+    return suffix(sqrt)(x*x + y*y);
+}

inc/fenv.h

@@ -1,26 +1,39 @@
 
 #pragma once
 
-typedef unsigned long fexcept_t;
+typedef unsigned fexcept_t;
-typedef struct fenv_t fenv_t;
+typedef unsigned fenv_t;
 
-#define FE_DIVBYZERO 0x04
+#define FE_INVALID   (1 << 0)
-#define FE_INEXACT   0x20
+#define FE_DIVBYZERO (1 << 2)
-#define FE_INVALID   0x01
+#define FE_OVERFLOW  (1 << 3)
-#define FE_OVERFLOW  0x08
+#define FE_UNDERFLOW (1 << 4)
-#define FE_UNDERFLOW 0x10
+#define FE_INEXACT   (1 << 5)
-#define FE_ALL_EXCEPT (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW)
+#define FE_ALL_EXCEPT   \
+    ( FE_INVALID        \
+    | FE_DIVBYZERO      \
+    | FE_OVERFLOW       \
+    | FE_UNDERFLOW      \
+    | FE_INEXACT        )
 
-#define FE_TONEAREST  0x000
+#define FE_TONEAREST  0x00
-#define FE_UPWARD     0x200
+#define FE_DOWNWARD   0x01
-#define FE_DOWNWARD   0x100
+#define FE_UPWARD     0x02
-#define FE_TOWARDZERO 0x300
+#define FE_TOWARDZERO 0x03
 
-// TODO: implement this
+extern fenv_t _fe_dfl_env;
-#define FE_DFL_ENV ((fenv_t*)0)
+#define FE_DFL_ENV (&_fe_dfl_env)
 
 int feclearexcept(int excepts);
 int fegetexceptflag(fexcept_t *flagp, int excepts);
 int feraiseexcept(int excepts);
 int fesetexceptflag(const fexcept_t *flagp, int excepts);
 int fetestexcept(int excepts);
+
+int fegetround(void);
+int fesetround(int round);
+
+int fegetenv(fenv_t *env);
+int fesetenv(fenv_t *env);
+int feholdexcept(fenv_t *envp);
+int feupdateenv(fenv_t const *envp);

test/test6.c

@@ -0,0 +1,47 @@
+#include <stdio.h>
+#include <math.h>
+#include <fenv.h>
+
+#include <emmintrin.h>
+
+#pragma STDC FENV_ACCESS ON
+
+static inline double rint (double const x) {
+    return (double)_mm_cvtsd_si32(_mm_load_sd(&x));
+}
+
+void show_fe_current_rounding_direction(void)
+{
+    printf("current rounding direction:  ");
+    switch (fegetround()) {
+           case FE_TONEAREST:  printf ("FE_TONEAREST");  break;
+           case FE_DOWNWARD:   printf ("FE_DOWNWARD");   break;
+           case FE_UPWARD:     printf ("FE_UPWARD");     break;
+           case FE_TOWARDZERO: printf ("FE_TOWARDZERO"); break;
+           default:            printf ("unknown");
+    };
+    printf("\n");
+}
+ 
+int main(void)
+{
+    /* Default rounding direction */
+    show_fe_current_rounding_direction();
+    printf("+11.5 -> %f\n", rint(+11.5)); /* midway between two integers */
+    printf("+12.5 -> %f\n", rint(+12.5)); /* midway between two integers */
+ 
+    /* Save current rounding direction. */
+    int curr_direction = fegetround();
+ 
+    /* Temporarily change current rounding direction. */
+    fesetround(FE_DOWNWARD);
+    show_fe_current_rounding_direction();
+    printf("+11.5 -> %f\n", rint(+11.5));
+    printf("+12.5 -> %f\n", rint(+12.5));
+ 
+    /* Restore default rounding direction. */
+    fesetround(curr_direction);
+    show_fe_current_rounding_direction(); 
+ 
+    return 0;
+}