math: add functions; update tests and special cases

Added special cases to comments for asin.go and fabs.go.
Added Trunc() to floor.go and floor_386.s.  Fixed formatting
error in hypot_386.s  Added new functions Acosh, Asinh,
Atanh, Copysign, Erf, Erfc, Expm1, and Log1p.  Added
386 FPU version of Fmod.  Added tests, benchmarks, and
precision to expected results in all_test.go.  Edited
makefile so it all compiles.

R=rsc
CC=golang-dev
https://golang.org/cl/195052

src/pkg/math/Makefile

@@ -15,6 +15,7 @@ OFILES_386=\
 	exp_386.$O\
 	fabs_386.$O\
 	floor_386.$O\
+	fmod_386.$O\
 	hypot_386.$O\
 	log_386.$O\
 	sin_386.$O\
@@ -25,17 +26,24 @@ OFILES=\
 	$(OFILES_$(GOARCH))
 
 ALLGOFILES=\
+	acosh.go\
 	asin.go\
+	asinh.go\
 	atan.go\
+	atanh.go\
 	atan2.go\
 	bits.go\
 	const.go\
+	copysign.go\
+	erf.go\
 	exp.go\
+	expm1.go\
 	fabs.go\
 	floor.go\
 	fmod.go\
 	hypot.go\
 	log.go\
+	log1p.go\
 	pow.go\
 	pow10.go\
 	sin.go\

src/pkg/math/acosh.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/e_acosh.c
+// and came with this notice.  The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// __ieee754_acosh(x)
+// Method :
+//	Based on
+//	        acosh(x) = log [ x + sqrt(x*x-1) ]
+//	we have
+//	        acosh(x) := log(x)+ln2,	if x is large; else
+//	        acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
+//	        acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
+//
+// Special cases:
+//	acosh(x) is NaN with signal if x<1.
+//	acosh(NaN) is NaN without signal.
+//
+
+// Acosh(x) calculates the inverse hyperbolic cosine of x.
+//
+// Special cases are:
+//	Acosh(x) = NaN if x < 1
+//	Acosh(NaN) = NaN
+func Acosh(x float64) float64 {
+	const (
+		Ln2   = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
+		Large = 1 << 28                    // 2^28
+	)
+	switch {
+	case x < 1 || IsNaN(x):
+		return NaN()
+	case x == 1:
+		return 0
+	case x >= Large:
+		return Log(x) + Ln2 // x > 2^28
+	case x > 2:
+		return Log(2*x - 1/(x+Sqrt(x*x-1))) // 2^28 > x > 2
+	}
+	t := x - 1
+	return Log1p(t + Sqrt(2*t+t*t)) // 2 >= x > 1
+}

src/pkg/math/asin.go

@@ -6,13 +6,16 @@ package math
 
 
 /*
-	Floating-point sine and cosine.
+	Floating-point arcsine and arccosine.
 
 	They are implemented by computing the arctangent
 	after appropriate range reduction.
 */
 
 // Asin returns the arcsine of x.
