update copyright file for recent contributions
[musl] / 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                 return 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