Home | History | Annotate | Download | only in msun
      1 
      2 /* @(#)e_log10.c 1.3 95/01/18 */
      3 /*
      4  * ====================================================
      5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      6  *
      7  * Developed at SunSoft, a Sun Microsystems, Inc. business.
      8  * Permission to use, copy, modify, and distribute this
      9  * software is freely granted, provided that this notice
     10  * is preserved.
     11  * ====================================================
     12  */
     13 
     14 #include <sys/cdefs.h>
     15 __FBSDID("$FreeBSD$");
     16 
     17 /*
     18  * Return the base 2 logarithm of x.  See e_log.c and k_log.h for most
     19  * comments.
     20  *
     21  * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
     22  * then does the combining and scaling steps
     23  *    log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
     24  * in not-quite-routine extra precision.
     25  */
     26 
     27 #include <float.h>
     28 
     29 #include "math.h"
     30 #include "math_private.h"
     31 #include "k_log.h"
     32 
     33 static const double
     34 two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
     35 ivln2hi    =  1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
     36 ivln2lo    =  1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
     37 
     38 static const double zero   =  0.0;
     39 static volatile double vzero = 0.0;
     40 
     41 double
     42 __ieee754_log2(double x)
     43 {
     44 	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
     45 	int32_t i,k,hx;
     46 	u_int32_t lx;
     47 
     48 	EXTRACT_WORDS(hx,lx,x);
     49 
     50 	k=0;
     51 	if (hx < 0x00100000) {			/* x < 2**-1022  */
     52 	    if (((hx&0x7fffffff)|lx)==0)
     53 		return -two54/vzero;		/* log(+-0)=-inf */
     54 	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
     55 	    k -= 54; x *= two54; /* subnormal number, scale up x */
     56 	    GET_HIGH_WORD(hx,x);
     57 	}
     58 	if (hx >= 0x7ff00000) return x+x;
     59 	if (hx == 0x3ff00000 && lx == 0)
     60 	    return zero;			/* log(1) = +0 */
     61 	k += (hx>>20)-1023;
     62 	hx &= 0x000fffff;
     63 	i = (hx+0x95f64)&0x100000;
     64 	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
     65 	k += (i>>20);
     66 	y = (double)k;
     67 	f = x - 1.0;
     68 	hfsq = 0.5*f*f;
     69 	r = k_log1p(f);
     70 
     71 	/*
     72 	 * f-hfsq must (for args near 1) be evaluated in extra precision
     73 	 * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
     74 	 * This is fairly efficient since f-hfsq only depends on f, so can
     75 	 * be evaluated in parallel with R.  Not combining hfsq with R also
     76 	 * keeps R small (though not as small as a true `lo' term would be),
     77 	 * so that extra precision is not needed for terms involving R.
     78 	 *
     79 	 * Compiler bugs involving extra precision used to break Dekker's
     80 	 * theorem for spitting f-hfsq as hi+lo, unless double_t was used
     81 	 * or the multi-precision calculations were avoided when double_t
     82 	 * has extra precision.  These problems are now automatically
     83 	 * avoided as a side effect of the optimization of combining the
     84 	 * Dekker splitting step with the clear-low-bits step.
     85 	 *
     86 	 * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
     87 	 * precision to avoid a very large cancellation when x is very near
     88 	 * these values.  Unlike the above cancellations, this problem is
     89 	 * specific to base 2.  It is strange that adding +-1 is so much
     90 	 * harder than adding +-ln2 or +-log10_2.
     91 	 *
     92 	 * This uses Dekker's theorem to normalize y+val_hi, so the
     93 	 * compiler bugs are back in some configurations, sigh.  And I
     94 	 * don't want to used double_t to avoid them, since that gives a
     95 	 * pessimization and the support for avoiding the pessimization
     96 	 * is not yet available.
     97 	 *
     98 	 * The multi-precision calculations for the multiplications are
     99 	 * routine.
    100 	 */
    101 	hi = f - hfsq;
    102 	SET_LOW_WORD(hi,0);
    103 	lo = (f - hi) - hfsq + r;
    104 	val_hi = hi*ivln2hi;
    105 	val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;
    106 
    107 	/* spadd(val_hi, val_lo, y), except for not using double_t: */
    108 	w = y + val_hi;
    109 	val_lo += (y - w) + val_hi;
    110 	val_hi = w;
    111 
    112 	return val_lo + val_hi;
    113 }
    114 
    115 #if (LDBL_MANT_DIG == 53)
    116 #define __weak_reference(sym,alias) \
    117     __asm__(".weak " #alias); \
    118     __asm__(".equ "  #alias ", " #sym)
    119 __weak_reference(log2, log2l);
    120 #endif
    121