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 <float.h> 13 #include <limits.h> 14 #include <math.h> 15 16 #include "fpmath.h" 17 18 long double 19 logbl(long double x) 20 { 21 union IEEEl2bits u; 22 unsigned long m; 23 int b; 24 25 u.e = x; 26 if (u.bits.exp == 0) { 27 if ((u.bits.manl | u.bits.manh) == 0) { /* x == 0 */ 28 u.bits.sign = 1; 29 return (1.0L / u.e); 30 } 31 /* denormalized */ 32 if (u.bits.manh == 0) { 33 m = 1lu << (LDBL_MANL_SIZE - 1); 34 for (b = LDBL_MANH_SIZE; !(u.bits.manl & m); m >>= 1) 35 b++; 36 } else { 37 m = 1lu << (LDBL_MANH_SIZE - 1); 38 for (b = 0; !(u.bits.manh & m); m >>= 1) 39 b++; 40 } 41 #ifdef LDBL_IMPLICIT_NBIT 42 b++; 43 #endif 44 return ((long double)(LDBL_MIN_EXP - b - 1)); 45 } 46 if (u.bits.exp < (LDBL_MAX_EXP << 1) - 1) /* normal */ 47 return ((long double)(u.bits.exp - LDBL_MAX_EXP + 1)); 48 else /* +/- inf or nan */ 49 return (x * x); 50 } 51