1*177668d1SDavid Schultz 2*177668d1SDavid Schultz /* @(#)e_log10.c 1.3 95/01/18 */ 3*177668d1SDavid Schultz /* 4*177668d1SDavid Schultz * ==================================================== 5*177668d1SDavid Schultz * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6*177668d1SDavid Schultz * 7*177668d1SDavid Schultz * Developed at SunSoft, a Sun Microsystems, Inc. business. 8*177668d1SDavid Schultz * Permission to use, copy, modify, and distribute this 9*177668d1SDavid Schultz * software is freely granted, provided that this notice 10*177668d1SDavid Schultz * is preserved. 11*177668d1SDavid Schultz * ==================================================== 12*177668d1SDavid Schultz */ 13*177668d1SDavid Schultz 14*177668d1SDavid Schultz #include <sys/cdefs.h> 15*177668d1SDavid Schultz __FBSDID("$FreeBSD$"); 16*177668d1SDavid Schultz 17*177668d1SDavid Schultz /* log2(x) 18*177668d1SDavid Schultz * Return the base 2 logarithm of x. 19*177668d1SDavid Schultz */ 20*177668d1SDavid Schultz 21*177668d1SDavid Schultz #include "math.h" 22*177668d1SDavid Schultz #include "math_private.h" 23*177668d1SDavid Schultz #include "k_log.h" 24*177668d1SDavid Schultz 25*177668d1SDavid Schultz static const double 26*177668d1SDavid Schultz two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ 27*177668d1SDavid Schultz ivln2hi = 0x1.71547652000p+0, 28*177668d1SDavid Schultz ivln2lo = 0x1.705fc2eefa2p-33; 29*177668d1SDavid Schultz 30*177668d1SDavid Schultz static const double zero = 0.0; 31*177668d1SDavid Schultz 32*177668d1SDavid Schultz double 33*177668d1SDavid Schultz __ieee754_log2(double x) 34*177668d1SDavid Schultz { 35*177668d1SDavid Schultz double f,hi,lo; 36*177668d1SDavid Schultz int32_t i,k,hx; 37*177668d1SDavid Schultz u_int32_t lx; 38*177668d1SDavid Schultz 39*177668d1SDavid Schultz EXTRACT_WORDS(hx,lx,x); 40*177668d1SDavid Schultz 41*177668d1SDavid Schultz k=0; 42*177668d1SDavid Schultz if (hx < 0x00100000) { /* x < 2**-1022 */ 43*177668d1SDavid Schultz if (((hx&0x7fffffff)|lx)==0) 44*177668d1SDavid Schultz return -two54/zero; /* log(+-0)=-inf */ 45*177668d1SDavid Schultz if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ 46*177668d1SDavid Schultz k -= 54; x *= two54; /* subnormal number, scale up x */ 47*177668d1SDavid Schultz GET_HIGH_WORD(hx,x); 48*177668d1SDavid Schultz } 49*177668d1SDavid Schultz if (hx >= 0x7ff00000) return x+x; 50*177668d1SDavid Schultz k += (hx>>20)-1023; 51*177668d1SDavid Schultz hx &= 0x000fffff; 52*177668d1SDavid Schultz i = (hx+0x95f64)&0x100000; 53*177668d1SDavid Schultz SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ 54*177668d1SDavid Schultz k += (i>>20); 55*177668d1SDavid Schultz f = __kernel_log(x); 56*177668d1SDavid Schultz hi = x = x - 1; 57*177668d1SDavid Schultz SET_LOW_WORD(hi,0); 58*177668d1SDavid Schultz lo = x - hi; 59*177668d1SDavid Schultz return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k; 60*177668d1SDavid Schultz } 61