1141296Sdas
2141296Sdas/* @(#)e_log10.c 1.3 95/01/18 */
32116Sjkh/*
42116Sjkh * ====================================================
52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
62116Sjkh *
7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business.
82116Sjkh * Permission to use, copy, modify, and distribute this
9141296Sdas * software is freely granted, provided that this notice
102116Sjkh * is preserved.
112116Sjkh * ====================================================
122116Sjkh */
132116Sjkh
14176451Sdas#include <sys/cdefs.h>
15176451Sdas__FBSDID("$FreeBSD$");
162116Sjkh
17216247Sdas/*
18216247Sdas * Return the base 2 logarithm of x. See k_log.c for details on the algorithm.
192116Sjkh */
202116Sjkh
212116Sjkh#include "math.h"
222116Sjkh#include "math_private.h"
23216211Sdas#include "k_log.h"
242116Sjkh
252116Sjkhstatic const double
262116Sjkhtwo54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
27216247Sdasivln2hi    =  1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
28216247Sdasivln2lo    =  1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
292116Sjkh
302116Sjkhstatic const double zero   =  0.0;
312116Sjkh
3297413Salfreddouble
33216211Sdas__ieee754_log2(double x)
342116Sjkh{
35216211Sdas	double f,hi,lo;
362116Sjkh	int32_t i,k,hx;
372116Sjkh	u_int32_t lx;
382116Sjkh
392116Sjkh	EXTRACT_WORDS(hx,lx,x);
402116Sjkh
412116Sjkh        k=0;
422116Sjkh        if (hx < 0x00100000) {                  /* x < 2**-1022  */
432116Sjkh            if (((hx&0x7fffffff)|lx)==0)
442116Sjkh                return -two54/zero;             /* log(+-0)=-inf */
452116Sjkh            if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
462116Sjkh            k -= 54; x *= two54; /* subnormal number, scale up x */
472116Sjkh	    GET_HIGH_WORD(hx,x);
482116Sjkh        }
492116Sjkh	if (hx >= 0x7ff00000) return x+x;
502116Sjkh	k += (hx>>20)-1023;
51216211Sdas	hx &= 0x000fffff;
52216211Sdas	i = (hx+0x95f64)&0x100000;
53216211Sdas	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
54216211Sdas	k += (i>>20);
55216211Sdas	f = __kernel_log(x);
56216211Sdas	hi = x = x - 1;
57216211Sdas	SET_LOW_WORD(hi,0);
58216211Sdas	lo = x - hi;
59216211Sdas	return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k;
602116Sjkh}
61