3 files changed, 164 insertions, 0 deletions
diff --git a/test/monniaux/math/exceptions.c b/test/monniaux/math/exceptions.c
new file mode 100644
index 00000000..72107066
--- /dev/null
+++ b/test/monniaux/math/exceptions.c
@@ -0,0 +1,84 @@
+#include <stdio.h>
+#include <fenv.h>
+#include <float.h>
+
+#pragma STDC FENV_ACCESS ON
+
+#if defined(__K1C__) && !defined(__COMPCERT__)
+int fetestexcept(int excepts) {
+  int mask = (K1_SFR_CS_IO_MASK | K1_SFR_CS_DZ_MASK | K1_SFR_CS_OV_MASK | K1_SFR_CS_UN_MASK | K1_SFR_CS_IN_MASK) & excepts;
+  unsigned long long cs = __builtin_k1_get(K1_SFR_CS);
+  return cs & mask;
+}
+
+int feclearexcept(int excepts) {
+  int mask = (K1_SFR_CS_IO_MASK | K1_SFR_CS_DZ_MASK | K1_SFR_CS_OV_MASK | K1_SFR_CS_UN_MASK | K1_SFR_CS_IN_MASK) & excepts;
+  __builtin_k1_wfxl(K1_SFR_CS, mask);
+  return 0;
+}
+#endif
+
+double add(double x, double y) {
+  return x+y;
+}
+
+double mul(double x, double y) {
+  return x*y;
+}
+
+float double2float(double x) {
+  return x;
+}
+
+float uint2float(unsigned x) {
+  return x;
+}
+
+double ulong2double(unsigned long x) {
+  return x;
+}
+
+unsigned double2uint(double x) {
+  return x;
+}
+
+int main() {
+  printf("%x\n", fetestexcept(FE_ALL_EXCEPT));
+
+  double v1 = add(3.0, 0.1);
+  printf("%x\n", fetestexcept(FE_ALL_EXCEPT));
+  feclearexcept(FE_INEXACT);
+
+  printf("%x\n", fetestexcept(FE_ALL_EXCEPT));
+  double v2 = mul(DBL_MAX, DBL_MAX);
+  printf("%g %x\n", v2, fetestexcept(FE_ALL_EXCEPT));
+  feclearexcept(FE_ALL_EXCEPT);
+
+  double v3 = mul(DBL_MIN, DBL_MIN);
+  printf("%g %x\n", v3, fetestexcept(FE_ALL_EXCEPT));
+  feclearexcept(FE_ALL_EXCEPT);
+
+  double v4 = double2float(DBL_MAX);
+  printf("%g %x\n", v4, fetestexcept(FE_ALL_EXCEPT));
+  feclearexcept(FE_ALL_EXCEPT);
+
+  float v5 = uint2float(0xC07FDFFFU);
+  printf("%g %x\n", v5, fetestexcept(FE_ALL_EXCEPT)); // BUG 0 should have INEXACT
+  feclearexcept(FE_ALL_EXCEPT);
+
+  double v6 = ulong2double(0x11217763AFF77C7CUL);
+  printf("%g %x\n", v6, fetestexcept(FE_ALL_EXCEPT)); // BUG 0 should have INEXACT
+  feclearexcept(FE_ALL_EXCEPT);
+
+  unsigned v7 = double2uint(-0.25); // softfloat says "0 and inexact" but here we have "0 and overflow" (due to negative input for unsigned?)
+  printf("%u %x\n", v7, fetestexcept(FE_ALL_EXCEPT));
+  feclearexcept(FE_ALL_EXCEPT);
+
+  // +41F.307672C5496EF
+  double d8 = 0x1.307672C5496EFp32;
+  unsigned v8 = double2uint(d8);
+  printf("%g %x %x\n", d8, v8, fetestexcept(FE_ALL_EXCEPT));
+  // BUG reports 307672C5 and inexact, but should report overflow
+  // bug comes from x86 https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89175
+  feclearexcept(FE_ALL_EXCEPT);
+}
diff --git a/test/monniaux/math/rounding.c b/test/monniaux/math/rounding.c
new file mode 100644
index 00000000..c2ce85e3
--- /dev/null
+++ b/test/monniaux/math/rounding.c
@@ -0,0 +1,36 @@
+#include <stdio.h>
+#include <fenv.h>
+
+#ifdef __K1C__
+#include <mppa_bare_runtime/k1c/registers.h>
+int fesetround(int rounding_mode) {
+  if (rounding_mode < 0 || rounding_mode > 3) return 1;
+  unsigned long long cs = __builtin_k1_get(K1_SFR_CS);
+  cs = (cs & ~(3 << 16)) | (rounding_mode << 16);
+  __builtin_k1_set(K1_SFR_CS, cs);
+  return 0;
+}
+
+int fegetround(void) {
+  unsigned long long cs = __builtin_k1_get(K1_SFR_CS);
+  return (cs >> 16) & 3;
+}
+#endif
+
+#pragma STDC FENV_ACCESS ON
+
+double add(double x) {
+  return x+0.1;
+}
+
+int main() {
+  if (fesetround(FE_DOWNWARD)) {
+    printf("fesetround(FE_DOWNWARD) unsuccessful\n");
+  }
+  double low = add(3.0);
+  if (fesetround(FE_UPWARD)) {
+    printf("fesetround(FE_UPWARD) unsuccessful\n");
+  }
+  double high = add(3.0);
+  printf("%d %a %d\n", low==high, high-low, fegetround());
+}
diff --git a/test/monniaux/math/test_gsl_pow.c b/test/monniaux/math/test_gsl_pow.c
new file mode 100644
index 00000000..54e53889
--- /dev/null
+++ b/test/monniaux/math/test_gsl_pow.c
@@ -0,0 +1,44 @@
+#include <math.h>
+#include <stdio.h>
+
+double gsl_pow_uint(double x, unsigned int n);
+
+double gsl_pow_int(double x, int n)
+{
+  unsigned int un;
+
+  if(n < 0) {
+    x = 1.0/x;
+    un = -n;
+  } else {
+    un = n;
+  }
+
+  return gsl_pow_uint(x, un);
+}
+
+double gsl_pow_uint(double x, unsigned int n)
+{
+  double value = 1.0;
+
+  /* repeated squaring method 
+   * returns 0.0^0 = 1.0, so continuous in x
+   */
+  do {
+     if(n & 1) value *= x;  /* for n odd */
+     n >>= 1;
+     x *= x;
+  } while (n);
+
+  return value;
+}
+
+int main() {
+  double y, y_expected;
+  for (int n = -9; n < 10; n++) {
+    y = gsl_pow_int (-3.14, n);
+    y_expected = pow (-3.14, n);
+    printf("%d %g %g\n", n, y, y_expected);
+  }
+  return 0;
+}