1*390c8400Smartynas /* $OpenBSD: s_logbl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $ */ 2*390c8400Smartynas /* 3*390c8400Smartynas * From: @(#)s_ilogb.c 5.1 93/09/24 4*390c8400Smartynas * ==================================================== 5*390c8400Smartynas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6*390c8400Smartynas * 7*390c8400Smartynas * Developed at SunPro, a Sun Microsystems, Inc. business. 8*390c8400Smartynas * Permission to use, copy, modify, and distribute this 9*390c8400Smartynas * software is freely granted, provided that this notice 10*390c8400Smartynas * is preserved. 11*390c8400Smartynas * ==================================================== 12*390c8400Smartynas */ 13*390c8400Smartynas 14*390c8400Smartynas #include <sys/types.h> 15*390c8400Smartynas #include <machine/ieee.h> 16*390c8400Smartynas #include <float.h> 17*390c8400Smartynas #include <limits.h> 18*390c8400Smartynas #include <math.h> 19*390c8400Smartynas 20*390c8400Smartynas long double logbl(long double x)21*390c8400Smartynaslogbl(long double x) 22*390c8400Smartynas { 23*390c8400Smartynas union { 24*390c8400Smartynas long double e; 25*390c8400Smartynas struct ieee_ext bits; 26*390c8400Smartynas } u; 27*390c8400Smartynas unsigned long m; 28*390c8400Smartynas int b; 29*390c8400Smartynas 30*390c8400Smartynas u.e = x; 31*390c8400Smartynas if (u.bits.ext_exp == 0) { 32*390c8400Smartynas if ((u.bits.ext_fracl 33*390c8400Smartynas #ifdef EXT_FRACLMBITS 34*390c8400Smartynas | u.bits.ext_fraclm 35*390c8400Smartynas #endif /* EXT_FRACLMBITS */ 36*390c8400Smartynas #ifdef EXT_FRACHMBITS 37*390c8400Smartynas | u.bits.ext_frachm 38*390c8400Smartynas #endif /* EXT_FRACHMBITS */ 39*390c8400Smartynas | u.bits.ext_frach) == 0) { /* x == 0 */ 40*390c8400Smartynas u.bits.ext_sign = 1; 41*390c8400Smartynas return (1.0L / u.e); 42*390c8400Smartynas } 43*390c8400Smartynas /* denormalized */ 44*390c8400Smartynas if (u.bits.ext_frach == 0 45*390c8400Smartynas #ifdef EXT_FRACHMBITS 46*390c8400Smartynas && u.bits.ext_frachm == 0 47*390c8400Smartynas #endif 48*390c8400Smartynas ) { 49*390c8400Smartynas m = 1lu << (EXT_FRACLBITS - 1); 50*390c8400Smartynas for (b = EXT_FRACHBITS; !(u.bits.ext_fracl & m); m >>= 1) 51*390c8400Smartynas b++; 52*390c8400Smartynas #if defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) 53*390c8400Smartynas m = 1lu << (EXT_FRACLMBITS - 1); 54*390c8400Smartynas for (b += EXT_FRACHMBITS; !(u.bits.ext_fraclm & m); 55*390c8400Smartynas m >>= 1) 56*390c8400Smartynas b++; 57*390c8400Smartynas #endif /* defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) */ 58*390c8400Smartynas } else { 59*390c8400Smartynas m = 1lu << (EXT_FRACHBITS - 1); 60*390c8400Smartynas for (b = 0; !(u.bits.ext_frach & m); m >>= 1) 61*390c8400Smartynas b++; 62*390c8400Smartynas #ifdef EXT_FRACHMBITS 63*390c8400Smartynas m = 1lu << (EXT_FRACHMBITS - 1); 64*390c8400Smartynas for (; !(u.bits.ext_frachm & m); m >>= 1) 65*390c8400Smartynas b++; 66*390c8400Smartynas #endif /* EXT_FRACHMBITS */ 67*390c8400Smartynas } 68*390c8400Smartynas #ifdef EXT_IMPLICIT_NBIT 69*390c8400Smartynas b++; 70*390c8400Smartynas #endif 71*390c8400Smartynas return ((long double)(LDBL_MIN_EXP - b - 1)); 72*390c8400Smartynas } 73*390c8400Smartynas if (u.bits.ext_exp < (LDBL_MAX_EXP << 1) - 1) /* normal */ 74*390c8400Smartynas return ((long double)(u.bits.ext_exp - LDBL_MAX_EXP + 1)); 75*390c8400Smartynas else /* +/- inf or nan */ 76*390c8400Smartynas return (x * x); 77*390c8400Smartynas } 78