xref: /freebsd/lib/msun/src/s_nexttowardf.c (revision 1e413cf93298b5b97441a21d9a50fdcd0ee9945e)
1 /*
2  * ====================================================
3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4  *
5  * Developed at SunPro, a Sun Microsystems, Inc. business.
6  * Permission to use, copy, modify, and distribute this
7  * software is freely granted, provided that this notice
8  * is preserved.
9  * ====================================================
10  */
11 
12 #ifndef lint
13 static char rcsid[] = "$FreeBSD$";
14 #endif
15 
16 #include <float.h>
17 
18 #include "fpmath.h"
19 #include "math.h"
20 #include "math_private.h"
21 
22 #define	LDBL_INFNAN_EXP	(LDBL_MAX_EXP * 2 - 1)
23 
24 float
25 nexttowardf(float x, long double y)
26 {
27 	union IEEEl2bits uy;
28 	volatile float t;
29 	int32_t hx,ix;
30 
31 	GET_FLOAT_WORD(hx,x);
32 	ix = hx&0x7fffffff;		/* |x| */
33 	uy.e = y;
34 
35 	if((ix>0x7f800000) ||
36 	   (uy.bits.exp == LDBL_INFNAN_EXP &&
37 	    ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
38 	   return x+y;	/* x or y is nan */
39 	if(x==y) return (float)y;		/* x=y, return y */
40 	if(ix==0) {				/* x == 0 */
41 	    SET_FLOAT_WORD(x,(uy.bits.sign<<31)|1);/* return +-minsubnormal */
42 	    t = x*x;
43 	    if(t==x) return t; else return x;	/* raise underflow flag */
44 	}
45 	if(hx>=0 ^ x < y)			/* x -= ulp */
46 	    hx -= 1;
47 	else					/* x += ulp */
48 	    hx += 1;
49 	ix = hx&0x7f800000;
50 	if(ix>=0x7f800000) return x+x;	/* overflow  */
51 	if(ix<0x00800000) {		/* underflow */
52 	    t = x*x;
53 	    if(t!=x) {		/* raise underflow flag */
54 	        SET_FLOAT_WORD(y,hx);
55 		return y;
56 	    }
57 	}
58 	SET_FLOAT_WORD(x,hx);
59 	return x;
60 }
61