105a0b428SJohn Marino /* @(#)s_nextafter.c 5.1 93/09/24 */
205a0b428SJohn Marino /*
305a0b428SJohn Marino * ====================================================
405a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
505a0b428SJohn Marino *
605a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
705a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
805a0b428SJohn Marino * software is freely granted, provided that this notice
905a0b428SJohn Marino * is preserved.
1005a0b428SJohn Marino * ====================================================
11*3d99825eSJohn Marino * Copyright (c) 2002, 2003 David Schultz <das@FreeBSD.ORG>
1205a0b428SJohn Marino */
1305a0b428SJohn Marino
1405a0b428SJohn Marino /* IEEE functions
1505a0b428SJohn Marino * nextafterl(x,y)
1605a0b428SJohn Marino * return the next machine floating-point number of x in the
1705a0b428SJohn Marino * direction toward y.
1805a0b428SJohn Marino * Special cases:
1905a0b428SJohn Marino */
2005a0b428SJohn Marino
2105a0b428SJohn Marino #include <math.h>
2274b7c7a8SJohn Marino #include <sys/cdefs.h>
2305a0b428SJohn Marino
2405a0b428SJohn Marino #include "math_private.h"
2505a0b428SJohn Marino
26*3d99825eSJohn Marino union IEEEl2bits {
27*3d99825eSJohn Marino long double e;
28*3d99825eSJohn Marino struct {
29*3d99825eSJohn Marino unsigned int manl :32;
30*3d99825eSJohn Marino unsigned int manh :32;
31*3d99825eSJohn Marino unsigned int exp :15;
32*3d99825eSJohn Marino unsigned int sign :1;
33*3d99825eSJohn Marino unsigned int junkl :16;
34*3d99825eSJohn Marino unsigned int junkh :32;
35*3d99825eSJohn Marino } bits;
36*3d99825eSJohn Marino struct {
37*3d99825eSJohn Marino unsigned long man :64;
38*3d99825eSJohn Marino unsigned int expsign :16;
39*3d99825eSJohn Marino unsigned long junk :48;
40*3d99825eSJohn Marino } xbits;
41*3d99825eSJohn Marino };
42*3d99825eSJohn Marino
43*3d99825eSJohn Marino #define LDBL_NBIT 0x80000000
44*3d99825eSJohn Marino #define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT)
45*3d99825eSJohn Marino
4605a0b428SJohn Marino long double
nextafterl(long double x,long double y)4705a0b428SJohn Marino nextafterl(long double x, long double y)
4805a0b428SJohn Marino {
49*3d99825eSJohn Marino volatile long double t;
50*3d99825eSJohn Marino union IEEEl2bits ux, uy;
5105a0b428SJohn Marino
52*3d99825eSJohn Marino ux.e = x;
53*3d99825eSJohn Marino uy.e = y;
5405a0b428SJohn Marino
55*3d99825eSJohn Marino if ((ux.bits.exp == 0x7fff &&
56*3d99825eSJohn Marino ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl) != 0) ||
57*3d99825eSJohn Marino (uy.bits.exp == 0x7fff &&
58*3d99825eSJohn Marino ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
59*3d99825eSJohn Marino return x+y; /* x or y is nan */
6005a0b428SJohn Marino if(x==y) return y; /* x=y, return y */
61*3d99825eSJohn Marino if(x==0.0) {
62*3d99825eSJohn Marino ux.bits.manh = 0; /* return +-minsubnormal */
63*3d99825eSJohn Marino ux.bits.manl = 1;
64*3d99825eSJohn Marino ux.bits.sign = uy.bits.sign;
65*3d99825eSJohn Marino t = ux.e*ux.e;
66*3d99825eSJohn Marino if(t==ux.e) return t; else return ux.e; /* raise underflow flag */
6705a0b428SJohn Marino }
68*3d99825eSJohn Marino if((x>0.0) ^ (x<y)) { /* x -= ulp */
69*3d99825eSJohn Marino if(ux.bits.manl==0) {
70*3d99825eSJohn Marino if ((ux.bits.manh&~LDBL_NBIT)==0)
71*3d99825eSJohn Marino ux.bits.exp -= 1;
72*3d99825eSJohn Marino ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT);
7305a0b428SJohn Marino }
74*3d99825eSJohn Marino ux.bits.manl -= 1;
75*3d99825eSJohn Marino } else { /* x += ulp */
76*3d99825eSJohn Marino ux.bits.manl += 1;
77*3d99825eSJohn Marino if(ux.bits.manl==0) {
78*3d99825eSJohn Marino ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT);
79*3d99825eSJohn Marino if ((ux.bits.manh&~LDBL_NBIT)==0)
80*3d99825eSJohn Marino ux.bits.exp += 1;
8105a0b428SJohn Marino }
8205a0b428SJohn Marino }
83*3d99825eSJohn Marino if(ux.bits.exp==0x7fff) return x+x; /* overflow */
84*3d99825eSJohn Marino if(ux.bits.exp==0) { /* underflow */
85*3d99825eSJohn Marino mask_nbit_l(ux);
86*3d99825eSJohn Marino t = ux.e * ux.e;
87*3d99825eSJohn Marino if(t!=ux.e) /* raise underflow flag */
88*3d99825eSJohn Marino return ux.e;
8905a0b428SJohn Marino }
90*3d99825eSJohn Marino return ux.e;
9105a0b428SJohn Marino }
9205a0b428SJohn Marino
9374b7c7a8SJohn Marino __strong_reference(nextafterl, nexttowardl);
94