1 /* s_nextafterl.c -- long double version of s_nextafter.c.
2 */
3
4 /*
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
7 *
8 * Developed at SunPro, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
11 * is preserved.
12 * ====================================================
13 */
14
15 #if defined(LIBM_SCCS) && !defined(lint)
16 static char rcsid[] = "$NetBSD: $";
17 #endif
18
19 /* IEEE functions
20 * nextafterl(x,y)
21 * return the next machine floating-point number of x in the
22 * direction toward y.
23 * Special cases:
24 */
25
26 #include <errno.h>
27 #include <float.h>
28 #include <math.h>
29 #include <math-barriers.h>
30 #include <math_private.h>
31 #include <math_ldbl_opt.h>
32
33 long double __nextafterl(long double x, long double y)
34 {
35 int64_t hx, hy, ihx, ihy, lx;
36 double xhi, xlo, yhi;
37
38 ldbl_unpack (x, &xhi, &xlo);
39 EXTRACT_WORDS64 (hx, xhi);
40 EXTRACT_WORDS64 (lx, xlo);
41 yhi = ldbl_high (y);
42 EXTRACT_WORDS64 (hy, yhi);
43 ihx = hx&0x7fffffffffffffffLL; /* |hx| */
44 ihy = hy&0x7fffffffffffffffLL; /* |hy| */
45
46 if((ihx>0x7ff0000000000000LL) || /* x is nan */
47 (ihy>0x7ff0000000000000LL)) /* y is nan */
48 return x+y; /* signal the nan */
49 if(x==y)
50 return y; /* x=y, return y */
51 if(ihx == 0) { /* x == 0 */
52 long double u; /* return +-minsubnormal */
53 hy = (hy & 0x8000000000000000ULL) | 1;
54 INSERT_WORDS64 (yhi, hy);
55 x = yhi;
56 u = math_opt_barrier (x);
57 u = u * u;
58 math_force_eval (u); /* raise underflow flag */
59 return x;
60 }
61
62 long double u;
63 if(x > y) { /* x > y, x -= ulp */
64 /* This isn't the largest magnitude correctly rounded
65 long double as you can see from the lowest mantissa
66 bit being zero. It is however the largest magnitude
67 long double with a 106 bit mantissa, and nextafterl
68 is insane with variable precision. So to make
69 nextafterl sane we assume 106 bit precision. */
70 if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL)) {
71 u = x+x; /* overflow, return -inf */
72 math_force_eval (u);
73 __set_errno (ERANGE);
74 return y;
75 }
76 if (hx >= 0x7ff0000000000000LL) {
77 u = 0x1.fffffffffffff7ffffffffffff8p+1023L;
78 return u;
79 }
80 if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */
81 u = math_opt_barrier (x);
82 x -= LDBL_TRUE_MIN;
83 if (ihx < 0x0360000000000000LL
84 || (hx > 0 && lx <= 0)
85 || (hx < 0 && lx > 1)) {
86 u = u * u;
87 math_force_eval (u); /* raise underflow flag */
88 __set_errno (ERANGE);
89 }
90 /* Avoid returning -0 in FE_DOWNWARD mode. */
91 if (x == 0.0L)
92 return 0.0L;
93 return x;
94 }
95 /* If the high double is an exact power of two and the low
96 double is the opposite sign, then 1ulp is one less than
97 what we might determine from the high double. Similarly
98 if X is an exact power of two, and positive, because
99 making it a little smaller will result in the exponent
100 decreasing by one and normalisation of the mantissa. */
101 if ((hx & 0x000fffffffffffffLL) == 0
102 && ((lx != 0 && (hx ^ lx) < 0)
103 || (lx == 0 && hx >= 0)))
104 ihx -= 1LL << 52;
105 if (ihx < (106LL << 52)) { /* ulp will denormal */
106 INSERT_WORDS64 (yhi, ihx & (0x7ffLL<<52));
107 u = yhi * 0x1p-105;
108 } else {
109 INSERT_WORDS64 (yhi, (ihx & (0x7ffLL<<52))-(105LL<<52));
110 u = yhi;
111 }
112 return x - u;
113 } else { /* x < y, x += ulp */
114 if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL)) {
115 u = x+x; /* overflow, return +inf */
116 math_force_eval (u);
117 __set_errno (ERANGE);
118 return y;
119 }
120 if ((uint64_t) hx >= 0xfff0000000000000ULL) {
121 u = -0x1.fffffffffffff7ffffffffffff8p+1023L;
122 return u;
123 }
124 if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */
125 u = math_opt_barrier (x);
126 x += LDBL_TRUE_MIN;
127 if (ihx < 0x0360000000000000LL
128 || (hx > 0 && lx < 0 && lx != 0x8000000000000001LL)
129 || (hx < 0 && lx >= 0)) {
130 u = u * u;
131 math_force_eval (u); /* raise underflow flag */
132 __set_errno (ERANGE);
133 }
134 if (x == 0.0L) /* handle negative LDBL_TRUE_MIN case */
135 x = -0.0L;
136 return x;
137 }
138 /* If the high double is an exact power of two and the low
139 double is the opposite sign, then 1ulp is one less than
140 what we might determine from the high double. Similarly
141 if X is an exact power of two, and negative, because
142 making it a little larger will result in the exponent
143 decreasing by one and normalisation of the mantissa. */
144 if ((hx & 0x000fffffffffffffLL) == 0
145 && ((lx != 0 && (hx ^ lx) < 0)
146 || (lx == 0 && hx < 0)))
147 ihx -= 1LL << 52;
148 if (ihx < (106LL << 52)) { /* ulp will denormal */
149 INSERT_WORDS64 (yhi, ihx & (0x7ffLL<<52));
150 u = yhi * 0x1p-105;
151 } else {
152 INSERT_WORDS64 (yhi, (ihx & (0x7ffLL<<52))-(105LL<<52));
153 u = yhi;
154 }
155 return x + u;
156 }
157 }
158 strong_alias (__nextafterl, __nexttowardl)
159 long_double_symbol (libm, __nextafterl, nextafterl);
160 long_double_symbol (libm, __nexttowardl, nexttowardl);