a3b42c99d197e4fc925653e4535c1ba95c74cb8b
[musl] / 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                 return 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 }