s_nextafter.c revision 140685 
119304Speter/* @(#)s_nextafter.c 5.1 93/09/24 */
219304Speter/*
319304Speter * ====================================================
419304Speter * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
519304Speter *
619304Speter * Developed at SunPro, a Sun Microsystems, Inc. business.
719304Speter * Permission to use, copy, modify, and distribute this
819304Speter * software is freely granted, provided that this notice
919304Speter * is preserved.
1019304Speter * ====================================================
1119304Speter */
1219304Speter
1319304Speter#ifndef lint
1419304Speterstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/s_nextafter.c 140685 2005-01-23 22:56:08Z das $";
1519304Speter#endif
1619304Speter
1719304Speter/* IEEE functions
1819304Speter *	nextafter(x,y)
1919304Speter *	return the next machine floating-point number of x in the
2019304Speter *	direction toward y.
2119304Speter *   Special cases:
2219304Speter */
2319304Speter
2419304Speter#include "math.h"
2519304Speter#include "math_private.h"
2619304Speter
2719304Speterdouble
2819304Speternextafter(double x, double y)
2919304Speter{
3019304Speter	int32_t hx,hy,ix,iy;
3119304Speter	u_int32_t lx,ly;
3219304Speter
33	EXTRACT_WORDS(hx,lx,x);
34	EXTRACT_WORDS(hy,ly,y);
35	ix = hx&0x7fffffff;		/* |x| */
36	iy = hy&0x7fffffff;		/* |y| */
37
38	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */
39	   ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0))     /* y is nan */
40	   return x+y;
41	if(x==y) return y;		/* x=y, return y */
42	if((ix|lx)==0) {			/* x == 0 */
43	    INSERT_WORDS(x,hy&0x80000000,1);	/* return +-minsubnormal */
44	    y = x*x;
45	    if(y==x) return y; else return x;	/* raise underflow flag */
46	}
47	if(hx>=0) {				/* x > 0 */
48	    if(hx>hy||((hx==hy)&&(lx>ly))) {	/* x > y, x -= ulp */
49		if(lx==0) hx -= 1;
50		lx -= 1;
51	    } else {				/* x < y, x += ulp */
52		lx += 1;
53		if(lx==0) hx += 1;
54	    }
55	} else {				/* x < 0 */
56	    if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */
57		if(lx==0) hx -= 1;
58		lx -= 1;
59	    } else {				/* x > y, x += ulp */
60		lx += 1;
61		if(lx==0) hx += 1;
62	    }
63	}
64	hy = hx&0x7ff00000;
65	if(hy>=0x7ff00000) return x+x;	/* overflow  */
66	if(hy<0x00100000) {		/* underflow */
67	    y = x*x;
68	    if(y!=x) {		/* raise underflow flag */
69	        INSERT_WORDS(y,hx,lx);
70		return y;
71	    }
72	}
73	INSERT_WORDS(x,hx,lx);
74	return x;
75}
76