1 /* 2 * ==================================================== 3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * 5 * Developed at SunPro, a Sun Microsystems, Inc. business. 6 * Permission to use, copy, modify, and distribute this 7 * software is freely granted, provided that this notice 8 * is preserved. 9 * ==================================================== 10 */ 11 12 #include <sys/cdefs.h> 13 #include <float.h> 14 #include <limits.h> 15 #include <math.h> 16 17 #include "fpmath.h" 18 19 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