nexttoward.c raw

   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
  43