lib/msun/src/e_log2.c

   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 k_log.c for details on the algorithm.
  19  */
  20
  21 #include "math.h"
  22 #include "math_private.h"
  23 #include "k_log.h"
  24
  25 static const double
  26 two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  27 ivln2hi    =  1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
  28 ivln2lo    =  1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
  29
  30 static const double zero   =  0.0;
  31
  32 double
  33 __ieee754_log2(double x)
  34 {
  35         double f,hi,lo;
  36         int32_t i,k,hx;
  37         u_int32_t lx;
  38
  39         EXTRACT_WORDS(hx,lx,x);
  40
  41         k=0;
  42         if (hx < 0x00100000) {                  /* x < 2**-1022  */
  43             if (((hx&0x7fffffff)|lx)==0)
  44                 return -two54/zero;             /* log(+-0)=-inf */
  45             if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
  46             k -= 54; x *= two54; /* subnormal number, scale up x */
  47             GET_HIGH_WORD(hx,x);
  48         }
  49         if (hx >= 0x7ff00000) return x+x;
  50         k += (hx>>20)-1023;
  51         hx &= 0x000fffff;
  52         i = (hx+0x95f64)&0x100000;
  53         SET_HIGH_WORD(x,hx|(i^0x3ff00000));     /* normalize x or x/2 */
  54         k += (i>>20);
  55         f = __kernel_log(x);
  56         hi = x = x - 1;
  57         SET_LOW_WORD(hi,0);
  58         lo = x - hi;
  59         return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k;
  60 }