Home | History | Annotate | Download | only in src
      1 /* @(#)s_nextafter.c 5.1 93/09/24 */
      2 /*
      3  * ====================================================
      4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      5  *
      6  * Developed at SunPro, a Sun Microsystems, Inc. business.
      7  * Permission to use, copy, modify, and distribute this
      8  * software is freely granted, provided that this notice
      9  * is preserved.
     10  * ====================================================
     11  */
     12 
     13 #include <sys/cdefs.h>
     14 __FBSDID("$FreeBSD$");
     15 
     16 /* IEEE functions
     17  *	nextafter(x,y)
     18  *	return the next machine floating-point number of x in the
     19  *	direction toward y.
     20  *   Special cases:
     21  */
     22 
     23 #include <float.h>
     24 
     25 #include "math.h"
     26 #include "math_private.h"
     27 
     28 double
     29 nextafter(double x, double y)
     30 {
     31 	volatile double t;
     32 	int32_t hx,hy,ix,iy;
     33 	u_int32_t lx,ly;
     34 
     35 	EXTRACT_WORDS(hx,lx,x);
     36 	EXTRACT_WORDS(hy,ly,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 y;		/* x=y, return y */
     44 	if((ix|lx)==0) {			/* x == 0 */
     45 	    INSERT_WORDS(x,hy&0x80000000,1);	/* return +-minsubnormal */
     46 	    t = x*x;
     47 	    if(t==x) return t; else return x;	/* raise underflow flag */
     48 	}
     49 	if(hx>=0) {				/* x > 0 */
     50 	    if(hx>hy||((hx==hy)&&(lx>ly))) {	/* x > y, x -= ulp */
     51 		if(lx==0) hx -= 1;
     52 		lx -= 1;
     53 	    } else {				/* x < y, x += ulp */
     54 		lx += 1;
     55 		if(lx==0) hx += 1;
     56 	    }
     57 	} else {				/* x < 0 */
     58 	    if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */
     59 		if(lx==0) hx -= 1;
     60 		lx -= 1;
     61 	    } else {				/* x > y, x += ulp */
     62 		lx += 1;
     63 		if(lx==0) hx += 1;
     64 	    }
     65 	}
     66 	hy = hx&0x7ff00000;
     67 	if(hy>=0x7ff00000) return x+x;	/* overflow  */
     68 	if(hy<0x00100000) {		/* underflow */
     69 	    t = x*x;
     70 	    if(t!=x) {		/* raise underflow flag */
     71 	        INSERT_WORDS(y,hx,lx);
     72 		return y;
     73 	    }
     74 	}
     75 	INSERT_WORDS(x,hx,lx);
     76 	return x;
     77 }
     78 
     79 #if (LDBL_MANT_DIG == 53)
     80 __weak_reference(nextafter, nexttoward);
     81 __weak_reference(nextafter, nexttowardl);
     82 __weak_reference(nextafter, nextafterl);
     83 #endif
     84