+//
+// Special case is:
+//	Asin(x) = NaN if x < -1 or x > 1
 func Asin(x float64) float64 {
 	sign := false
 	if x < 0 {
@@ -37,4 +40,7 @@ func Asin(x float64) float64 {
 }
 
 // Acos returns the arccosine of x.
+//
+// Special case is:
+//	Acos(x) = NaN if x < -1 or x > 1
 func Acos(x float64) float64 { return Pi/2 - Asin(x) }

src/pkg/math/asinh.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/s_asinh.c
+// and came with this notice.  The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// asinh(x)
+// Method :
+//	Based on
+//	        asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
+//	we have
+//	asinh(x) := x  if  1+x*x=1,
+//	         := sign(x)*(log(x)+ln2)) for large |x|, else
+//	         := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
+//	         := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
+//
+
+// Asinh(x) calculates the inverse hyperbolic sine of x.
+//
+// Special cases are:
+//	Asinh(+Inf) = +Inf
+//	Asinh(-Inf) = -Inf
+//	Asinh(NaN) = NaN
+func Asinh(x float64) float64 {
+	const (
+		Ln2      = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
+		NearZero = 1.0 / (1 << 28)            // 2^-28
+		Large    = 1 << 28                    // 2^28
+	)
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	// special cases
+	if x != x || x > MaxFloat64 || x < -MaxFloat64 { // IsNaN(x) || IsInf(x, 0)
+		return x
+	}
+	sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+	var temp float64
+	switch {
+	case x > Large:
+		temp = Log(x) + Ln2 // |x| > 2^28
+	case x > 2:
+		temp = Log(2*x + 1/(Sqrt(x*x+1)+x)) // 2^28 > |x| > 2.0
+	case x < NearZero:
+		temp = x // |x| < 2^-28
+	default:
+		temp = Log1p(x + x*x/(1+Sqrt(1+x*x))) // 2.0 > |x| > 2^-28
+	}
+	if sign {
+		temp = -temp
+	}
+	return temp
+}

src/pkg/math/atanh.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/e_atanh.c
+// and came with this notice.  The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// __ieee754_atanh(x)
+// Method :
+//	1. Reduce x to positive by atanh(-x) = -atanh(x)
+//	2. For x>=0.5
+//	            1              2x                          x
+//	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+//	            2             1 - x                      1 - x
+//
+//	For x<0.5
+//	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+//
+// Special cases:
+//	atanh(x) is NaN if |x| > 1 with signal;
+//	atanh(NaN) is that NaN with no signal;
+//	atanh(+-1) is +-INF with signal.
+//
+
+// Atanh(x) calculates the inverse hyperbolic tangent of x.
+//
+// Special cases are:
+//	Atanh(x) = NaN if x < -1 or x > 1
+//	Atanh(1) = +Inf
+//	Atanh(-1) = -Inf
+//	Atanh(NaN) = NaN
+func Atanh(x float64) float64 {
+	const NearZero = 1.0 / (1 << 28) // 2^-28
+	// TODO(rsc): Remove manual inlining of IsNaN
+	// when compiler does it for us
+	// special cases
+	switch {
+	case x < -1 || x > 1 || x != x: // x < -1 || x > 1 || IsNaN(x):
+		return NaN()
+	case x == 1:
+		return Inf(1)
+	case x == -1:
+		return Inf(-1)
+	}
+	sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+	var temp float64
+	switch {
+	case x < NearZero:
+		temp = x
+	case x < 0.5:
+		temp = x + x
+		temp = 0.5 * Log1p(temp+temp*x/(1-x))
+	default:
+		temp = 0.5 * Log1p((x+x)/(1-x))
+	}
+	if sign {
+		temp = -temp
+	}
+	return temp
+}

src/pkg/math/copysign.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// Copysign(x, y) returns a value with the magnitude
+// of x and the sign of y.
+func Copysign(x, y float64) float64 {
+	if x < 0 && y > 0 || x > 0 && y < 0 {
+		return -x
+	}
+	return x
+}

src/pkg/math/exp.go

@@ -86,7 +86,7 @@ package math
 // Special cases are:
 //	Exp(+Inf) = +Inf
 //	Exp(NaN) = NaN
-// Very large values overflow to -Inf or +Inf.
+// Very large values overflow to 0 or +Inf.
 // Very small values underflow to 1.
 func Exp(x float64) float64 {
 	const (

src/pkg/math/expm1.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/s_expm1.c
+// and came with this notice.  The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// expm1(x)
+// Returns exp(x)-1, the exponential of x minus 1.
+//
+// Method
+//   1. Argument reduction:
+//      Given x, find r and integer k such that
+//
+//               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658
+//
+//      Here a correction term c will be computed to compensate
+//      the error in r when rounded to a floating-point number.
+//
+//   2. Approximating expm1(r) by a special rational function on
+//      the interval [0,0.34658]:
+//      Since
+//          r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
+//      we define R1(r*r) by
+//          r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
+//      That is,
+//          R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
+//                   = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
+//                   = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
+//      We use a special Reme algorithm on [0,0.347] to generate
+//      a polynomial of degree 5 in r*r to approximate R1. The
+//      maximum error of this polynomial approximation is bounded
+//      by 2**-61. In other words,
+//          R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
+//      where   Q1  =  -1.6666666666666567384E-2,
+//              Q2  =   3.9682539681370365873E-4,
+//              Q3  =  -9.9206344733435987357E-6,
+//              Q4  =   2.5051361420808517002E-7,
+//              Q5  =  -6.2843505682382617102E-9;
+//      (where z=r*r, and the values of Q1 to Q5 are listed below)
+//      with error bounded by
+//          |                  5           |     -61
+//          | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
+//          |                              |
+//
+//      expm1(r) = exp(r)-1 is then computed by the following
+//      specific way which minimize the accumulation rounding error:
+//                             2     3
+//                            r     r    [ 3 - (R1 + R1*r/2)  ]
+//            expm1(r) = r + --- + --- * [--------------------]
+//                            2     2    [ 6 - r*(3 - R1*r/2) ]
+//
+//      To compensate the error in the argument reduction, we use
+//              expm1(r+c) = expm1(r) + c + expm1(r)*c
+//                         ~ expm1(r) + c + r*c
+//      Thus c+r*c will be added in as the correction terms for
+//      expm1(r+c). Now rearrange the term to avoid optimization
+//      screw up:
+//                      (      2                                    2 )
+//                      ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
+//       expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
+//                      ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
+//                      (                                             )
+//
+//                 = r - E
+//   3. Scale back to obtain expm1(x):
+//      From step 1, we have
+//         expm1(x) = either 2^k*[expm1(r)+1] - 1
+//                  = or     2^k*[expm1(r) + (1-2^-k)]
+//   4. Implementation notes:
+//      (A). To save one multiplication, we scale the coefficient Qi
+//           to Qi*2^i, and replace z by (x^2)/2.
+//      (B). To achieve maximum accuracy, we compute expm1(x) by
+//        (i)   if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
+//        (ii)  if k=0, return r-E
+//        (iii) if k=-1, return 0.5*(r-E)-0.5
+//        (iv)  if k=1 if r < -0.25, return 2*((r+0.5)- E)
+//                     else          return  1.0+2.0*(r-E);
+//        (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
+//        (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
+//        (vii) return 2^k(1-((E+2^-k)-r))
+//
+// Special cases:
+//      expm1(INF) is INF, expm1(NaN) is NaN;
+//      expm1(-INF) is -1, and
+//      for finite argument, only expm1(0)=0 is exact.
+//
+// Accuracy:
+//      according to an error analysis, the error is always less than
+//      1 ulp (unit in the last place).
+//
+// Misc. info.
+//      For IEEE double
+//          if x >  7.09782712893383973096e+02 then expm1(x) overflow
+//
+// Constants:
+// The hexadecimal values are the intended ones for the following
+// constants. The decimal values may be used, provided that the
+// compiler will convert from decimal to binary accurately enough
+// to produce the hexadecimal values shown.
+//
+
+// Expm1 returns e^x - 1, the base-e exponential of x minus 1.
+// It is more accurate than Exp(x) - 1 when x is near zero.
+//
+// Special cases are:
+//      Expm1(+Inf) = +Inf
+//      Expm1(-Inf) = -1
+//	Expm1(NaN) = NaN
+// Very large values overflow to -1 or +Inf.
+func Expm1(x float64) float64 {
+	const (
+		Othreshold = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF
+		Ln2X56     = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1
+		Ln2HalfX3  = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73
+		Ln2Half    = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef
+		Ln2Hi      = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
+		Ln2Lo      = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
+		InvLn2     = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
+		Tiny       = 1.0 / (1 << 54)            // 2^-54 = 0x3c90000000000000
+		// scaled coefficients related to expm1
+		Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4
+		Q2 = 1.58730158725481460165e-03  // 0x3F5A01A019FE5585
+		Q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7
+		Q4 = 4.00821782732936239552e-06  // 0x3ED0CFCA86E65239
+		Q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D
+	)
+
+	// special cases
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	switch {
+	case x > MaxFloat64 || x != x: // IsInf(x, 1) || IsNaN(x):
+		return x
+	case x < -MaxFloat64: // IsInf(x, -1):
+		return -1
+	}
+
+	absx := x
+	sign := false
+	if x < 0 {
+		absx = -absx
+		sign = true
+	}
+
+	// filter out huge argument
+	if absx >= Ln2X56 { // if |x| >= 56 * ln2
+		if absx >= Othreshold { // if |x| >= 709.78...
+			return Inf(1) // overflow
+		}
+		if sign {
+			return -1 // x < -56*ln2, return -1.0
+		}
+	}
+
+	// argument reduction
+	var c float64
+	var k int
+	if absx > Ln2Half { // if  |x| > 0.5 * ln2
+		var hi, lo float64
+		if absx < Ln2HalfX3 { // and |x| < 1.5 * ln2
+			if !sign {
+				hi = x - Ln2Hi
+				lo = Ln2Lo
+				k = 1
+			} else {
+				hi = x + Ln2Hi
+				lo = -Ln2Lo
+				k = -1
+			}
+		} else {
+			if !sign {
+				k = int(InvLn2*x + 0.5)
+			} else {
+				k = int(InvLn2*x - 0.5)
+			}
+			t := float64(k)
+			hi = x - t*Ln2Hi // t * Ln2Hi is exact here
+			lo = t * Ln2Lo
+		}
+		x = hi - lo
+		c = (hi - x) - lo
+	} else if absx < Tiny { // when |x| < 2^-54, return x
+		return x
+	} else {
+		k = 0
+	}
+
+	// x is now in primary range
+	hfx := 0.5 * x
+	hxs := x * hfx
+	r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))))
+	t := 3 - r1*hfx
+	e := hxs * ((r1 - t) / (6.0 - x*t))
+	if k != 0 {
+		e = (x*(e-c) - c)
+		e -= hxs
+		switch {
+		case k == -1:
+			return 0.5*(x-e) - 0.5
+		case k == 1:
+			if x < -0.25 {
+				return -2 * (e - (x + 0.5))
+			}
+			return 1 + 2*(x-e)
+		case k <= -2 || k > 56: // suffice to return exp(x)-1
+			y := 1 - (e - x)
+			y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
+			return y - 1
+		}
+		if k < 20 {
+			t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2^-k
+			y := t - (e - x)
+			y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
+			return y
+		}
+		t := Float64frombits(uint64((0x3ff - k) << 52)) // 2^-k
+		y := x - (e + t)
+		y += 1
+		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
+		return y
+	}
+	return x - (x*e - hxs) // c is 0
+}

