Source file src/pkg/math/log.go

     1	// Copyright 2009 The Go Authors. All rights reserved.
     2	// Use of this source code is governed by a BSD-style
     3	// license that can be found in the LICENSE file.
     4	
     5	package math
     6	
     7	/*
     8		Floating-point logarithm.
     9	*/
    10	
    11	// The original C code, the long comment, and the constants
    12	// below are from FreeBSD's /usr/src/lib/msun/src/e_log.c
    13	// and came with this notice.  The go code is a simpler
    14	// version of the original C.
    15	//
    16	// ====================================================
    17	// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    18	//
    19	// Developed at SunPro, a Sun Microsystems, Inc. business.
    20	// Permission to use, copy, modify, and distribute this
    21	// software is freely granted, provided that this notice
    22	// is preserved.
    23	// ====================================================
    24	//
    25	// __ieee754_log(x)
    26	// Return the logarithm of x
    27	//
    28	// Method :
    29	//   1. Argument Reduction: find k and f such that
    30	//			x = 2**k * (1+f),
    31	//	   where  sqrt(2)/2 < 1+f < sqrt(2) .
    32	//
    33	//   2. Approximation of log(1+f).
    34	//	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
    35	//		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
    36	//	     	 = 2s + s*R
    37	//      We use a special Reme algorithm on [0,0.1716] to generate
    38	//	a polynomial of degree 14 to approximate R.  The maximum error
    39	//	of this polynomial approximation is bounded by 2**-58.45. In
    40	//	other words,
    41	//		        2      4      6      8      10      12      14
    42	//	    R(z) ~ L1*s +L2*s +L3*s +L4*s +L5*s  +L6*s  +L7*s
    43	//	(the values of L1 to L7 are listed in the program) and
    44	//	    |      2          14          |     -58.45
    45	//	    | L1*s +...+L7*s    -  R(z) | <= 2
    46	//	    |                             |
    47	//	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
    48	//	In order to guarantee error in log below 1ulp, we compute log by
    49	//		log(1+f) = f - s*(f - R)		(if f is not too large)
    50	//		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
    51	//
    52	//	3. Finally,  log(x) = k*Ln2 + log(1+f).
    53	//			    = k*Ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*Ln2_lo)))
    54	//	   Here Ln2 is split into two floating point number:
    55	//			Ln2_hi + Ln2_lo,
    56	//	   where n*Ln2_hi is always exact for |n| < 2000.
    57	//
    58	// Special cases:
    59	//	log(x) is NaN with signal if x < 0 (including -INF) ;
    60	//	log(+INF) is +INF; log(0) is -INF with signal;
    61	//	log(NaN) is that NaN with no signal.
    62	//
    63	// Accuracy:
    64	//	according to an error analysis, the error is always less than
    65	//	1 ulp (unit in the last place).
    66	//
    67	// Constants:
    68	// The hexadecimal values are the intended ones for the following
    69	// constants. The decimal values may be used, provided that the
    70	// compiler will convert from decimal to binary accurately enough
    71	// to produce the hexadecimal values shown.
    72	
    73	// Log returns the natural logarithm of x.
    74	//
    75	// Special cases are:
    76	//	Log(+Inf) = +Inf
    77	//	Log(0) = -Inf
    78	//	Log(x < 0) = NaN
    79	//	Log(NaN) = NaN
    80	func Log(x float64) float64
    81	
    82	func log(x float64) float64 {
    83		const (
    84			Ln2Hi = 6.93147180369123816490e-01 /* 3fe62e42 fee00000 */
    85			Ln2Lo = 1.90821492927058770002e-10 /* 3dea39ef 35793c76 */
    86			L1    = 6.666666666666735130e-01   /* 3FE55555 55555593 */
    87			L2    = 3.999999999940941908e-01   /* 3FD99999 9997FA04 */
    88			L3    = 2.857142874366239149e-01   /* 3FD24924 94229359 */
    89			L4    = 2.222219843214978396e-01   /* 3FCC71C5 1D8E78AF */
    90			L5    = 1.818357216161805012e-01   /* 3FC74664 96CB03DE */
    91			L6    = 1.531383769920937332e-01   /* 3FC39A09 D078C69F */
    92			L7    = 1.479819860511658591e-01   /* 3FC2F112 DF3E5244 */
    93		)
    94	
    95		// special cases
    96		switch {
    97		case IsNaN(x) || IsInf(x, 1):
    98			return x
    99		case x < 0:
   100			return NaN()
   101		case x == 0:
   102			return Inf(-1)
   103		}
   104	
   105		// reduce
   106		f1, ki := Frexp(x)
   107		if f1 < Sqrt2/2 {
   108			f1 *= 2
   109			ki--
   110		}
   111		f := f1 - 1
   112		k := float64(ki)
   113	
   114		// compute
   115		s := f / (2 + f)
   116		s2 := s * s
   117		s4 := s2 * s2
   118		t1 := s2 * (L1 + s4*(L3+s4*(L5+s4*L7)))
   119		t2 := s4 * (L2 + s4*(L4+s4*L6))
   120		R := t1 + t2
   121		hfsq := 0.5 * f * f
   122		return k*Ln2Hi - ((hfsq - (s*(hfsq+R) + k*Ln2Lo)) - f)
   123	}