2010-02-01 23:21:40 -07:00
|
|
|
// 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
|
2016-03-01 16:21:55 -07:00
|
|
|
// and came with this notice. The go code is a simplified
|
2010-02-01 23:21:40 -07:00
|
|
|
// 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
|
2010-04-09 15:37:33 -06:00
|
|
|
// 1+x = 2**k * (1+f),
|
2010-02-01 23:21:40 -07:00
|
|
|
// 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)
|
2011-11-27 19:01:08 -07:00
|
|
|
// and
|
2010-02-01 23:21:40 -07:00
|
|
|
// | 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
|
2011-11-27 19:01:08 -07:00
|
|
|
// Log1p(±0) = ±0
|
2010-02-01 23:21:40 -07:00
|
|
|
// Log1p(-1) = -Inf
|
|
|
|
// Log1p(x < -1) = NaN
|
|
|
|
// Log1p(NaN) = NaN
|
math: regularize build
This will be nicer to the automatic tools.
It requires a few more assembly stubs
but fewer Go files.
There are a few instances where it looks like
there are new blobs of code, but they are just
being copied out of deleted files.
There is no new code here.
Suppose you have a portable implementation for Sin
and a 386-specific assembly one. The old way to
do this was to write three files
sin_decl.go
func Sin(x float64) float64 // declaration only
sin_386.s
assembly implementation
sin_port.go
func Sin(x float64) float64 { ... } // pure-Go impl
and then link in either sin_decl.go+sin_386.s or
just sin_port.go. The Makefile actually did the magic
of linking in only the _port.go files for those without
assembly and only the _decl.go files for those with
assembly, or at least some of that magic.
The biggest problem with this, beyond being hard
to explain to the build system, is that once you do
explain it to the build system, godoc knows which
of sin_port.go or sin_decl.go are involved on a given
architecture, and it (correctly) ignores the other.
That means you have to put identical doc comments
in both files.
The new approach, which is more like what we did
in the later packages math/big and sync/atomic,
is to have
sin.go
func Sin(x float64) float64 // decl only
func sin(x float64) float64 {...} // pure-Go impl
sin_386.s
// assembly for Sin (ignores sin)
sin_amd64.s
// assembly for Sin: jmp sin
sin_arm.s
// assembly for Sin: jmp sin
Once we abandon Makefiles we can put all the assembly
stubs in one source file, so the number of files will
actually go down.
Chris asked whether the branches cost anything.
Given that they are branching to pure-Go implementations
that are not typically known for their speed, the single
direct branch is not going to be noticeable. That is,
it's on the slow path.
An alternative would have been to preserve the old
"only write assembly files when there's an implementation"
and still have just one copy of the declaration of Sin
(and thus one doc comment) by doing:
sin.go
func Sin(x float64) float64 { return sin(x) }
sin_decl.go
func sin(x float64) float64 // declaration only
sin_386.s
// assembly for sin
sin_port.go
func sin(x float64) float64 { portable code }
In this version everyone would link in sin.go and
then either sin_decl.go+sin_386.s or sin_port.go.
This has an extra function call on all paths, including
the "fast path" to get to assembly, and it triples the
number of Go files involved compared to what I did
in this CL. On the other hand you don't have to
write assembly stubs. After starting down this path
I decided that the assembly stubs were the easier
approach.
As for generating the assembly stubs on the fly, much
of the goal here is to eliminate magic from the build
process, so that zero-configuration tools like goinstall
or the new go tool can handle this package.
R=golang-dev, r, cw, iant, r
CC=golang-dev
https://golang.org/cl/5488057
2011-12-13 13:20:12 -07:00
|
|
|
func Log1p(x float64) float64
|
|
|
|
|
|
|
|
func log1p(x float64) float64 {
|
2010-02-01 23:21:40 -07:00
|
|
|
const (
|
|
|
|
Sqrt2M1 = 4.142135623730950488017e-01 // Sqrt(2)-1 = 0x3fda827999fcef34
|
|
|
|
Sqrt2HalfM1 = -2.928932188134524755992e-01 // Sqrt(2)/2-1 = 0xbfd2bec333018866
|
2010-04-09 15:37:33 -06:00
|
|
|
Small = 1.0 / (1 << 29) // 2**-29 = 0x3e20000000000000
|
|
|
|
Tiny = 1.0 / (1 << 54) // 2**-54
|
|
|
|
Two53 = 1 << 53 // 2**53
|
2010-02-01 23:21:40 -07:00
|
|
|
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
|
|
|
|
switch {
|
2012-02-01 08:08:31 -07:00
|
|
|
case x < -1 || IsNaN(x): // includes -Inf
|
2010-02-01 23:21:40 -07:00
|
|
|
return NaN()
|
|
|
|
case x == -1:
|
|
|
|
return Inf(-1)
|
2012-02-01 08:08:31 -07:00
|
|
|
case IsInf(x, 1):
|
2010-02-01 23:21:40 -07:00
|
|
|
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
|
2010-04-09 15:37:33 -06:00
|
|
|
if absx < Small { // |x| < 2**-29
|
|
|
|
if absx < Tiny { // |x| < 2**-54
|
2010-02-01 23:21:40 -07:00
|
|
|
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
|
2010-04-09 15:37:33 -06:00
|
|
|
if iu == 0 { // |f| < 2**-20
|
2010-02-01 23:21:40 -07:00
|
|
|
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)
|
|
|
|
}
|