src/math/nexttoward.c

   1 #include "libm.h"
   2
   3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
   4 double nexttoward(double x, long double y)
   5 {
   6         return nextafter(x, y);
   7 }
   8 #else
   9 double nexttoward(double x, long double y)
  10 {
  11         union {double f; uint64_t i;} ux = {x};
  12         int e;
  13
  14         if (isnan(x) || isnan(y))
  15                 return x + y;
  16         if (x == y)
  17                 return y;
  18         if (x == 0) {
  19                 ux.i = 1;
  20                 if (signbit(y))
  21                         ux.i |= 1ULL<<63;
  22         } else if (x < y) {
  23                 if (signbit(x))
  24                         ux.i--;
  25                 else
  26                         ux.i++;
  27         } else {
  28                 if (signbit(x))
  29                         ux.i++;
  30                 else
  31                         ux.i--;
  32         }
  33         e = ux.i>>52 & 0x7ff;
  34         /* raise overflow if ux.f is infinite and x is finite */
  35         if (e == 0x7ff)
  36                 FORCE_EVAL(x+x);
  37         /* raise underflow if ux.f is subnormal or zero */
  38         if (e == 0)
  39                 FORCE_EVAL(x*x + ux.f*ux.f);
  40         return ux.f;
  41 }
  42 #endif