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); |
161 | |