1 /* @(#)s_nextafter.c 5.1 93/09/24 */
2 /*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 * Copyright (c) 2002, 2003 David Schultz <das@FreeBSD.ORG>
12 */
13
14 /* IEEE functions
15 * nextafterl(x,y)
16 * return the next machine floating-point number of x in the
17 * direction toward y.
18 * Special cases:
19 */
20
21 #include <math.h>
22 #include <sys/cdefs.h>
23
24 #include "math_private.h"
25
26 union IEEEl2bits {
27 long double e;
28 struct {
29 unsigned int manl :32;
30 unsigned int manh :32;
31 unsigned int exp :15;
32 unsigned int sign :1;
33 unsigned int junkl :16;
34 unsigned int junkh :32;
35 } bits;
36 struct {
37 unsigned long man :64;
38 unsigned int expsign :16;
39 unsigned long junk :48;
40 } xbits;
41 };
42
43 #define LDBL_NBIT 0x80000000
44 #define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT)
45
46 long double
nextafterl(long double x,long double y)47 nextafterl(long double x, long double y)
48 {
49 volatile long double t;
50 union IEEEl2bits ux, uy;
51
52 ux.e = x;
53 uy.e = y;
54
55 if ((ux.bits.exp == 0x7fff &&
56 ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl) != 0) ||
57 (uy.bits.exp == 0x7fff &&
58 ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
59 return x+y; /* x or y is nan */
60 if(x==y) return y; /* x=y, return y */
61 if(x==0.0) {
62 ux.bits.manh = 0; /* return +-minsubnormal */
63 ux.bits.manl = 1;
64 ux.bits.sign = uy.bits.sign;
65 t = ux.e*ux.e;
66 if(t==ux.e) return t; else return ux.e; /* raise underflow flag */
67 }
68 if((x>0.0) ^ (x<y)) { /* x -= ulp */
69 if(ux.bits.manl==0) {
70 if ((ux.bits.manh&~LDBL_NBIT)==0)
71 ux.bits.exp -= 1;
72 ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT);
73 }
74 ux.bits.manl -= 1;
75 } else { /* x += ulp */
76 ux.bits.manl += 1;
77 if(ux.bits.manl==0) {
78 ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT);
79 if ((ux.bits.manh&~LDBL_NBIT)==0)
80 ux.bits.exp += 1;
81 }
82 }
83 if(ux.bits.exp==0x7fff) return x+x; /* overflow */
84 if(ux.bits.exp==0) { /* underflow */
85 mask_nbit_l(ux);
86 t = ux.e * ux.e;
87 if(t!=ux.e) /* raise underflow flag */
88 return ux.e;
89 }
90 return ux.e;
91 }
92
93 __strong_reference(nextafterl, nexttowardl);
94