1 /* $OpenBSD: s_logbl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $ */ 2 /* 3 * From: @(#)s_ilogb.c 5.1 93/09/24 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/types.h> 15 #include <machine/ieee.h> 16 #include <float.h> 17 #include <limits.h> 18 #include <math.h> 19 20 long double 21 logbl(long double x) 22 { 23 union { 24 long double e; 25 struct ieee_ext bits; 26 } u; 27 unsigned long m; 28 int b; 29 30 u.e = x; 31 if (u.bits.ext_exp == 0) { 32 if ((u.bits.ext_fracl 33 #ifdef EXT_FRACLMBITS 34 | u.bits.ext_fraclm 35 #endif /* EXT_FRACLMBITS */ 36 #ifdef EXT_FRACHMBITS 37 | u.bits.ext_frachm 38 #endif /* EXT_FRACHMBITS */ 39 | u.bits.ext_frach) == 0) { /* x == 0 */ 40 u.bits.ext_sign = 1; 41 return (1.0L / u.e); 42 } 43 /* denormalized */ 44 if (u.bits.ext_frach == 0 45 #ifdef EXT_FRACHMBITS 46 && u.bits.ext_frachm == 0 47 #endif 48 ) { 49 m = 1lu << (EXT_FRACLBITS - 1); 50 for (b = EXT_FRACHBITS; !(u.bits.ext_fracl & m); m >>= 1) 51 b++; 52 #if defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) 53 m = 1lu << (EXT_FRACLMBITS - 1); 54 for (b += EXT_FRACHMBITS; !(u.bits.ext_fraclm & m); 55 m >>= 1) 56 b++; 57 #endif /* defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) */ 58 } else { 59 m = 1lu << (EXT_FRACHBITS - 1); 60 for (b = 0; !(u.bits.ext_frach & m); m >>= 1) 61 b++; 62 #ifdef EXT_FRACHMBITS 63 m = 1lu << (EXT_FRACHMBITS - 1); 64 for (; !(u.bits.ext_frachm & m); m >>= 1) 65 b++; 66 #endif /* EXT_FRACHMBITS */ 67 } 68 #ifdef EXT_IMPLICIT_NBIT 69 b++; 70 #endif 71 return ((long double)(LDBL_MIN_EXP - b - 1)); 72 } 73 if (u.bits.ext_exp < (LDBL_MAX_EXP << 1) - 1) /* normal */ 74 return ((long double)(u.bits.ext_exp - LDBL_MAX_EXP + 1)); 75 else /* +/- inf or nan */ 76 return (x * x); 77 } 78