1 /* $NetBSD: s_nextafterl.c,v 1.5 2014/01/31 19:38:47 matt 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.5 2014/01/31 19:38:47 matt 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
33 /*
34 * IEEE functions
35 * nextafterl(x,y)
36 * return the next machine floating-point number of x in the
37 * direction toward y.
38 * Special cases:
39 * If x == y, y shall be returned
40 * If x or y is NaN, a NaN shall be returned
41 */
42 long double
nextafterl(long double x,long double y)43 nextafterl(long double x, long double y)
44 {
45 volatile long double t;
46 union ieee_ext_u ux, uy;
47
48 ux.extu_ld = x;
49 uy.extu_ld = y;
50
51 if ((ux.extu_exp == EXT_EXP_INFNAN &&
52 ((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) ||
53 (uy.extu_exp == EXT_EXP_INFNAN &&
54 ((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0))
55 return x+y; /* x or y is nan */
56
57 if (x == y) return y; /* x=y, return y */
58
59 if (x == 0.0) {
60 ux.extu_frach = 0; /* return +-minsubnormal */
61 ux.extu_fracl = 1;
62 ux.extu_sign = uy.extu_sign;
63 t = ux.extu_ld * ux.extu_ld;
64 if (t == ux.extu_ld)
65 return t;
66 else
67 return ux.extu_ld; /* raise underflow flag */
68 }
69
70 if ((x>0.0) ^ (x<y)) { /* x -= ulp */
71 if (ux.extu_fracl == 0) {
72 if ((ux.extu_frach & ~LDBL_NBIT) == 0)
73 ux.extu_exp -= 1;
74 ux.extu_frach = (ux.extu_frach - 1) |
75 (ux.extu_frach & LDBL_NBIT);
76 }
77 ux.extu_fracl -= 1;
78 } else { /* x += ulp */
79 ux.extu_fracl += 1;
80 if (ux.extu_fracl == 0) {
81 ux.extu_frach = (ux.extu_frach + 1) |
82 (ux.extu_frach & LDBL_NBIT);
83 if ((ux.extu_frach & ~LDBL_NBIT) == 0)
84 ux.extu_exp += 1;
85 }
86 }
87
88 if (ux.extu_exp == EXT_EXP_INFNAN)
89 return x+x; /* overflow */
90
91 if (ux.extu_exp == 0) { /* underflow */
92 #ifndef LDBL_IMPLICIT_NBIT
93 mask_nbit_l(ux);
94 #endif
95 t = ux.extu_ld * ux.extu_ld;
96 if (t != ux.extu_ld) /* raise underflow flag */
97 return ux.extu_ld;
98 }
99
100 return ux.extu_ld;
101 }
102 #endif
103
104 #endif /* __HAVE_LONG_DOUBLE */
105