1 /* $NetBSD: s_nexttowardf.c,v 1.3 2013/02/09 23:14:44 christos Exp $ */ 2 3 /* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunPro, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 */ 13 14 #include <sys/cdefs.h> 15 #if 0 16 __FBSDID("$FreeBSD: src/lib/msun/src/s_nexttowardf.c,v 1.3 2011/02/10 07:38:38 das Exp $"); 17 #else 18 __RCSID("$NetBSD: s_nexttowardf.c,v 1.3 2013/02/09 23:14:44 christos Exp $"); 19 #endif 20 21 #include <string.h> 22 #include <float.h> 23 #include <machine/ieee.h> 24 25 #include "math.h" 26 #include "math_private.h" 27 28 #ifdef EXT_EXP_INFNAN 29 float 30 nexttowardf(float x, long double y) 31 { 32 volatile float t; 33 int32_t hx,ix; 34 union ieee_ext_u uy; 35 36 GET_FLOAT_WORD(hx,x); 37 ix = hx&0x7fffffff; /* |x| */ 38 39 memset(&uy, 0, sizeof(uy)); 40 uy.extu_ld = y; 41 uy.extu_ext.ext_frach &= ~0x80000000; 42 43 if((ix>0x7f800000) || 44 (uy.extu_ext.ext_exp == EXT_EXP_INFNAN && 45 (uy.extu_ext.ext_frach | uy.extu_ext.ext_fracl) != 0)) 46 return x+y; /* x or y is nan */ 47 if(x==y) return (float)y; /* x=y, return y */ 48 if(ix==0) { /* x == 0 */ 49 SET_FLOAT_WORD(x,(uy.extu_ext.ext_sign<<31)|1);/* return +-minsubnormal */ 50 t = x*x; 51 if(t==x) return t; else return x; /* raise underflow flag */ 52 } 53 if((hx >= 0) ^ (x < y)) /* x -= ulp */ 54 hx -= 1; 55 else /* x += ulp */ 56 hx += 1; 57 ix = hx&0x7f800000; 58 if(ix>=0x7f800000) return x+x; /* overflow */ 59 if(ix<0x00800000) { /* underflow */ 60 t = x*x; 61 if(t!=x) { /* raise underflow flag */ 62 SET_FLOAT_WORD(x,hx); 63 return x; 64 } 65 } 66 SET_FLOAT_WORD(x,hx); 67 return x; 68 } 69 #endif 70