1 /* $NetBSD: s_nextafterl.c,v 1.6 2019/04/27 23:04:32 kamil Exp $ */
2
3 /* @(#)s_nextafter.c 5.1 93/09/24 */
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 #include <sys/cdefs.h>
16 __RCSID("$NetBSD: s_nextafterl.c,v 1.6 2019/04/27 23:04:32 kamil Exp $");
17
18 #include <float.h>
19 #include <math.h>
20 #include <machine/ieee.h>
21
22 #ifdef __HAVE_LONG_DOUBLE
23
24 #ifdef EXT_EXP_INFNAN
25 #if LDBL_MAX_EXP != 0x4000
26 #error "Unsupported long double format"
27 #endif
28
29 #ifdef LDBL_IMPLICIT_NBIT
30 #define LDBL_NBIT 0
31 #endif
32
__strong_alias(nexttowardl,nextafterl)33 __strong_alias(nexttowardl, nextafterl)
34
35 /*
36 * IEEE functions
37 * nextafterl(x,y)
38 * return the next machine floating-point number of x in the
39 * direction toward y.
40 * Special cases:
41 * If x == y, y shall be returned
42 * If x or y is NaN, a NaN shall be returned
43 */
44 long double
45 nextafterl(long double x, long double y)
46 {
47 volatile long double t;
48 union ieee_ext_u ux, uy;
49
50 ux.extu_ld = x;
51 uy.extu_ld = y;
52
53 if ((ux.extu_exp == EXT_EXP_INFNAN &&
54 ((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) ||
55 (uy.extu_exp == EXT_EXP_INFNAN &&
56 ((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0))
57 return x+y; /* x or y is nan */
58
59 if (x == y) return y; /* x=y, return y */
60
61 if (x == 0.0) {
62 ux.extu_frach = 0; /* return +-minsubnormal */
63 ux.extu_fracl = 1;
64 ux.extu_sign = uy.extu_sign;
65 t = ux.extu_ld * ux.extu_ld;
66 if (t == ux.extu_ld)
67 return t;
68 else
69 return ux.extu_ld; /* raise underflow flag */
70 }
71
72 if ((x>0.0) ^ (x<y)) { /* x -= ulp */
73 if (ux.extu_fracl == 0) {
74 if ((ux.extu_frach & ~LDBL_NBIT) == 0)
75 ux.extu_exp -= 1;
76 ux.extu_frach = (ux.extu_frach - 1) |
77 (ux.extu_frach & LDBL_NBIT);
78 }
79 ux.extu_fracl -= 1;
80 } else { /* x += ulp */
81 ux.extu_fracl += 1;
82 if (ux.extu_fracl == 0) {
83 ux.extu_frach = (ux.extu_frach + 1) |
84 (ux.extu_frach & LDBL_NBIT);
85 if ((ux.extu_frach & ~LDBL_NBIT) == 0)
86 ux.extu_exp += 1;
87 }
88 }
89
90 if (ux.extu_exp == EXT_EXP_INFNAN)
91 return x+x; /* overflow */
92
93 if (ux.extu_exp == 0) { /* underflow */
94 #ifndef LDBL_IMPLICIT_NBIT
95 mask_nbit_l(ux);
96 #endif
97 t = ux.extu_ld * ux.extu_ld;
98 if (t != ux.extu_ld) /* raise underflow flag */
99 return ux.extu_ld;
100 }
101
102 return ux.extu_ld;
103 }
104 #endif
105
106 #endif /* __HAVE_LONG_DOUBLE */
107