Possible bug in Tcl or Windows or Tcl on Windows
in [TCL]
|
From: Gerry Snyder on 6 Sep 2009 19:53 Robert Heller wrote: > At Sat, 05 Sep 2009 08:11:48 -0700 Gerry Snyder <mesmerizerfan(a)gmail.com> wrote: > >> Errr,... >> >> A math library that says 15 / 3 = 4.999999 or sqrt(10201) = 100.9999999 >> has an error that should be fixed. > > You are assuming that the inputs are 15 and 3 (or 10201). If one does a > long string of fp calculations, it is quite possible to end up with: > > 14.999997 / 3 = 4.999999 > or > sqrt(10200.9999798000) = 100.9999999 The discussion was not about a long string of fp calculations. The original problem was based on squaring two integers and taking the square root of the sum. A very different situation. Yes, I am aware of the dangers of expecting integer results from floating point math. I was first bitten by them over 45 years ago. But here I think the 10201 should be exact and its root should be exact. Gerry
From: Kevin Kenny on 7 Sep 2009 13:30 Gerry Snyder wrote: > But here I think the 10201 should be exact and its root should be exact. What part of "the bug is in the Microsoft Visual C++ runtime" didn't you understand? The C hypot() function is a little bit more complicated than sqrt(a*a + b*b), since it's required to take precautions against overflow and underflow in the intermediate results. The usual calculation (assuming that a <= b) looks something like: abs(a) * sqrt(1.0 + (abs(b)/abs(a))**2) with code around it to avoid divisions by zero and underflows if a and b are too different in magnitude (including if a is zero, of course). The imprecision comes from the division of abs(b)/abs(a); this division can yield a number that cannot be represented exactly. This approach is taken in Ancient Unix: http://unix-tree.huihoo.org/V7/usr/src/libm/hypot.c.html and in the Netlib 'fn' package: http://www.netlib.org/fn/cabs.f A similar approach, with Newton-Raphson iteration for sqrt(1 + p) is due to Moler and Morrison: http://www.koders.com/c/fid7D3C8841ADC384A5F8DE0D081C88331E3909BF3A.aspx I suspect, without source diving, that the Microsoft implementation is similar, since the Ancient Unix, Netlib, and Moler implementations all get the same result that it does. The BSD math library has a better implementation, which is hardly surprising, since it came from Kahan's lab and was written at about the same time as Kahan was working on the IEEE floating point standard. Both BSD and glibc return correctly rounded results for our problem child. It is clear that the MSVC library (like the Ancient Unix, Netlib and Moler ones) exhibits less-than-optimal rounding behaviour. It's less clear whether this is a severe enough problem to justify having the Tcl maintainers take on the maintenance of an independent implementation of hypot. As I said, I'm getting tired of carrying around code that serves no other purpose than to work around platform deficiencies in the standard C libraries. -- 73 de ke9tv/2, Kevin
From: Jackson on 8 Sep 2009 04:03 Try http://docs.sun.com/source/806-3568/ncg_goldberg.html for more than you ever wanted to know about floating point numbers. Jackson
From: Donald Arseneau on 8 Sep 2009 23:30 On Sep 5, 4:43 am, Robert Heller <hel...(a)deepsoft.com> wrote: > At Sat, 5 Sep 2009 01:43:14 -0700 (PDT) Donald Arseneau <a...(a)triumf.ca> wrote: > > On Sep 4, 5:12 pm, Robert Heller <hel...(a)deepsoft.com> wrote: > > > > > > why do you think 101 would have an exact representation in binary? > > > > > Because it is 1100101 ? > > > > That is an integer. It is not a floating point binary fraction. > > > Let me turn it around then ... why do you think it doesn't? > > > (I *know* it does.) > > > Do you think integers are not a subset of fp numbers for > > some reason? > > Sure they are. It is just not always the case that a fp math library > function will in fact return an integer result. Then why didn't you say "Are you sure the number you have is exactly 101?" > I.e. the match function > is not return 101, it is returning 100.9999999999999 or something like > that. Yes, the demonstration was that the hypot function was returning 100.99999999999999 instead of 101. I don't necessarily hold the view that all built-in mathematical functions must be precise to the last bit (but it would sure be nice). However it is not helpful to blur the distinction between what results were calculated and what numbers are representable. > to guard bit jitter and round off, etc. One should never assume > mathenetically exact results and almost *never* should use an equality > test, even when mathematically called for. To an extent, what is "mathematically" called for doesn't matter. Given expr { 100.0+(1.0/3.0)+(1.0/3.0)+(1.0/3.0) } the fp result *should be* 100.99999999999999 instead of 101... 101 would be an *error*! I will join in encouraging fuzzy comparisons for floating point numbers, but am still disappointed that hypot is inexact. (Knowing the floating- point hardware in use would clarify the deficiency... most fp processors today carry extra precision to eliminate round-off error in the final answer.) Donald Arseneau
From: Kevin Kenny on 9 Sep 2009 08:58
Donald Arseneau wrote: > I will join in encouraging fuzzy comparisons for floating point numbers, > but am still disappointed that hypot is inexact. (Knowing the floating- > point hardware in use would clarify the deficiency... most fp processors > today carry extra precision to eliminate round-off error in the final > answer.) Reiterating what I've said before: The hardware is any x86. The library is MS Visual C++ runtime, any version up to VC2008. As I said before, hypot(a,b) isn't a simple sqrt(a*a+b*b), because it's required to take precautions against overflow of the intermediate results. A common na�ve implementation is hypot(a,b) = abs(b)*sqrt(1.0+(abs(a)/abs(b))**2) (assuming without loss of generality that 0 < a <= b.) It's easy to see where the roundoff problem comes from in the formula above. The MSVC library has this bug, as do the example hypot implementations in netlib and in Moler and Morrison. The BSD and glibc implementations use a somewhat more complex implementation but give errors that are thought to be <1 ulp in all cases. Tcl could get around the problem - WHICH IS IN THE MSVC RUNTIME, NOT IN TCL ITSELF (sorry for shouting, but this thread has gone on too long) - by carrying around its own implementation of hypot(). The decision on whether to do so involves balancing the value of having a bit-precise hypot() function with the cost of the additional work for the Tcl maintainers. (My read: Bit-precise floating point *conversions* are worthwhile, because they restored Everything Is A String. Bit-precise hypot() is less so. On the other hand, the cost to the maintainers of answering the ignorant postings in this thread is beginning to approach the amortized cost of maintaining our own hypot() in perpetuity.) -- 73 de ke9tv/2, Kevin |