math: Log1p may not be quite correct

[omion]

What version of Go are you using (go version)?

$ go version
go version go1.11 windows/amd64

Does this issue reproduce with the latest release?

Presumably; the file math/log1p hasn't changed in 2 years.

What operating system and processor architecture are you using (go env)?

go env Output

$ go env
set GOARCH=amd64
set GOBIN=
set GOCACHE=C:\Users\Omion\AppData\Local\go-build
set GOEXE=.exe
set GOFLAGS=
set GOHOSTARCH=amd64
set GOHOSTOS=windows
set GOOS=windows
set GOPATH=D:\Documents\ownCloud\go
set GOPROXY=
set GORACE=
set GOROOT=E:\Programs\Go
set GOTMPDIR=
set GOTOOLDIR=E:\Programs\Go\pkg\tool\windows_amd64
set GCCGO=gccgo
set CC=gcc
set CXX=g++
set CGO_ENABLED=1
set GOMOD=
set CGO_CFLAGS=-g -O2
set CGO_CPPFLAGS=
set CGO_CXXFLAGS=-g -O2
set CGO_FFLAGS=-g -O2
set CGO_LDFLAGS=-g -O2
set PKG_CONFIG=pkg-config
set GOGCCFLAGS=-m64 -mthreads -fno-caret-diagnostics -Qunused-arguments -fmessage-length=0 -fdebug-prefix-map=C:\Users\Omion\AppData\Local\Temp\go-build811762018=/tmp/go-build -gno-record-gcc-switches

What did you do?

Use the math.Log1p function

What did you expect to see?

An error of less than 1 ULP, as the log1p.go file's "Accuracy:" comment states

What did you see instead?

Very large errors for a small set of inputs.

First of all, I would like to apologize for not creating a pull request - I don't really know how git works...

I have determined that the problem is in the file math/log1p.go. Lines 154 - 159 (as of this writing) are:

			if k > 0 {
				c = 1.0 - (u - x)
			} else {
				c = x - (u - 1.0) // correction term
				c /= u
			}

whereas they should be:

			if k > 0 {
				c = 1.0 - (u - x)
			} else {
				c = x - (u - 1.0) // correction term
			}
			c /= u

The comment at the top of the file states:

//      ... 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.

however the code only does the proper c/u step if k < 0 (that's when the input to logp1 is between -1 and -0.29ish). For inputs between (2^p)-1 and 2^p for p from 1 to 52, the variable c is not calculated correctly. The maximum absolute error I have found is 0.5.

The playground shows the result of this issue. The inputs listed are in order, yet the log1p of the middle one is 0.5 larger than the other two. The error is around 1%.

I have also found the file from FreeBSD that is mentioned in the log1p.go header. It contains the following C:

	        c  = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
		c /= u;

The line c /= u is run no matter the result of (k>0) is, unlike the Go code.

[odeke-em]

Thank you for filing this issue @omion and welcome to the Go project!

I'll kindly page some experts your way @griesemer @bmkessler.

[Konstantin8105]

Playground and results are:

0.00000000000000e+00
5.00000000000000e-01
0.00000000000000e+00

All results must be zero.
Check by scilab:

-->log1p(4503599627370495.0)-log(4503599627370495.0+1)
 ans  =
 
    0. 


-->log1p(4503599627370495.5)-log(4503599627370495.5+1)
 ans  =
 
    0. 


-->log1p(4503599627370496.0)-log(4503599627370496.0+1)
 ans  =
 
    0.  

In according to https://code.woboq.org/userspace/glibc/sysdeps/ieee754/flt-32/s_log1pf.c.html#72
see line 72:

c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0);
c /= u;

Implementation in that part of code is different.
Location of problem is correct.

[gopherbot]

Change https://golang.org/cl/156120 mentions this issue: math: fix ternary correction statement in Log1p