1177668d1SDavid Schultz 2177668d1SDavid Schultz /* @(#)e_log10.c 1.3 95/01/18 */ 3177668d1SDavid Schultz /* 4177668d1SDavid Schultz * ==================================================== 5177668d1SDavid Schultz * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6177668d1SDavid Schultz * 7177668d1SDavid Schultz * Developed at SunSoft, a Sun Microsystems, Inc. business. 8177668d1SDavid Schultz * Permission to use, copy, modify, and distribute this 9177668d1SDavid Schultz * software is freely granted, provided that this notice 10177668d1SDavid Schultz * is preserved. 11177668d1SDavid Schultz * ==================================================== 12177668d1SDavid Schultz */ 13177668d1SDavid Schultz 14177668d1SDavid Schultz #include <sys/cdefs.h> 15177668d1SDavid Schultz __FBSDID("$FreeBSD$"); 16177668d1SDavid Schultz 17*63687c8bSDavid Schultz /* 18*63687c8bSDavid Schultz * Return the base 2 logarithm of x. See k_log.c for details on the algorithm. 19177668d1SDavid Schultz */ 20177668d1SDavid Schultz 21177668d1SDavid Schultz #include "math.h" 22177668d1SDavid Schultz #include "math_private.h" 23177668d1SDavid Schultz #include "k_log.h" 24177668d1SDavid Schultz 25177668d1SDavid Schultz static const double 26177668d1SDavid Schultz two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ 27*63687c8bSDavid Schultz ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ 28*63687c8bSDavid Schultz ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ 29177668d1SDavid Schultz 30177668d1SDavid Schultz static const double zero = 0.0; 31177668d1SDavid Schultz 32177668d1SDavid Schultz double 33177668d1SDavid Schultz __ieee754_log2(double x) 34177668d1SDavid Schultz { 35177668d1SDavid Schultz double f,hi,lo; 36177668d1SDavid Schultz int32_t i,k,hx; 37177668d1SDavid Schultz u_int32_t lx; 38177668d1SDavid Schultz 39177668d1SDavid Schultz EXTRACT_WORDS(hx,lx,x); 40177668d1SDavid Schultz 41177668d1SDavid Schultz k=0; 42177668d1SDavid Schultz if (hx < 0x00100000) { /* x < 2**-1022 */ 43177668d1SDavid Schultz if (((hx&0x7fffffff)|lx)==0) 44177668d1SDavid Schultz return -two54/zero; /* log(+-0)=-inf */ 45177668d1SDavid Schultz if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ 46177668d1SDavid Schultz k -= 54; x *= two54; /* subnormal number, scale up x */ 47177668d1SDavid Schultz GET_HIGH_WORD(hx,x); 48177668d1SDavid Schultz } 49177668d1SDavid Schultz if (hx >= 0x7ff00000) return x+x; 50177668d1SDavid Schultz k += (hx>>20)-1023; 51177668d1SDavid Schultz hx &= 0x000fffff; 52177668d1SDavid Schultz i = (hx+0x95f64)&0x100000; 53177668d1SDavid Schultz SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ 54177668d1SDavid Schultz k += (i>>20); 55177668d1SDavid Schultz f = __kernel_log(x); 56177668d1SDavid Schultz hi = x = x - 1; 57177668d1SDavid Schultz SET_LOW_WORD(hi,0); 58177668d1SDavid Schultz lo = x - hi; 59177668d1SDavid Schultz return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k; 60177668d1SDavid Schultz } 61