1 /* $NetBSD: s_nextafterl.c,v 1.5 2014/01/31 19:38:47 matt 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.5 2014/01/31 19:38:47 matt Exp $"); 17 18 #include <float.h> 19 #include <math.h> 20 #include <machine/ieee.h> 21 22 #ifdef __HAVE_LONG_DOUBLE 23 24 #ifdef EXT_EXP_INFNAN 25 #if LDBL_MAX_EXP != 0x4000 26 #error "Unsupported long double format" 27 #endif 28 29 #ifdef LDBL_IMPLICIT_NBIT 30 #define LDBL_NBIT 0 31 #endif 32 33 /* 34 * IEEE functions 35 * nextafterl(x,y) 36 * return the next machine floating-point number of x in the 37 * direction toward y. 38 * Special cases: 39 * If x == y, y shall be returned 40 * If x or y is NaN, a NaN shall be returned 41 */ 42 long double 43 nextafterl(long double x, long double y) 44 { 45 volatile long double t; 46 union ieee_ext_u ux, uy; 47 48 ux.extu_ld = x; 49 uy.extu_ld = y; 50 51 if ((ux.extu_exp == EXT_EXP_INFNAN && 52 ((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) || 53 (uy.extu_exp == EXT_EXP_INFNAN && 54 ((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0)) 55 return x+y; /* x or y is nan */ 56 57 if (x == y) return y; /* x=y, return y */ 58 59 if (x == 0.0) { 60 ux.extu_frach = 0; /* return +-minsubnormal */ 61 ux.extu_fracl = 1; 62 ux.extu_sign = uy.extu_sign; 63 t = ux.extu_ld * ux.extu_ld; 64 if (t == ux.extu_ld) 65 return t; 66 else 67 return ux.extu_ld; /* raise underflow flag */ 68 } 69 70 if ((x>0.0) ^ (x<y)) { /* x -= ulp */ 71 if (ux.extu_fracl == 0) { 72 if ((ux.extu_frach & ~LDBL_NBIT) == 0) 73 ux.extu_exp -= 1; 74 ux.extu_frach = (ux.extu_frach - 1) | 75 (ux.extu_frach & LDBL_NBIT); 76 } 77 ux.extu_fracl -= 1; 78 } else { /* x += ulp */ 79 ux.extu_fracl += 1; 80 if (ux.extu_fracl == 0) { 81 ux.extu_frach = (ux.extu_frach + 1) | 82 (ux.extu_frach & LDBL_NBIT); 83 if ((ux.extu_frach & ~LDBL_NBIT) == 0) 84 ux.extu_exp += 1; 85 } 86 } 87 88 if (ux.extu_exp == EXT_EXP_INFNAN) 89 return x+x; /* overflow */ 90 91 if (ux.extu_exp == 0) { /* underflow */ 92 #ifndef LDBL_IMPLICIT_NBIT 93 mask_nbit_l(ux); 94 #endif 95 t = ux.extu_ld * ux.extu_ld; 96 if (t != ux.extu_ld) /* raise underflow flag */ 97 return ux.extu_ld; 98 } 99 100 return ux.extu_ld; 101 } 102 #endif 103 104 #endif /* __HAVE_LONG_DOUBLE */ 105