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 #define SIGN ((uint64_t)1<<63)
  10
  11 double nexttoward(double x, long double y)
  12 {
  13         union dshape ux;
  14         int e;
  15
  16         if (isnan(x) || isnan(y))
  17                 return x + y;
  18         if (x == y)
  19                 return y;
  20         ux.value = x;
  21         if (x == 0) {
  22                 ux.bits = 1;
  23                 if (signbit(y))
  24                         ux.bits |= SIGN;
  25         } else if (x < y) {
  26                 if (signbit(x))
  27                         ux.bits--;
  28                 else
  29                         ux.bits++;
  30         } else {
  31                 if (signbit(x))
  32                         ux.bits++;
  33                 else
  34                         ux.bits--;
  35         }
  36         e = ux.bits>>52 & 0x7ff;
  37         /* raise overflow if ux.value is infinite and x is finite */
  38         if (e == 0x7ff)
  39                 FORCE_EVAL(x+x);
  40         /* raise underflow if ux.value is subnormal or zero */
  41         if (e == 0)
  42                 FORCE_EVAL(x*x + ux.value*ux.value);
  43         return ux.value;
  44 }
  45 #endif