round half-integers to even instead of away from zero

This changes the behaviour to match r7rs (round x) instead of C round(x).

An answer to https://stackoverflow.com/questions/32746523/ieee-754-compliant-round-half-to-even
suggests using remainder(). The following will work if FE_TONEAREST is defined, but C11
requires FE_TONEAREST to be defined if and only if the implemenetation supports it in
fegetround() and fesetround() [Draft N1570]. On the other hand, remainder() must be defined.
C23 will have roundeven(), but this is not yet available on all platforms.
The behaviour of remainder is described in Draft N1570, page 254, footnote 239.

Alternative implementation:

  double round_to_nearest_even(double x)
  {
  #pragma STDC FENV_ACCESS ON
    int mode;
    double nearest;
    mode = fegetround();
    fesetround(FE_TONEAREST);
    nearest = nearbyint(x);
    fesetround(mode);
  #pragma STDC FENV_ACCESS OFF
    return nearest;
  }

3 files changed, 8 insertions, 2 deletions

diff --git a/include/cyclone/runtime.h b/include/cyclone/runtime.h
index 82608b30..432be667 100644
--- a/include/cyclone/runtime.h
+++ b/include/cyclone/runtime.h
@@ -504,6 +504,7 @@ int Cyc_have_mstreams();
   } \
   return_closcall1(data, cont, &d)
 
+double round_to_nearest_even(double);
 void Cyc_exact(void *data, object cont, object z);
 object Cyc_exact_no_cps(void *data, object ptr, object z);
 
diff --git a/runtime.c b/runtime.c
index 3448cecb..59143b47 100644
--- a/runtime.c
+++ b/runtime.c
@@ -8765,6 +8765,11 @@ int num2ratio(double x, double *numerator, double *denominator)
   return 0;
 }
 
+double round_to_nearest_even(double x)
+{
+ return x-remainder(x,1.0);
+}
+
 /**
  * Receive a Scheme number and pass requested portion of a rational number to
  * the continuation `cont`. Pass numerator if `numerator` is true, else the
diff --git a/scheme/base.sld b/scheme/base.sld
index 669b8cde..ae585953 100644
--- a/scheme/base.sld
+++ b/scheme/base.sld
@@ -1372,9 +1372,9 @@
     " return_double_op_no_cps(data, ptr, trunc, z);")
   (define-c round
     "(void *data, int argc, closure _, object k, object z)"
-    " return_double_op(data, k, round, z); "
+    " return_double_op(data, k, round_to_nearest_even, z); "
     "(void *data, object ptr, object z)"
-    " return_double_op_no_cps(data, ptr, round, z);")
+    " return_double_op_no_cps(data, ptr, round_to_nearest_even, z);")
   (define-c exact
     "(void *data, int argc, closure _, object k, object z)"
     " Cyc_exact(data, k, z); "