math package accuracy

[minux]

Hi gophers,

There are at least three issues about the accuracy of our math package:

#8909 about Hypot, #9545 about Asinh and #9536 about Log.

I've made some investigation into the last issue about math.Log.

The reporter claimed that our math.Log is wrong for one particular input,

however, our result is still within 1ulp of the actual value.

http://play.golang.org/p/F5CW3_8tmv

The problem is just that our result is not the closest representable float64

value of the actual value. (i.e. our error is less than 1ulp as documented

in the code (http://golang.org/src/math/log.go#L63), but bigger than 0.5ulp.)

The question is, to what degree of accuracy do we want for our math

package? Is 1ulp enough? Or do we need 0.5ulp, or guarantee the closest

representable value?

As an example, I've tested popular open source implementation of log on

float64 in libm. Most libraries (*BSD, newlib) are based on Sun's fdlibm, as

our code, and that only guarantees <1ulp error.

Glibc uses "IBM Accurate Mathematical Library", see

https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/e_log.c;h=7a498e83,

and uses multiple precision in the last reduction step to guarantee the "correct

rounded (to nearest) value of log(x)"

The only other piece code that computes the correct result for log is based

on http://dl.acm.org/citation.cfm?id=98267.98294, and needs 2KB of lookup

table, but no multiple precision floats.

(All the code except the glibc version is in:

https://github.com/minux/gobugs/tree/master/issue9546_math.Log_accuracy)

so the question is:

1. are we content with <1ulp accuracy for our math package? If this is the

case, we can close all three issues with working as we've intended.

2. if not, are we willing to use multiple precision floats or 2KB lookup table

to make sure our math.Log always return the most correct answer?

Thanks.

[Ian Lance Taylor]

On Sun, Jan 25, 2015 at 6:37 PM, minux <mi...@golang.org> wrote:
>
> so the question is:
> 1. are we content with <1ulp accuracy for our math package? If this is the
> case, we can close all three issues with working as we've intended.
>
> 2. if not, are we willing to use multiple precision floats or 2KB lookup
> table
> to make sure our math.Log always return the most correct answer?

I suppose I would vote for multi-precision floats. I don't mind
having the current code, but given that somebody has written a more
precise algorithmsp I can't think of a good reason to avoid using it.

Ian

[Michael Jones]

I have a doubled precision math library in Go. Doubled in that it uses two doubles in tandem, with exponents adjusted to keep the fraction bits contiguous and also uses care to do perfect rounding. It could be used here and there as needed in math primitives.

--
You received this message because you are subscribed to the Google Groups "golang-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to golang-dev+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Michael T. Jones | Chief Technology Advocate, Google  | m...@google.com |  +1 650-335-5765

[Brendan Tracey]

I don’t have enough experience with very nitty-gritty numerical work to know how much a ULP or two matters. I defer to the judgement of others. However, I do want to say that I use real Go programs where math.Tanh is 10 - 20% of the total runtime, and I have written Matlab programs where Log was a similarly large fraction. I can’t say what the correct decision is, but I just wanted to point out that there is code that would be affected by significant changes in execution time of the core math routines.

[Michael Jones]

The doubled math libraries are about 4x the cost in terms of execution time.

[Charlie Dorian]

When writing many of the math routines, I tried to strike a balance between speed (which I felt was very important - I used assembly when I could get away with it) and reasonable accuracy. Russ established a relative accuracy of 4x10**-16 as veryclose, which I regarded as good enough. Sometimes the accuracy was relaxed (Acos, Cosh, Exp, Exp2, Gamma, J0, J1, Jn, Lgamma, Pow, Sinh, Y0, Y1, Yn, and large argument Cos, Sin and Tan), essentially where exponential functions are involved. Look in math/all_test.go for details of the accuracy tests.

In general, I would find examples of code with compatible licenses, convert it to go, compare the speeds, and submit the fastest. From time-to-time, I still search for code that is at least as fast, more accurate, and useable with the license.

To answer your questions, 1) if possible, 1ulp is great, but other functions are less accurate; 2) math.Log is really not the only problem - a usable source of better algorithms is needed.

[minux]

I should have mentioned that for amd64, at least we can use FPU

to calculate Log more accurately. (And I made a typo in the initial

post, the issue for math.Log should be #9546, not #9536.)

see https://github.com/golang/go/issues/9546#issuecomment-70945323

for working code.

My benchmark showed that it's ~50% slower than our SSE2 implementation.

benchmark        old ns/op     new ns/op     delta

BenchmarkLog     11.5          16.8          +46.09%

(my original reply on the issue tracker reported no difference, apparently

I did something wrong. Now I've double checked this result, and also

disassembled the binary to make sure.)

[Charlie Dorian]

Using FPU instructions would not fix the non-Intel implementation (i.e., the pure go code) of math.Log; math/log_amd64.s is just a straight assembly version of math/log.go  . However, allowing use of the FPU in amd64 mode would speed up the amd64 math.Mod and math.Remainder greatly!

[Robert Griesemer]

There are a lot of applications where utmost accuracy is not needed (e.g. graphics applications where the error won't be visible), and typically there performance is paramount. In those applications, often special-cases versions of the relevant algorithms are used.

For the std math library it would be nice to return a result with an error <1ulp (ideally <= 0.5ulp). I'd probably vote for lookup tables (presumably they are faster). Higher-precision float arithmetic in software is not cheap, even if specialized to say 64bit mantissa (vs 53), if that's enough.

Ideally we should have an expert in the field look into this, or somebody willing to put in the time to understand the existing algorithms in depth so that they can be analyzed.

- gri

[Michael Jones]

I am not that expert, but know enough to have one expert observation. If one can use a single-rounding, higher precision floating point multiply accumulate (FPMA, aka QMAF from IBM RIOS days) then a same or nearly always same-speed library could be built that had the 0.5 ULP property. The FPMA also allows retention of more precision in argument reduction, which helps on the back side in the final rounding.

[Charlie Dorian]

Issue #681 is "ponies: fused multiply add". Its most recent comment is "rsc self-assigned this on May 18, 2014".

[minux]

On Mon, Jan 26, 2015 at 2:12 AM, Robert Griesemer <g...@golang.org> wrote:

There are a lot of applications where utmost accuracy is not needed (e.g. graphics applications where the error won't be visible), and typically there performance is paramount. In those applications, often special-cases versions of the relevant algorithms are used.

Right. I think the question is really for Go's math package, where should we draw the line between

ultimate accuracy and adequate performance (and what's adequate?) 

For the std math library it would be nice to return a result with an error <1ulp (ideally <= 0.5ulp). I'd probably vote for lookup tables (presumably they are faster). Higher-precision float arithmetic in software is not cheap, even if specialized to say 64bit mantissa (vs 53), if that's enough.

I did a preliminary benchmark using GCC 5 -O3 (glibc is built using 4.8.4, but that doesn't matter,

because its log is written in assembly)

All test run log(0.6023203277884887), which presumably is not easy to compute, 500M times:

glibc: 8.78s, fdlibm: 14.32s, FreeBSD's log (optimized fdlibm): 5.80s, openlibm (the table lookup method): 11.39s

The conclusion is that we might have a 2x slow-down if we switched from fdlibm to the more precise

table lookup method, unless, for course, we use hand written assembly.

I will try to extract glibc's portable C version and report back the result. (This seems non-trivial, as glibc build

process is quite convoluted that I haven't fully understand yet.)

Ideally we should have an expert in the field look into this, or somebody willing to put in the time to understand the existing algorithms in depth so that they can be analyzed.

Is there any good reference material recommended on this topic? I find material on this very

sparse, perhaps people regarded this as too low-level.

[Ian Lance Taylor]

On Mon, Jan 26, 2015 at 6:14 PM, minux <mi...@golang.org> wrote:
>
>
> I did a preliminary benchmark using GCC 5 -O3 (glibc is built using 4.8.4,
> but that doesn't matter,
> because its log is written in assembly)
>
> All test run log(0.6023203277884887), which presumably is not easy to
> compute, 500M times:
>
> glibc: 8.78s, fdlibm: 14.32s, FreeBSD's log (optimized fdlibm): 5.80s,
> openlibm (the table lookup method): 11.39s
>
> The conclusion is that we might have a 2x slow-down if we switched from
> fdlibm to the more precise
> table lookup method, unless, for course, we use hand written assembly.

Does the slowdown happen for all cases, or just a few?

> I will try to extract glibc's portable C version and report back the result.
> (This seems non-trivial, as glibc build
> process is quite convoluted that I haven't fully understand yet.)

The key idea for glibc is that each system has a list of directories,
and the overall build looks for specific base file names. The first
file with the desired base name in the list is the one that is used.

Ian

[Ulrich Kunitz]

The book "IA-64 and Elementary Functions" from Peter Markstein discusses the implementation of elementary functions. It gives the mathematical background for the C implementations he is providing.

[Russ Cox]

In general, for me, the answer to whether there is a problem depends on how complex the solutions are. If there is an easy way to get the correctly-rounded answer, great. If not, then I think it's fine to keep using the off-by-1-ulp code until an easy way presents itself. If by multiple precision floats you mean the kind of "float64 pair" math that mtj referred to, then that might be fine. Certainly a handful of the fdlibm routines already take this approach. If you mean arbitrary precision then I think that's not worth doing. And 2KB tables are probably not worth doing either.

For these specific cases, if we're on par with fdlibm, I think that's a fine bar, and you're welcome to close the issues with a note to that effect.

Russ