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 28 double 29 nexttoward(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