src/pkg/math/fabs.go

@@ -5,6 +5,11 @@
 package math
 
 // Fabs returns the absolute value of x.
+//
+// Special cases are:
+//	Fabs(+Inf) = +Inf
+//	Fabs(-Inf) = +Inf
+//	Fabs(NaN) = NaN
 func Fabs(x float64) float64 {
 	if x < 0 {
 		return -x

src/pkg/math/floor.go

-// Copyright 2009 The Go Authors. All rights reserved.
+// Copyright 2009-2010 The Go Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
@@ -12,6 +12,8 @@ package math
 //	Floor(-Inf) = -Inf
 //	Floor(NaN) = NaN
 func Floor(x float64) float64 {
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
 	if x != x || x > MaxFloat64 || x < -MaxFloat64 { // IsNaN(x) || IsInf(x, 0)
 		return x
 	}
@@ -33,3 +35,19 @@ func Floor(x float64) float64 {
 //	Ceil(-Inf) = -Inf
 //	Ceil(NaN) = NaN
 func Ceil(x float64) float64 { return -Floor(-x) }
+
+// Trunc returns the integer value of x.
+//
+// Special cases are:
+//	Trunc(+Inf) = +Inf
+//	Trunc(-Inf) = -Inf
+//	Trunc(NaN) = NaN
+func Trunc(x float64) float64 {
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	if x != x || x > MaxFloat64 || x < -MaxFloat64 { // IsNaN(x) || IsInf(x, 0)
+		return x
+	}
+	d, _ := Modf(x)
+	return d
+}

src/pkg/math/floor_386.s

@@ -25,7 +25,20 @@ TEXT ·Floor(SB),7,$0
 	ORW     $0x0400, AX  // Rounding Control set to -Inf
 	MOVW    AX, -4(SP)   // store new Control Word
 	FLDCW   -4(SP)       // load new Control Word
-	FRNDINT              // F0=floor(x)
+	FRNDINT              // F0=Floor(x)
+	FLDCW   -2(SP)       // load old Control Word
+	FMOVDP  F0, r+8(FP)
+	RET
+
+// func Trunc(x float64) float64
+TEXT ·Trunc(SB),7,$0
+	FMOVD   x+0(FP), F0  // F0=x
+	FSTCW   -2(SP)       // save old Control Word
+	MOVW    -2(SP), AX
+	ORW     $0x0c00, AX  // Rounding Control set to truncate
+	MOVW    AX, -4(SP)   // store new Control Word
+	FLDCW   -4(SP)       // load new Control Word
+	FRNDINT              // F0=Trunc(x)
 	FLDCW   -2(SP)       // load old Control Word
 	FMOVDP  F0, r+8(FP)
 	RET

src/pkg/math/floor_decl.go

@@ -6,3 +6,4 @@ package math
 
 func Ceil(x float64) float64
 func Floor(x float64) float64
+func Trunc(x float64) float64

src/pkg/math/fmod_386.s

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// func Fmod(x, y float64) float64
+TEXT ·Fmod(SB),7,$0
+	FMOVD   y+8(FP), F0  // F0=y
+	FMOVD   x+0(FP), F0  // F0=x, F1=y
+	FPREM                // F0=reduced_x, F1=y
+	FSTSW   AX           // AX=status word
+	ANDW    $0x0400, AX
+	JNE     -3(PC)       // jump if reduction incomplete
+	FMOVDP  F0, F1       // F0=x-q*y
+	FMOVDP  F0, r+16(FP)
+	RET

src/pkg/math/fmod_decl.go

+// Copyright 2010 The Go Authors.  All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+func Fmod(x, y float64) float64

src/pkg/math/log1p.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+
+// The original C code, the long comment, and the constants
+// below are from FreeBSD's /usr/src/lib/msun/src/s_log1p.c
+// and came with this notice.  The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+//
+// double log1p(double x)
+//
+// Method :
+//   1. Argument Reduction: find k and f such that
+//                      1+x = 2^k * (1+f),
+//         where  sqrt(2)/2 < 1+f < sqrt(2) .
+//
+//      Note. If k=0, then f=x is exact. However, if k!=0, then f
+//      may not be representable exactly. In that case, a correction
+//      term is need. Let u=1+x rounded. Let c = (1+x)-u, then
+//      log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
+//      and add back the correction term c/u.
+//      (Note: when x > 2**53, one can simply return log(x))
+//
+//   2. Approximation of log1p(f).
+//      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+//               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+//               = 2s + s*R
+//      We use a special Reme algorithm on [0,0.1716] to generate
+//      a polynomial of degree 14 to approximate R The maximum error
+//      of this polynomial approximation is bounded by 2**-58.45. In
+//      other words,
+//                      2      4      6      8      10      12      14
+//          R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s  +Lp6*s  +Lp7*s
+//      (the values of Lp1 to Lp7 are listed in the program)
+//      a-0.2929nd
+//          |      2          14          |     -58.45
+//          | Lp1*s +...+Lp7*s    -  R(z) | <= 2
+//          |                             |
+//      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+//      In order to guarantee error in log below 1ulp, we compute log
+//      by
+//              log1p(f) = f - (hfsq - s*(hfsq+R)).
+//
+//   3. Finally, log1p(x) = k*ln2 + log1p(f).
+//                        = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+//      Here ln2 is split into two floating point number:
+//                   ln2_hi + ln2_lo,
+//      where n*ln2_hi is always exact for |n| < 2000.
+//
+// Special cases:
+//      log1p(x) is NaN with signal if x < -1 (including -INF) ;
+//      log1p(+INF) is +INF; log1p(-1) is -INF with signal;
+//      log1p(NaN) is that NaN with no signal.
+//
+// Accuracy:
+//      according to an error analysis, the error is always less than
+//      1 ulp (unit in the last place).
+//
+// Constants:
+// The hexadecimal values are the intended ones for the following
+// constants. The decimal values may be used, provided that the
+// compiler will convert from decimal to binary accurately enough
+// to produce the hexadecimal values shown.
+//
+// Note: Assuming log() return accurate answer, the following
+//       algorithm can be used to compute log1p(x) to within a few ULP:
+//
+//              u = 1+x;
+//              if(u==1.0) return x ; else
+//                         return log(u)*(x/(u-1.0));
+//
+//       See HP-15C Advanced Functions Handbook, p.193.
+
+// Log1p returns the natural logarithm of 1 plus its argument x.
+// It is more accurate than Log(1 + x) when x is near zero.
+//
+// Special cases are:
+//	Log1p(+Inf) = +Inf
+//	Log1p(-1) = -Inf
+//	Log1p(x < -1) = NaN
+//	Log1p(NaN) = NaN
+func Log1p(x float64) float64 {
+	const (
+		Sqrt2M1     = 4.142135623730950488017e-01  // Sqrt(2)-1 = 0x3fda827999fcef34
+		Sqrt2HalfM1 = -2.928932188134524755992e-01 // Sqrt(2)/2-1 = 0xbfd2bec333018866
+		Small       = 1.0 / (1 << 29)              // 2^-29 = 0x3e20000000000000
+		Tiny        = 1.0 / (1 << 54)              // 2^-54
+		Two53       = 1 << 53                      // 2^53
+		Ln2Hi       = 6.93147180369123816490e-01   // 3fe62e42fee00000
+		Ln2Lo       = 1.90821492927058770002e-10   // 3dea39ef35793c76
+		Lp1         = 6.666666666666735130e-01     // 3FE5555555555593
+		Lp2         = 3.999999999940941908e-01     // 3FD999999997FA04
+		Lp3         = 2.857142874366239149e-01     // 3FD2492494229359
+		Lp4         = 2.222219843214978396e-01     // 3FCC71C51D8E78AF
+		Lp5         = 1.818357216161805012e-01     // 3FC7466496CB03DE
+		Lp6         = 1.531383769920937332e-01     // 3FC39A09D078C69F
+		Lp7         = 1.479819860511658591e-01     // 3FC2F112DF3E5244
+	)
+
+	// special cases
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	switch {
+	case x < -1 || x != x: // x < -1 || IsNaN(x): // includes -Inf
+		return NaN()
+	case x == -1:
+		return Inf(-1)
+	case x > MaxFloat64: // IsInf(x, 1):
+		return Inf(1)
+	}
+
+	absx := x
+	if absx < 0 {
+		absx = -absx
+	}
+
+	var f float64
+	var iu uint64
+	k := 1
+	if absx < Sqrt2M1 { //  |x| < Sqrt(2)-1
+		if absx < Small { // |x| < 2^-29
+			if absx < Tiny { // |x| < 2^-54
+				return x
+			}
+			return x - x*x*0.5
+		}
+		if x > Sqrt2HalfM1 { // Sqrt(2)/2-1 < x
+			// (Sqrt(2)/2-1) < x < (Sqrt(2)-1)
+			k = 0
+			f = x
+			iu = 1
+		}
+	}
+	var c float64
+	if k != 0 {
+		var u float64
+		if absx < Two53 { // 1<<53
+			u = 1.0 + x
+			iu = Float64bits(u)
+			k = int((iu >> 52) - 1023)
+			if k > 0 {
+				c = 1.0 - (u - x)
+			} else {
+				c = x - (u - 1.0) // correction term
+				c /= u
+			}
+		} else {
+			u = x
+			iu = Float64bits(u)
+			k = int((iu >> 52) - 1023)
+			c = 0
+		}
+		iu &= 0x000fffffffffffff
+		if iu < 0x0006a09e667f3bcd { // mantissa of Sqrt(2)
+			u = Float64frombits(iu | 0x3ff0000000000000) // normalize u
+		} else {
+			k += 1
+			u = Float64frombits(iu | 0x3fe0000000000000) // normalize u/2
+			iu = (0x0010000000000000 - iu) >> 2
+		}
+		f = u - 1.0 // Sqrt(2)/2 < u < Sqrt(2)
+	}
+	hfsq := 0.5 * f * f
+	var s, R, z float64
+	if iu == 0 { // |f| < 2^-20
+		if f == 0 {
+			if k == 0 {
+				return 0
+			} else {
+				c += float64(k) * Ln2Lo
+				return float64(k)*Ln2Hi + c
+			}
+		}
+		R = hfsq * (1.0 - 0.66666666666666666*f) // avoid division
+		if k == 0 {
+			return f - R
+		}
+		return float64(k)*Ln2Hi - ((R - (float64(k)*Ln2Lo + c)) - f)
+	}
+	s = f / (2.0 + f)
+	z = s * s
+	R = z * (Lp1 + z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))))
+	if k == 0 {
+		return f - (hfsq - s*(hfsq+R))
+	}
+	return float64(k)*Ln2Hi - ((hfsq - (s*(hfsq+R) + (float64(k)*Ln2Lo + c))) - f)
+}