pintos-os.org Git - pspp/blob - lib/round.c

   1 /* Round toward nearest, breaking ties away from zero.
   2    Copyright (C) 2007, 2010 Free Software Foundation, Inc.
   3
   4    This program is free software; you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 2, or (at your option)
   7    any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License along
  15    with this program; if not, write to the Free Software Foundation,
  16    Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
  17
  18 /* Written by Ben Pfaff <blp@gnu.org>, 2007.
  19    Based heavily on code by Bruno Haible. */
  20
  21 #include <config.h>
  22
  23 #include <float.h>
  24 #include <math.h>
  25
  26 #ifdef USE_LONG_DOUBLE
  27 # define ROUND roundl
  28 # define FLOOR floorl
  29 # define CEIL ceill
  30 # define DOUBLE long double
  31 # define MANT_DIG LDBL_MANT_DIG
  32 # define MIN LDBL_MIN
  33 # define L_(literal) literal##L
  34 # define HAVE_FLOOR_AND_CEIL HAVE_FLOORL_AND_CEILL
  35 #elif ! defined USE_FLOAT
  36 # define ROUND round
  37 # define FLOOR floor
  38 # define CEIL ceil
  39 # define DOUBLE double
  40 # define MANT_DIG DBL_MANT_DIG
  41 # define MIN DBL_MIN
  42 # define L_(literal) literal
  43 # define HAVE_FLOOR_AND_CEIL 1
  44 #else /* defined USE_FLOAT */
  45 # define ROUND roundf
  46 # define FLOOR floorf
  47 # define CEIL ceilf
  48 # define DOUBLE float
  49 # define MANT_DIG FLT_MANT_DIG
  50 # define MIN FLT_MIN
  51 # define L_(literal) literal##f
  52 # define HAVE_FLOOR_AND_CEIL HAVE_FLOORF_AND_CEILF
  53 #endif
  54
  55 /* -0.0.  See minus-zero.h.  */
  56 #if defined __hpux || defined __sgi || defined __ICC
  57 # define MINUS_ZERO (-MIN * MIN)
  58 #else
  59 # define MINUS_ZERO L_(-0.0)
  60 #endif
  61
  62 /* If we're being included from test-round2[f].c, it already defined names for
  63    our round implementations.  Otherwise, pick the preferred implementation for
  64    this machine. */
  65 #if !defined FLOOR_BASED_ROUND && !defined FLOOR_FREE_ROUND
  66 # if HAVE_FLOOR_AND_CEIL
  67 #  define FLOOR_BASED_ROUND ROUND
  68 # else
  69 #  define FLOOR_FREE_ROUND ROUND
  70 # endif
  71 #endif
  72
  73 #ifdef FLOOR_BASED_ROUND
  74 /* An implementation of the C99 round function based on floor and ceil.  We use
  75    this when floor and ceil are available, on the assumption that they are
  76    faster than the open-coded versions below. */
  77 DOUBLE
  78 FLOOR_BASED_ROUND (DOUBLE x)
  79 {
  80   if (x >= L_(0.0))
  81     {
  82       DOUBLE y = FLOOR (x);
  83       if (x - y >= L_(0.5))
  84         y += L_(1.0);
  85       return y;
  86     }
  87   else
  88     {
  89       DOUBLE y = CEIL (x);
  90       if (y - x >= L_(0.5))
  91         y -= L_(1.0);
  92       return y;
  93     }
  94 }
  95 #endif /* FLOOR_BASED_ROUND */
  96
  97 #ifdef FLOOR_FREE_ROUND
  98 /* An implementation of the C99 round function without floor or ceil.
  99    We use this when floor or ceil is missing. */
 100 DOUBLE
 101 FLOOR_FREE_ROUND (DOUBLE x)
 102 {
 103   /* 2^(MANT_DIG-1).  */
 104   static const DOUBLE TWO_MANT_DIG =
 105     /* Assume MANT_DIG <= 5 * 31.
 106        Use the identity
 107        n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5).  */
 108     (DOUBLE) (1U << ((MANT_DIG - 1) / 5))
 109     * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5))
 110     * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5))
 111     * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5))
 112     * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5));
 113
 114   /* The use of 'volatile' guarantees that excess precision bits are dropped at
 115      each addition step and before the following comparison at the caller's
 116      site.  It is necessary on x86 systems where double-floats are not IEEE
 117      compliant by default, to avoid that the results become platform and
 118      compiler option dependent.  'volatile' is a portable alternative to gcc's
 119      -ffloat-store option.  */
 120   volatile DOUBLE y = x;
 121   volatile DOUBLE z = y;
 122
 123   if (z > L_(0.0))
 124     {
 125       /* Avoid rounding error for x = 0.5 - 2^(-MANT_DIG-1).  */
 126       if (z < L_(0.5))
 127         z = L_(0.0);
 128       /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1.  */
 129       else if (z < TWO_MANT_DIG)
 130         {
 131           /* Add 0.5 to the absolute value.  */
 132           y = z += L_(0.5);
 133           /* Round to the next integer (nearest or up or down, doesn't
 134              matter).  */
 135           z += TWO_MANT_DIG;
 136           z -= TWO_MANT_DIG;
 137           /* Enforce rounding down.  */
 138           if (z > y)
 139             z -= L_(1.0);
 140         }
 141     }
 142   else if (z < L_(0.0))
 143     {
 144       /* Avoid rounding error for x = -(0.5 - 2^(-MANT_DIG-1)).  */
 145       if (z > - L_(0.5))
 146         z = MINUS_ZERO;
 147       /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1.  */
 148       else if (z > -TWO_MANT_DIG)
 149         {
 150           /* Add 0.5 to the absolute value.  */
 151           y = z -= L_(0.5);
 152           /* Round to the next integer (nearest or up or down, doesn't
 153              matter).  */
 154           z -= TWO_MANT_DIG;
 155           z += TWO_MANT_DIG;
 156           /* Enforce rounding up.  */
 157           if (z < y)
 158             z += L_(1.0);
 159         }
 160     }
 161   return z;
 162 }
 163 #endif /* FLOOR_FREE_ROUND */
 164