missing/nextafter.c

   1 #include "ruby/missing.h"
   2
   3 #include <math.h>
   4 #include <float.h>
   5
   6 /* This function doesn't set errno.  It should on POSIX, though. */
   7
   8 double
   9 nextafter(double x, double y)
  10 {
  11     double x1, x2, d;
  12     int e;
  13
  14     if (isnan(x))
  15         return x;
  16     if (isnan(y))
  17         return y;
  18
  19     if (x == y)
  20         return y;
  21
  22     if (x == 0) {
  23         /* the minimum "subnormal" float */
  24         x1 = ldexp(0.5, DBL_MIN_EXP - DBL_MANT_DIG + 1);
  25         if (x1 == 0)
  26             x1 = DBL_MIN; /* the minimum "normal" float */
  27         if (0 < y)
  28             return x1;
  29         else
  30             return -x1;
  31     }
  32
  33     if (x < 0) {
  34         if (isinf(x))
  35             return -DBL_MAX;
  36         if (x == -DBL_MAX && y < 0 && isinf(y))
  37             return y;
  38     }
  39     else {
  40         if (isinf(x))
  41             return DBL_MAX;
  42         if (x == DBL_MAX && 0 < y && isinf(y))
  43             return y;
  44     }
  45
  46     x1 = frexp(x, &e);
  47
  48     if (x < y) {
  49         d = DBL_EPSILON/2;
  50         if (x1 == -0.5) {
  51             x1 *= 2;
  52             e--;
  53         }
  54     }
  55     else {
  56         d = -DBL_EPSILON/2;
  57         if (x1 == 0.5) {
  58             x1 *= 2;
  59             e--;
  60         }
  61     }
  62
  63     if (e < DBL_MIN_EXP) {
  64         d = ldexp(d, DBL_MIN_EXP-e);
  65     }
  66
  67     x2 = x1 + d;
  68
  69     if (x2 == 0.0) {
  70         if (x1 < 0)
  71             return -0.0;
  72         else
  73             return +0.0;
  74     }
  75
  76     return ldexp(x2, e);
  77 }