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