1 /* $NetBSD: s_nextafterl.c,v 1.3 2013/02/14 08:56:56 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.3 2013/02/14 08:56:56 matt 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 #ifdef LDBL_IMPLICIT_NBIT 28 #define LDBL_NBIT 0 29 #endif 30 31 /* 32 * IEEE functions 33 * nextafterl(x,y) 34 * return the next machine floating-point number of x in the 35 * direction toward y. 36 * Special cases: 37 * If x == y, y shall be returned 38 * If x or y is NaN, a NaN shall be returned 39 */ 40 long double 41 nextafterl(long double x, long double y) 42 { 43 volatile long double t; 44 union ieee_ext_u ux, uy; 45 46 ux.extu_ld = x; 47 uy.extu_ld = y; 48 49 if ((ux.extu_exp == EXT_EXP_NAN && 50 ((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) || 51 (uy.extu_exp == EXT_EXP_NAN && 52 ((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0)) 53 return x+y; /* x or y is nan */ 54 55 if (x == y) return y; /* x=y, return y */ 56 57 if (x == 0.0) { 58 ux.extu_frach = 0; /* return +-minsubnormal */ 59 ux.extu_fracl = 1; 60 ux.extu_sign = uy.extu_sign; 61 t = ux.extu_ld * ux.extu_ld; 62 if (t == ux.extu_ld) 63 return t; 64 else 65 return ux.extu_ld; /* raise underflow flag */ 66 } 67 68 if ((x>0.0) ^ (x<y)) { /* x -= ulp */ 69 if (ux.extu_fracl == 0) { 70 if ((ux.extu_frach & ~LDBL_NBIT) == 0) 71 ux.extu_exp -= 1; 72 ux.extu_frach = (ux.extu_frach - 1) | 73 (ux.extu_frach & LDBL_NBIT); 74 } 75 ux.extu_fracl -= 1; 76 } else { /* x += ulp */ 77 ux.extu_fracl += 1; 78 if (ux.extu_fracl == 0) { 79 ux.extu_frach = (ux.extu_frach + 1) | 80 (ux.extu_frach & LDBL_NBIT); 81 if ((ux.extu_frach & ~LDBL_NBIT) == 0) 82 ux.extu_exp += 1; 83 } 84 } 85 86 if (ux.extu_exp == EXT_EXP_INF) 87 return x+x; /* overflow */ 88 89 if (ux.extu_exp == 0) { /* underflow */ 90 #ifndef LDBL_IMPLICIT_NBIT 91 mask_nbit_l(ux); 92 #endif 93 t = ux.extu_ld * ux.extu_ld; 94 if (t != ux.extu_ld) /* raise underflow flag */ 95 return ux.extu_ld; 96 } 97 98 return ux.extu_ld; 99 } 100 #endif 101