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 #ifndef lint 14 static char rcsid[] = "$Id: s_nextafter.c,v 1.4 1994/03/03 17:04:43 jtc Exp $"; 15 #endif 16 17 /* IEEE functions 18 * nextafter(x,y) 19 * return the next machine floating-point number of x in the 20 * direction toward y. 21 * Special cases: 22 */ 23 24 #include <math.h> 25 #include <machine/endian.h> 26 27 #if BYTE_ORDER == LITTLE_ENDIAN 28 #define n0 1 29 #define n1 0 30 #else 31 #define n0 0 32 #define n1 1 33 #endif 34 35 #ifdef __STDC__ 36 double nextafter(double x, double y) 37 #else 38 double nextafter(x,y) 39 double x,y; 40 #endif 41 { 42 int hx,hy,ix,iy; 43 unsigned lx,ly; 44 45 hx = *( n0 + (int*)&x); /* high word of x */ 46 lx = *( n1 + (int*)&x); /* low word of x */ 47 hy = *( n0 + (int*)&y); /* high word of y */ 48 ly = *( n1 + (int*)&y); /* low word of y */ 49 ix = hx&0x7fffffff; /* |x| */ 50 iy = hy&0x7fffffff; /* |y| */ 51 52 if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */ 53 ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0)) /* y is nan */ 54 return x+y; 55 if(x==y) return x; /* x=y, return x */ 56 if((ix|lx)==0) { /* x == 0 */ 57 *(n0+(int*)&x) = hy&0x80000000; /* return +-minsubnormal */ 58 *(n1+(int*)&x) = 1; 59 y = x*x; 60 if(y==x) return y; else return x; /* raise underflow flag */ 61 } 62 if(hx>=0) { /* x > 0 */ 63 if(hx>hy||((hx==hy)&&(lx>ly))) { /* x > y, x -= ulp */ 64 if(lx==0) hx -= 1; 65 lx -= 1; 66 } else { /* x < y, x += ulp */ 67 lx += 1; 68 if(lx==0) hx += 1; 69 } 70 } else { /* x < 0 */ 71 if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */ 72 if(lx==0) hx -= 1; 73 lx -= 1; 74 } else { /* x > y, x += ulp */ 75 lx += 1; 76 if(lx==0) hx += 1; 77 } 78 } 79 hy = hx&0x7ff00000; 80 if(hy>=0x7ff00000) return x+x; /* overflow */ 81 if(hy<0x00100000) { /* underflow */ 82 y = x*x; 83 if(y!=x) { /* raise underflow flag */ 84 *(n0+(int*)&y) = hx; *(n1+(int*)&y) = lx; 85 return y; 86 } 87 } 88 *(n0+(int*)&x) = hx; *(n1+(int*)&x) = lx; 89 return x; 90 } 91