scm_round certain wrong results

[Kevin Ryde]

        * numbers.c (scm_round): Test for x already an integer, to avoid bad
        rounding in x+0.5 when x is a big value already an integer.  In
        various modes x+0.5 could give an adjacent integer, leading to that as
        the result, when we really just wanted x itself.

        * standalone/test-round.c: New file, exercising scm_round.
        * standalone/Makefile.am: Add it.

        * configure.in (AC_CHECK_FUNCS): Add fesetround.
        (AC_CHECK_HEADERS): Add fenv.h.

The test program exercises the bad cases.  In the usual default
nearest-even hardware rounding I think x = +/- (2^53-1) is the only
bad case.

But if an application has gone to for instance round-upwards mode then
lots of other big values come out wrong too.  I believe it's worth
scm_round ensuring it's not influenced by the hardware rounding mode,
even if that's not something guile normally fiddles with.

This would be for the 1.6 branch too.  scm_round is used for `round'
there, you can see a problem for instance with

        (define x #i#x1FFFFFFFFFFFFF)
        (integer? x)     => #t
        (= x (round x))  => #f

If you're doing this on an i386 or m68k you need numbers.c compiled
unoptimized, since optimized will hold values in the coprocessor,
giving enough extra precision to avoid this case (but not others).  On
other cpus I think it will happen always.

--- numbers.c.~1.236.~  2004-04-18 07:55:31.000000000 +1000
+++ numbers.c   2004-04-18 08:29:42.000000000 +1000
@@ -4884,11 +4884,42 @@
 #endif
 }
 
+/* scm_round is done using floor(x+0.5) to round to nearest with half-way
+   (ie. when x is an integer plus 0.5) cases going upwards.  Then half-way
+   cases are identified and adjusted down if the round-upwards didn't give
+   the desired even integer.
+
+   "plus_half == result" identifies a half-way case.  If plus_half, which is
+   x + 0.5, is an integer then x must be an integer plus 0.5.
+
+   An odd "result" value is identified with result/2 != floor(result/2).
+   This is done with plus_half, since that value is ready for use sooner in
+   a pipelined cpu, and we're already requiring plus_half == result.
+
+   Note however that we need to be careful when x is big and already an
+   integer.  In that case "x+0.5" may round to an adjacent integer, causing
+   us to return such a value, incorrectly.  For instance if the hardware is
+   in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
+   (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
+   returned.  Or if the hardware is in round-upwards mode, then other bigger
+   values like say x == 2^128 will see x+0.5 rounding up to the next higher
+   representable value, 2^128+2^76 (or whatever), again incorrect.
+
+   These bad roundings of x+0.5 are avoided by testing at the start whether
+   x is already an integer.  If it is then clearly that's the desired result
+   already.  And if it's not then the exponent must be small enough to allow
+   an 0.5 to be represented, and hence added without a bad rounding.  */
+
 double
 scm_round (double x)
 {
-  double plus_half = x + 0.5;
-  double result = floor (plus_half);
+  double plus_half, result;
+
+  if (x == floor (x))
+    return x;
+
+  plus_half = x + 0.5;
+  result = floor (plus_half);
   /* Adjust so that the scm_round is towards even.  */
   return ((plus_half == result && plus_half / 2 != floor (plus_half / 2))
          ? result - 1

/* Copyright (C) 2004 Free Software Foundation, Inc.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */

#include "config.h"

#include <assert.h>
#include <math.h>
#include <stdio.h>

#if HAVE_FENV_H
#include <fenv.h>
#endif

#include "libguile.h"


#define numberof(x)  (sizeof (x) / sizeof ((x)[0]))

static void
test_scm_round ()
{
  /* FE constants are defined only where supported, in particular for
     instance some ARM systems have been seen with only a couple of modes */
  static const int modes[] = {
    0,
#ifdef FE_TONEAREST
    FE_TONEAREST,
#endif
#ifdef FE_UPWARD
    FE_UPWARD,
#endif
#ifdef FE_DOWNWARD
    FE_DOWNWARD,
#endif
#ifdef FE_TOWARDZERO
    FE_TOWARDZERO,
#endif
  };

  double  x, want;
  int  i;

  for (i = 0; i < numberof (modes); i++)
    {
      /* First iteration is the default rounding mode, ie. no call to
         fesetround.  Subsequent iterations are the FE modes from the
         table.  */
      if (i != 0)
        {
#if HAVE_FESETROUND
          fesetround (modes[i]);
#endif
        }

      assert (scm_round (0.0) == 0.0);
      assert (scm_round (1.0) == 1.0);
      assert (scm_round (-1.0) == -1.0);

      assert (scm_round (0.5) == 0.0);
      assert (scm_round (1.5) == 2.0);
      assert (scm_round (-1.5) == -2.0);
      assert (scm_round (2.5) == 2.0);
      assert (scm_round (-2.5) == -2.0);
      assert (scm_round (3.5) == 4.0);
      assert (scm_round (-3.5) == -4.0);

      /* 2^(DBL_MANT_DIG-1)-1+0.5 */
      x = ldexp (1.0, DBL_MANT_DIG - 1) - 1.0 + 0.5;
      want = ldexp (1.0, DBL_MANT_DIG - 1);
      assert (scm_round (x) == want);

      /* -(2^(DBL_MANT_DIG-1)-1+0.5) */
      x = - (ldexp (1.0, DBL_MANT_DIG - 1) - 1.0 + 0.5);
      want = - ldexp (1.0, DBL_MANT_DIG - 1);
      assert (scm_round (x) == want);

      /* 2^DBL_MANT_DIG-1
         In the past scm_round had incorrectly incremented this value, due
         to the way that x+0.5 would round upwards (in the usual default
         nearest-even mode on most systems).  */
      x = ldexp (1.0, DBL_MANT_DIG) - 1.0;
      assert (x == floor (x));      /* should be an integer already */
      assert (scm_round (x) == x);  /* scm_round should return it unchanged */

      /* -(2^DBL_MANT_DIG-1) */
      x = - (ldexp (1.0, DBL_MANT_DIG) - 1.0);
      assert (x == floor (x));      /* should be an integer already */
      assert (scm_round (x) == x);  /* scm_round should return it unchanged */

      /* 2^64 */
      x = ldexp (1.0, 64);
      assert (scm_round (x) == x);

      /* -2^64
         In the past scm_round had incorrectely incremented this value in
         any mode except FE_NEAREST, due to x+0.5 round either up or down to
         the next representable value (an integer).  */
      x = - ldexp (1.0, 64);
      assert (scm_round (x) == x);
    }
}

int
main (int argc, char *argv[])
{
  scm_init_guile();
  test_scm_round ();
  return 0;
}