1 /* 2 * From: @(#)s_ilogb.c 5.1 93/09/24 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 13 #include <float.h> 14 #include <limits.h> 15 #include <openlibm.h> 16 #include "math_private.h" 17 #include "fpmath.h" 18 19 DLLEXPORT long double 20 logbl(long double x) 21 { 22 union IEEEl2bits u; 23 unsigned long m; 24 int b; 25 26 u.e = x; 27 if (u.bits.exp == 0) { 28 if ((u.bits.manl | u.bits.manh) == 0) { /* x == 0 */ 29 u.bits.sign = 1; 30 return (1.0L / u.e); 31 } 32 /* denormalized */ 33 if (u.bits.manh == 0) { 34 m = 1lu << (LDBL_MANL_SIZE - 1); 35 for (b = LDBL_MANH_SIZE; !(u.bits.manl & m); m >>= 1) 36 b++; 37 } else { 38 m = 1lu << (LDBL_MANH_SIZE - 1); 39 for (b = 0; !(u.bits.manh & m); m >>= 1) 40 b++; 41 } 42 #ifdef LDBL_IMPLICIT_NBIT 43 b++; 44 #endif 45 return ((long double)(LDBL_MIN_EXP - b - 1)); 46 } 47 if (u.bits.exp < (LDBL_MAX_EXP << 1) - 1) /* normal */ 48 return ((long double)(u.bits.exp - LDBL_MAX_EXP + 1)); 49 else /* +/- inf or nan */ 50 return (x * x); 51 } 52