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
13static char rcsid[] = "$FreeBSD: src/lib/msun/src/s_nexttowardf.c,v 1.1 2005/03/07 04:57:38 das Exp $";
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
24float
25nexttowardf(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