granicus.if.org Git - postgresql/blob - src/port/rint.c

   1 /*-------------------------------------------------------------------------
   2  *
   3  * rint.c
   4  *        rint() implementation
   5  *
   6  * By Pedro Gimeno Fortea, donated to the public domain
   7  *
   8  * IDENTIFICATION
   9  *        src/port/rint.c
  10  *
  11  *-------------------------------------------------------------------------
  12  */
  13 #include "c.h"
  14
  15 #include <float.h>
  16 #include <math.h>
  17
  18 /*
  19  * Round to nearest integer, with halfway cases going to the nearest even.
  20  */
  21 double
  22 rint(double x)
  23 {
  24         double          x_orig;
  25         double          r;
  26
  27         /* Per POSIX, NaNs must be returned unchanged. */
  28         if (isnan(x))
  29                 return x;
  30
  31         if (x <= 0.0)
  32         {
  33                 /* Both positive and negative zero should be returned unchanged. */
  34                 if (x == 0.0)
  35                         return x;
  36
  37                 /*
  38                  * Subtracting 0.5 from a number very close to -0.5 can round to
  39                  * exactly -1.0, producing incorrect results, so we take the opposite
  40                  * approach: add 0.5 to the negative number, so that it goes closer to
  41                  * zero (or at most to +0.5, which is dealt with next), avoiding the
  42                  * precision issue.
  43                  */
  44                 x_orig = x;
  45                 x += 0.5;
  46
  47                 /*
  48                  * Be careful to return minus zero when input+0.5 >= 0, as that's what
  49                  * rint() should return with negative input.
  50                  */
  51                 if (x >= 0.0)
  52                         return -0.0;
  53
  54                 /*
  55                  * For very big numbers the input may have no decimals.  That case is
  56                  * detected by testing x+0.5 == x+1.0; if that happens, the input is
  57                  * returned unchanged.  This also covers the case of minus infinity.
  58                  */
  59                 if (x == x_orig + 1.0)
  60                         return x_orig;
  61
  62                 /* Otherwise produce a rounded estimate. */
  63                 r = floor(x);
  64
  65                 /*
  66                  * If the rounding did not produce exactly input+0.5 then we're done.
  67                  */
  68                 if (r != x)
  69                         return r;
  70
  71                 /*
  72                  * The original fractional part was exactly 0.5 (since
  73                  * floor(input+0.5) == input+0.5).  We need to round to nearest even.
  74                  * Dividing input+0.5 by 2, taking the floor and multiplying by 2
  75                  * yields the closest even number.  This part assumes that division by
  76                  * 2 is exact, which should be OK because underflow is impossible
  77                  * here: x is an integer.
  78                  */
  79                 return floor(x * 0.5) * 2.0;
  80         }
  81         else
  82         {
  83                 /*
  84                  * The positive case is similar but with signs inverted and using
  85                  * ceil() instead of floor().
  86                  */
  87                 x_orig = x;
  88                 x -= 0.5;
  89                 if (x <= 0.0)
  90                         return 0.0;
  91                 if (x == x_orig - 1.0)
  92                         return x_orig;
  93                 r = ceil(x);
  94                 if (r != x)
  95                         return r;
  96                 return ceil(x * 0.5) * 2.0;
  97         }
  98 }