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