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 * Copyright (c) 2002, 2003 David Schultz <das@FreeBSD.ORG> 12 */ 13 14 /* IEEE functions 15 * nextafterl(x,y) 16 * return the next machine floating-point number of x in the 17 * direction toward y. 18 * Special cases: 19 */ 20 21 #include <math.h> 22 #include <sys/cdefs.h> 23 24 #include "math_private.h" 25 26 union IEEEl2bits { 27 long double e; 28 struct { 29 unsigned int manl :32; 30 unsigned int manh :32; 31 unsigned int exp :15; 32 unsigned int sign :1; 33 unsigned int junkl :16; 34 unsigned int junkh :32; 35 } bits; 36 struct { 37 unsigned long man :64; 38 unsigned int expsign :16; 39 unsigned long junk :48; 40 } xbits; 41 }; 42 43 #define LDBL_NBIT 0x80000000 44 #define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT) 45 46 long double 47 nextafterl(long double x, long double y) 48 { 49 volatile long double t; 50 union IEEEl2bits ux, uy; 51 52 ux.e = x; 53 uy.e = y; 54 55 if ((ux.bits.exp == 0x7fff && 56 ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl) != 0) || 57 (uy.bits.exp == 0x7fff && 58 ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0)) 59 return x+y; /* x or y is nan */ 60 if(x==y) return y; /* x=y, return y */ 61 if(x==0.0) { 62 ux.bits.manh = 0; /* return +-minsubnormal */ 63 ux.bits.manl = 1; 64 ux.bits.sign = uy.bits.sign; 65 t = ux.e*ux.e; 66 if(t==ux.e) return t; else return ux.e; /* raise underflow flag */ 67 } 68 if((x>0.0) ^ (x<y)) { /* x -= ulp */ 69 if(ux.bits.manl==0) { 70 if ((ux.bits.manh&~LDBL_NBIT)==0) 71 ux.bits.exp -= 1; 72 ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT); 73 } 74 ux.bits.manl -= 1; 75 } else { /* x += ulp */ 76 ux.bits.manl += 1; 77 if(ux.bits.manl==0) { 78 ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT); 79 if ((ux.bits.manh&~LDBL_NBIT)==0) 80 ux.bits.exp += 1; 81 } 82 } 83 if(ux.bits.exp==0x7fff) return x+x; /* overflow */ 84 if(ux.bits.exp==0) { /* underflow */ 85 mask_nbit_l(ux); 86 t = ux.e * ux.e; 87 if(t!=ux.e) /* raise underflow flag */ 88 return ux.e; 89 } 90 return ux.e; 91 } 92 93 __strong_reference(nextafterl, nexttowardl); 94