1
2/* @(#)s_nextafter.c 1.3 95/01/18 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14/* IEEE functions
15 *	nextafter(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 "fdlibm.h"
22
23#ifdef __STDC__
24	double ieee_nextafter(double x, double y)
25#else
26	double ieee_nextafter(x,y)
27	double x,y;
28#endif
29{
30	int	hx,hy,ix,iy;
31	unsigned lx,ly;
32
33	hx = __HI(x);		/* high word of x */
34	lx = __LO(x);		/* low  word of x */
35	hy = __HI(y);		/* high word of y */
36	ly = __LO(y);		/* low  word of y */
37	ix = hx&0x7fffffff;		/* |x| */
38	iy = hy&0x7fffffff;		/* |y| */
39
40	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */
41	   ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0))     /* y is nan */
42	   return x+y;
43	if(x==y) return x;		/* x=y, return x */
44	if((ix|lx)==0) {			/* x == 0 */
45	    __HI(x) = hy&0x80000000;	/* return +-minsubnormal */
46	    __LO(x) = 1;
47	    y = x*x;
48	    if(y==x) return y; else return x;	/* raise underflow flag */
49	}
50	if(hx>=0) {				/* x > 0 */
51	    if(hx>hy||((hx==hy)&&(lx>ly))) {	/* x > y, x -= ulp */
52		if(lx==0) hx -= 1;
53		lx -= 1;
54	    } else {				/* x < y, x += ulp */
55		lx += 1;
56		if(lx==0) hx += 1;
57	    }
58	} else {				/* x < 0 */
59	    if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */
60		if(lx==0) hx -= 1;
61		lx -= 1;
62	    } else {				/* x > y, x += ulp */
63		lx += 1;
64		if(lx==0) hx += 1;
65	    }
66	}
67	hy = hx&0x7ff00000;
68	if(hy>=0x7ff00000) return x+x;	/* overflow  */
69	if(hy<0x00100000) {		/* underflow */
70	    y = x*x;
71	    if(y!=x) {		/* raise underflow flag */
72		__HI(y) = hx; __LO(y) = lx;
73		return y;
74	    }
75	}
76	__HI(x) = hx; __LO(x) = lx;
77	return x;
78}
79