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 /* IEEE functions 13 * nextafter(x,y) 14 * return the next machine floating-point number of x in the 15 * direction toward y. 16 * Special cases: 17 */ 18 19 #include <float.h> 20 21 #include "fpmath.h" 22 #include "math.h" 23 #include "math_private.h" 24 25 #if LDBL_MAX_EXP != 0x4000 26 #error "Unsupported long double format" 27 #endif 28 29 long double 30 nextafterl(long double x, long double y) 31 { 32 volatile long double t; 33 union IEEEl2bits ux, uy; 34 35 ux.e = x; 36 uy.e = y; 37 38 if ((ux.bits.exp == 0x7fff && 39 ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl) != 0) || 40 (uy.bits.exp == 0x7fff && 41 ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0)) 42 return x+y; /* x or y is nan */ 43 if(x==y) return y; /* x=y, return y */ 44 if(x==0.0) { 45 ux.bits.manh = 0; /* return +-minsubnormal */ 46 ux.bits.manl = 1; 47 ux.bits.sign = uy.bits.sign; 48 t = ux.e*ux.e; 49 if(t==ux.e) return t; else return ux.e; /* raise underflow flag */ 50 } 51 if(x>0.0 ^ x<y) { /* x -= ulp */ 52 if(ux.bits.manl==0) { 53 if ((ux.bits.manh&~LDBL_NBIT)==0) 54 ux.bits.exp -= 1; 55 ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT); 56 } 57 ux.bits.manl -= 1; 58 } else { /* x += ulp */ 59 ux.bits.manl += 1; 60 if(ux.bits.manl==0) { 61 ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT); 62 if ((ux.bits.manh&~LDBL_NBIT)==0) 63 ux.bits.exp += 1; 64 } 65 } 66 if(ux.bits.exp==0x7fff) return x+x; /* overflow */ 67 if(ux.bits.exp==0) { /* underflow */ 68 mask_nbit_l(ux); 69 t = ux.e * ux.e; 70 if(t!=ux.e) /* raise underflow flag */ 71 return ux.e; 72 } 73 return ux.e; 74 } 75 76 __strong_reference(nextafterl, nexttowardl); 77