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
17219360Sdas/*
18219360Sdas * Return the base 10 logarithm of x. See k_log.c for details on the algorithm.
192116Sjkh */
202116Sjkh
212116Sjkh#include "math.h"
222116Sjkh#include "math_private.h"
23219360Sdas#include "k_log.h"
242116Sjkh
252116Sjkhstatic const double
262116Sjkhtwo54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
27219360Sdasivln10hi   =  4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
28219360Sdasivln10lo   =  2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
292116Sjkhlog10_2hi  =  3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
302116Sjkhlog10_2lo  =  3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
312116Sjkh
322116Sjkhstatic const double zero   =  0.0;
332116Sjkh
3497413Salfreddouble
35117912Speter__ieee754_log10(double x)
362116Sjkh{
37219360Sdas	double f,hi,lo,y,z;
382116Sjkh	int32_t i,k,hx;
392116Sjkh	u_int32_t lx;
402116Sjkh
412116Sjkh	EXTRACT_WORDS(hx,lx,x);
422116Sjkh
432116Sjkh        k=0;
442116Sjkh        if (hx < 0x00100000) {                  /* x < 2**-1022  */
452116Sjkh            if (((hx&0x7fffffff)|lx)==0)
462116Sjkh                return -two54/zero;             /* log(+-0)=-inf */
472116Sjkh            if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
482116Sjkh            k -= 54; x *= two54; /* subnormal number, scale up x */
492116Sjkh	    GET_HIGH_WORD(hx,x);
502116Sjkh        }
512116Sjkh	if (hx >= 0x7ff00000) return x+x;
522116Sjkh	k += (hx>>20)-1023;
53219360Sdas	hx &= 0x000fffff;
54219360Sdas	i = (hx+0x95f64)&0x100000;
55219360Sdas	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
56219360Sdas	k += (i>>20);
57219360Sdas	y = (double)k;
58219360Sdas	f = __kernel_log(x);
59219360Sdas	hi = x = x - 1;
60219360Sdas	SET_LOW_WORD(hi,0);
61219360Sdas	lo = x - hi;
62219360Sdas	z = y*log10_2lo + (x+f)*ivln10lo + (lo+f)*ivln10hi + hi*ivln10hi;
632116Sjkh	return  z+y*log10_2hi;
642116Sjkh}
65