src/math/nextafter.c

   1 #include "libm.h"
   2
   3 #define SIGN ((uint64_t)1<<63)
   4
   5 double nextafter(double x, double y)
   6 {
   7         union dshape ux, uy;
   8         uint64_t ax, ay;
   9         int e;
  10
  11         if (isnan(x) || isnan(y))
  12                 return x + y;
  13         ux.value = x;
  14         uy.value = y;
  15         if (ux.bits == uy.bits)
  16                 return y;
  17         ax = ux.bits & ~SIGN;
  18         ay = uy.bits & ~SIGN;
  19         if (ax == 0) {
  20                 if (ay == 0)
  21                         return y;
  22                 ux.bits = (uy.bits & SIGN) | 1;
  23         } else if (ax > ay || ((ux.bits ^ uy.bits) & SIGN))
  24                 ux.bits--;
  25         else
  26                 ux.bits++;
  27         e = ux.bits >> 52 & 0x7ff;
  28         /* raise overflow if ux.value is infinite and x is finite */
  29         if (e == 0x7ff)
  30                 FORCE_EVAL(x+x);
  31         /* raise underflow if ux.value is subnormal or zero */
  32         if (e == 0)
  33                 FORCE_EVAL(x*x + ux.value*ux.value);
  34         return ux.value;
  35 }