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/*
13 * We assume that a long double has a 15-bit exponent.  On systems
14 * where long double is the same as double, nexttoward() is an alias
15 * for nextafter(), so we don't use this routine.
16 */
17
18#include <float.h>
19
20#include "fpmath.h"
21#include "math.h"
22#include "math_private.h"
23
24#if LDBL_MAX_EXP != 0x4000
25#error "Unsupported long double format"
26#endif
27
28double
29nexttoward(double x, long double y)
30{
31	union IEEEl2bits uy;
32	volatile double t;
33	int32_t hx,ix;
34	u_int32_t lx;
35
36	EXTRACT_WORDS(hx,lx,x);
37	ix = hx&0x7fffffff;		/* |x| */
38	uy.e = y;
39
40	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||
41	    (uy.bits.exp == 0x7fff &&
42	     ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
43	   return x+y;	/* x or y is nan */
44	if(x==y) return (double)y;		/* x=y, return y */
45	if(x==0.0) {
46	    INSERT_WORDS(x,uy.bits.sign<<31,1);	/* return +-minsubnormal */
47	    t = x*x;
48	    if(t==x) return t; else return x;	/* raise underflow flag */
49	}
50	if(hx>0.0 ^ x < y) {			/* x -= ulp */
51	    if(lx==0) hx -= 1;
52	    lx -= 1;
53	} else {				/* x += ulp */
54	    lx += 1;
55	    if(lx==0) hx += 1;
56	}
57	ix = hx&0x7ff00000;
58	if(ix>=0x7ff00000) return x+x;	/* overflow  */
59	if(ix<0x00100000) {		/* underflow */
60	    t = x*x;
61	    if(t!=x) {		/* raise underflow flag */
62	        INSERT_WORDS(x,hx,lx);
63		return x;
64	    }
65	}
66	INSERT_WORDS(x,hx,lx);
67	return x;
68}
69