Closed gopherbot closed 9 years ago
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.
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.
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?
Paging Mr. @griesemer ...
Isn't that Dr. @griesemer? ;)
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 https://github.com/golang/go/issues/8909#issuecomment-69102666.
Will there be any license issues with translating glibc code to Go and include into Go standard packages?
Seems like yes, there would be. Unless those routines are dual-licensed under something that's not GPL.
Oh, glibc is LGPL.
[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 https://github.com/golang/go/issues/8909#issuecomment-69109626.
Mmm, it's old Sun code, we already did that with Sqrt: http://golang.org/src/math/sqrt.go
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
glibc is LGPL so it should be fine. Also the code in glibc is like 99.5% equal to the freebsd one.
We've never promised to get the last bit right.
We certainly CANNOT put LGPL code into this library.
by lukeusmaximus39: