1 /* $NetBSD: s_nexttoward.c,v 1.3 2024/05/11 02:07:53 riastradh 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_nexttoward.c,v 1.3 2024/05/11 02:07:53 riastradh Exp $");
17
18 /*
19 * We assume that a long double has a 15-bit exponent. On systems
20 * where long double is the same as double, nexttoward() is an alias
21 * for nextafter(), so we don't use this routine.
22 */
23 #include <float.h>
24
25 #include <machine/ieee.h>
26 #include "math.h"
27 #include "math_private.h"
28
29 #if LDBL_MAX_EXP != 0x4000
30 #error "Unsupported long double format"
31 #endif
32
33 #ifdef LDBL_IMPLICIT_NBIT
34 #define LDBL_NBIT 0
35 #endif
36
37 /*
38 * The nexttoward() function is equivalent to nextafter() function,
39 * except that the second parameter shall have type long double and
40 * the functions shall return y converted to the type of the function
41 * if x equals y.
42 *
43 * Special cases: XXX
44 */
45 double
nexttoward(double x,long double y)46 nexttoward(double x, long double y)
47 {
48 union ieee_ext_u uy;
49 volatile double t;
50 int32_t hx, ix;
51 uint32_t lx;
52
53 EXTRACT_WORDS(hx, lx, x);
54 ix = hx & 0x7fffffff; /* |x| */
55 uy.extu_ld = y;
56
57 if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) ||
58 (uy.extu_exp == 0x7fff &&
59 ((uy.extu_frach & ~LDBL_NBIT) | uy.extu_fracl) != 0))
60 return x+y; /* x or y is nan */
61
62 if (x == y)
63 return (double)y; /* x=y, return y */
64
65 if (x == 0.0) {
66 INSERT_WORDS(x, uy.extu_sign<<31, 1); /* return +-minsubnormal */
67 t = x*x;
68 if (t == x)
69 return t;
70 else
71 return x; /* raise underflow flag */
72 }
73
74 if ((hx >= 0) ^ (x < y)) { /* x -= ulp */
75 if (lx == 0) hx -= 1;
76 lx -= 1;
77 } else { /* x += ulp */
78 lx += 1;
79 if (lx == 0) hx += 1;
80 }
81 ix = hx & 0x7ff00000;
82 if (ix >= 0x7ff00000) return x+x; /* overflow */
83 if (ix < 0x00100000) { /* underflow */
84 t = x*x;
85 if (t != x) { /* raise underflow flag */
86 INSERT_WORDS(y, hx, lx);
87 return y;
88 }
89 }
90 INSERT_WORDS(x, hx, lx);
91
92 return x;
93 }
94