math: math.Hypot returns incorrect result

[gopherbot]

by lukeusmaximus39:

Before filing a bug, please check whether it has been fixed since the
latest release. Search the issue tracker and check that you're running the
latest version of Go:

Run "go version" and compare against
http://golang.org/doc/devel/release.html  If a newer version of Go exists,
install it and retry what you did to reproduce the problem.

Thanks.

What does 'go version' print?
go version go1.3.3 linux/amd64

What steps reproduce the problem?
Compute the function math.Hypot(0.08879849107482636, -0.07597714093596328)

If possible, include a link to a program on play.golang.org.
http://play.golang.org/p/4jsrWchSln

What happened?
The value 0.11686615404799307 was returned

What should have happened instead?
It should have returned the correct value of 0.11686615404799305

Please provide any additional information below.

Java and Python both give the result 0.11686615404799305
I would suggest that this is a problem with square root accuracy/precision.

[ianlancetaylor]

Comment 1:

Labels changed: added repo-main, release-go1.5.

[josharian]

Comment 2:

Running this through a high precision calculator yields 0.1168661540479930566903, so the
Java/Python results do look more accurate. In the case of Python at least, it is calling
through to libc's hypot.
From a quick poke, it looks like when the arguments are close in magnitude (p <= 2q)
and reasonable in size, you get more accurate results from
return math.Sqrt((p-q)*(p-q) + 2*p*q)
than from what we currently use
q = q / p
return p * math.Sqrt(1+q*q)
Cf http://play.golang.org/p/yI4-ENaoB7
Input from someone with expertise would be welcomed.

Status changed to Accepted.

[ALTree]

Comment 3:

G. W. Stewart in "Matrix Algorithms: Volume 1", SIAM 1998, pages 139 and 144, cited by
B. Einarsson 
in "Accuracy and Reliability in Scientific Computing", SIAM 2005, page 48, suggest using
s = p + q
hypot(p, q) = s * math.Sqrt((p/s)*(p/s) + (q/s)*(q/s))
so, scaling with p+q instead of max{|p|, |q|}, because (I cite from the former)
"If the scale factor s in the algorithm [..] is replaced by max{|a|, |b|}, the results
may be less accurate on a hexadecimal machine. The reason is that the number
sqrt((a/s)^2 + (b/s)^2) [equivalent to our q = q / p;  return p * math.Sqrt(1+q*q); ndr]
is a little bit greater than one [iff min{a,b} is small? I think, ndr] so that the
leading three bits 
in its representation are zero."
Scaling with s = a + b returns the correct number in this case:
http://play.golang.org/p/P0GduYlgQ-
and also provides protection against [over|under]flow.
If we don't want to port the (rather complicated) libc hypot(), maybe the suggested
alternative 
expression would be a better "simple but good enough" implementation than the current
one.
This claim should be tested, though.

[ALTree]

I ported libc hypot to go and it seems to be accurate (tested on million random float64, 100% of the time the result was bit-for-bit equal to the one from glibc).

Unfortunatly the go code is slower (at least on x86_64) than the current (inaccurate) assembly implementation. And understably so: the guarantee that the error on the result will be below 1 ulp requires quite some work - more than 100 lines of C - while the current assembly implementation uses like a dozen of machine instructions.

Some quick benchmarks (go1.4 linux/amd64). Two small floats:

BenchmarkHypotOld   50000000            24.9 ns/op
BenchmarkHypotNew   20000000            86.3 ns/op

A small float and a huge one

BenchmarkHypotOld   50000000            24.9 ns/op
BenchmarkHypotNew   30000000            42.5 ns/op

Thoughts?
Too slow? Is it worth a CL?

[bradfitz]

Paging Mr. @griesemer ...

[cespare]

Isn't that Dr. @griesemer? ;)

[randall77]

Herr Griesemer?

On Wed, Jan 7, 2015 at 2:17 PM, Caleb Spare notifications@github.com
wrote:

Isn't that Dr. @griesemer https://github.com/griesemer? ;)

—
Reply to this email directly or view it on GitHub
#8909 (comment).

[minux]

Will there be any license issues with translating glibc code to Go and
include into Go standard packages?

[bradfitz]

Seems like yes, there would be. Unless those routines are dual-licensed under something that's not GPL.

[bradfitz]

Oh, glibc is LGPL.

[davecheney]

[insert dramatic chipmonk]
On 8 Jan 2015 10:11, "Brad Fitzpatrick" notifications@github.com wrote:

Oh, glibc is LGPL.

—
Reply to this email directly or view it on GitHub
#8909 (comment).

[ALTree]

Mmm, it's old Sun code, we already did that with Sqrt: http://golang.org/src/math/sqrt.go

[minux]

Then translate the real msun code, without looking at the glibc source.

https://raw.githubusercontent.com/freebsd/freebsd/master/lib/msun/src/e_hypot.c

[ALTree]

glibc is LGPL so it should be fine. Also the code in glibc is like 99.5% equal to the freebsd one.

[rsc]

We've never promised to get the last bit right.

We certainly CANNOT put LGPL code into this library.