1 /* $NetBSD: s_nextafterl.c,v 1.2 2010/09/17 20:39:39 christos Exp $ */ 2 3 /* @(#)s_nextafter.c 5.1 93/09/24 */ 4 /* 5 * ==================================================== 6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 7 * 8 * Developed at SunPro, a Sun Microsystems, Inc. business. 9 * Permission to use, copy, modify, and distribute this 10 * software is freely granted, provided that this notice 11 * is preserved. 12 * ==================================================== 13 */ 14 15 #include <sys/cdefs.h> 16 __RCSID("$NetBSD: s_nextafterl.c,v 1.2 2010/09/17 20:39:39 christos Exp $"); 17 18 #include <float.h> 19 #include <math.h> 20 #include <machine/ieee.h> 21 22 #ifdef EXT_EXP_INFNAN 23 #if LDBL_MAX_EXP != 0x4000 24 #error "Unsupported long double format" 25 #endif 26 27 /* 28 * IEEE functions 29 * nextafterl(x,y) 30 * return the next machine floating-point number of x in the 31 * direction toward y. 32 * Special cases: 33 * If x == y, y shall be returned 34 * If x or y is NaN, a NaN shall be returned 35 */ 36 long double 37 nextafterl(long double x, long double y) 38 { 39 volatile long double t; 40 union ieee_ext_u ux, uy; 41 42 ux.extu_ld = x; 43 uy.extu_ld = y; 44 45 if ((ux.extu_exp == EXT_EXP_NAN && 46 ((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) || 47 (uy.extu_exp == EXT_EXP_NAN && 48 ((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0)) 49 return x+y; /* x or y is nan */ 50 51 if (x == y) return y; /* x=y, return y */ 52 53 if (x == 0.0) { 54 ux.extu_frach = 0; /* return +-minsubnormal */ 55 ux.extu_fracl = 1; 56 ux.extu_sign = uy.extu_sign; 57 t = ux.extu_ld * ux.extu_ld; 58 if (t == ux.extu_ld) 59 return t; 60 else 61 return ux.extu_ld; /* raise underflow flag */ 62 } 63 64 if ((x>0.0) ^ (x<y)) { /* x -= ulp */ 65 if (ux.extu_fracl == 0) { 66 if ((ux.extu_frach & ~LDBL_NBIT) == 0) 67 ux.extu_exp -= 1; 68 ux.extu_frach = (ux.extu_frach - 1) | 69 (ux.extu_frach & LDBL_NBIT); 70 } 71 ux.extu_fracl -= 1; 72 } else { /* x += ulp */ 73 ux.extu_fracl += 1; 74 if (ux.extu_fracl == 0) { 75 ux.extu_frach = (ux.extu_frach + 1) | 76 (ux.extu_frach & LDBL_NBIT); 77 if ((ux.extu_frach & ~LDBL_NBIT) == 0) 78 ux.extu_exp += 1; 79 } 80 } 81 82 if (ux.extu_exp == EXT_EXP_INF) 83 return x+x; /* overflow */ 84 85 if (ux.extu_exp == 0) { /* underflow */ 86 mask_nbit_l(ux); 87 t = ux.extu_ld * ux.extu_ld; 88 if (t != ux.extu_ld) /* raise underflow flag */ 89 return ux.extu_ld; 90 } 91 92 return ux.extu_ld; 93 } 94 #endif 95