f



Possible bug in Tcl or Windows or Tcl on Windows

Hi,

There seems to be a bug in the way numbers are compared in Tcl.
Consider the below script for calculating Pythagorean triplets.
For hypotenuse upto a value of 100, there should have been 63
unique triplets.

On Windows XP the script detects only 62. The script doesn't detect
the case where c=99, b=20 ==> a=101.

However running the same script under Tcl 8.4.1 in Cygwin detects
63 triplets.

I don't have a Linux machine at hand to test it there.

Following is the script and relevant output. Could anyone shed some
light
on the cause of this. Maybe it has something to do with how the
numbers
are represented internally?

Running the script for N>100 shows up many more such missed values.
An equivalent program in C runs correctly on the same machine. C code
was compiled using both gcc and VC++6.0.


#########################################################################
# a^2 = b^2 + c^2
proc pythag {MAX} {
	set i 0
	for {set c 2} {$c <= $MAX} {incr c} {
		for {set b 1} {$b < $c} {incr b} {
			set a [expr hypot($c, $b)]	;# Calc. Hypot

			if { ($c == 99) && ($b == 20)} { ;# <<<<<<<<
				puts ">> [expr round($a)] == $a"
			}

			if {[expr round($a)] == $a} {
				puts "$a : $b : $c"
				incr i
			}
		}
	}
	return $i
}

if {$argc == 1} {
	set MAX [lindex $argv 0]
} else {
	puts stderr "Usage: tclsh $argv0 N"
	exit
}

puts [pythag $MAX]

############# OUTPUT ################
Tcl 8.4.1 (Cygwin)
------------------------
5.0 : 3 : 4
10.0 : 6 : 8
<snip>
104.0 : 40 : 96
120.0 : 72 : 96
>> 101 == 101.0
101.0 : 20 : 99
125.0 : 75 : 100
63

Tcl 8.4.18 (Activestate Dist.)
---------------------
5.0 : 3 : 4
10.0 : 6 : 8
<snip>
104.0 : 40 : 96
120.0 : 72 : 96
>> 101 == 101.0
125.0 : 75 : 100
62

Tcl 8.5.2 (Activestate Dist.)
----------------------
5.0 : 3 : 4
10.0 : 6 : 8
<snip>
104.0 : 40 : 96
120.0 : 72 : 96
>> 101 == 100.99999999999999
125.0 : 75 : 100
62
###########################################


TIA
--
0
9/4/2009 4:58:19 PM
comp.lang.tcl 23429 articles. 2 followers. Post Follow

24 Replies
1345 Views

Similar Articles

[PageSpeed] 18

On Sep 4, 12:58=A0pm, solar <solaradmin2...@gmail.com> wrote:

>> 101 =3D=3D 100.99999999999999

From what I can tell, looking at tcl/generic/tclBasic.c and
tclExecute.c, along with tclWinPort.h , Tcl is calling what appears to
be the underlying C library hypot function directly. At least on Unix,
that function is defined to take 2 type double arguments and returns a
double return code.

Here's the peculiar thing. I modified your example to ensure that the
2 arguments to hypot were truncated to be an integer. So the values
passed to hypot should in my case anyways be exactly 20 and 99. The
underlying function, it would seem to me, should be calculating the
result as the square root of 10201. One would think that, even with
the fuzziness of floating point, that the answer to that question
would be 101. Apparently not, however.

http://wiki.tcl.tk/879 is one of I suspect several wiki pages that
discusses the dilemma of dealing with real numbers.



0
lvirden (1938)
9/4/2009 5:41:39 PM
Larry W. Virden wrote:
> On Sep 4, 12:58 pm, solar <solaradmin2...@gmail.com> wrote:
> 
>>> 101 == 100.99999999999999
> 
> From what I can tell, looking at tcl/generic/tclBasic.c and
> tclExecute.c, along with tclWinPort.h , Tcl is calling what appears to
> be the underlying C library hypot function directly. At least on Unix,
> that function is defined to take 2 type double arguments and returns a
> double return code.
> 
> Here's the peculiar thing. I modified your example to ensure that the
> 2 arguments to hypot were truncated to be an integer. So the values
> passed to hypot should in my case anyways be exactly 20 and 99. The
> underlying function, it would seem to me, should be calculating the
> result as the square root of 10201. 

doesn't matter if the inputs are integers - the result is a double

>  One would think that, even with
> the fuzziness of floating point, that the answer to that question
> would be 101. Apparently not, however.

why do you think 101 would have an exact representation in binary?

> 
> http://wiki.tcl.tk/879 is one of I suspect several wiki pages that
> discusses the dilemma of dealing with real numbers.
> 
> 
> 

yes the may be variances in behavior accross versions of tcl and various 
OS's but the only "bug" I see is in the code of the script itself that
is doing an exact comparison of floating point numbers.

a better check would be to get the rounded result - cast to integer
and then square that integer and compare it to the sum of squares of
the input. (since you end up needing the sum of squares anyway, you
could calc that then just use sqrt (then round, then cast, then square)
instead of using hypot.

or you could just change you comparison to use a suitable epsilon.

Bruce
0
9/4/2009 6:15:20 PM
solar wrote:
> Hi,
> 
> There seems to be a bug in the way numbers are compared in Tcl.
> Consider the below script for calculating Pythagorean triplets.
> For hypotenuse upto a value of 100, there should have been 63
> unique triplets.

It appears that the MSVC library generates a result with an error
of 1 unit in the least significant place for hypot(99.0, 20.0).
There is very little that Tcl can do about bugs in the underlying
C library.

I'm getting seriously tired of reimplementing C library functions
because the vendors can't be troubled to make them actually work.

By the way, you *do* realize that there are much better ways to
generate Pythagorean triples?

If m and n are relatively prime, m > n, and exactly one of m and
n is even, then

a = mn, b = m**2 - n**2, c = m**2 + n**2

is a primitive Pythagorean triple, and all three sides can be
multiplied by any integer to yield more Pythagorean triples.
No square rooting or other floating point machination needed.

-- 
73 de ke9tv/2, Kevin
0
kennykb (564)
9/4/2009 6:40:23 PM
At Fri, 4 Sep 2009 09:58:19 -0700 (PDT) solar <solaradmin2009@gmail.com> wrote:

> 
> Hi,
> 
> There seems to be a bug in the way numbers are compared in Tcl.
> Consider the below script for calculating Pythagorean triplets.
> For hypotenuse upto a value of 100, there should have been 63
> unique triplets.
> 
> On Windows XP the script detects only 62. The script doesn't detect
> the case where c=99, b=20 ==> a=101.
> 
> However running the same script under Tcl 8.4.1 in Cygwin detects
> 63 triplets.
> 
> I don't have a Linux machine at hand to test it there.
> 
> Following is the script and relevant output. Could anyone shed some
> light
> on the cause of this. Maybe it has something to do with how the
> numbers
> are represented internally?

Not likely.

I get 63 both with native linux (tcl_patchLevel = 8.4.7) and using a
Tclkit.exe running under Wine (tcl_patchLevel = 8.4.12).

> 
> Running the script for N>100 shows up many more such missed values.
> An equivalent program in C runs correctly on the same machine. C code
> was compiled using both gcc and VC++6.0.
> 
> 
> #########################################################################
> # a^2 = b^2 + c^2
> proc pythag {MAX} {
> 	set i 0
> 	for {set c 2} {$c <= $MAX} {incr c} {
> 		for {set b 1} {$b < $c} {incr b} {
> 			set a [expr hypot($c, $b)]	;# Calc. Hypot
> 
> 			if { ($c == 99) && ($b == 20)} { ;# <<<<<<<<
> 				puts ">> [expr round($a)] == $a"
> 			}
> 
> 			if {[expr round($a)] == $a} {
				*I* would suspect this line.
				You are comparing with *equality* a
				floating point number. That is often bad
				news.  Remember, the the *computer* is using
				finite presision *binary* floating point numbers
				the results are not *mathematically*
				what you might expect...  All it would
				take would be some guard bit jitter to
				completely ruin your day... Maybe you
				want (I'm suspecting that 101.0 is
				really 101.00000000001, which is != 101:
			if {[expr round($a)] == [expr {int(a)}]} {
> 				puts "$a : $b : $c"
> 				incr i
> 			}
> 		}
> 	}
> 	return $i
> }
> 
> if {$argc == 1} {
> 	set MAX [lindex $argv 0]
> } else {
> 	puts stderr "Usage: tclsh $argv0 N"
> 	exit
> }
> 
> puts [pythag $MAX]
> 
> ############# OUTPUT ################
> Tcl 8.4.1 (Cygwin)
> ------------------------
> 5.0 : 3 : 4
> 10.0 : 6 : 8
> <snip>
> 104.0 : 40 : 96
> 120.0 : 72 : 96
> >> 101 == 101.0
> 101.0 : 20 : 99
> 125.0 : 75 : 100
> 63
> 
> Tcl 8.4.18 (Activestate Dist.)
> ---------------------
> 5.0 : 3 : 4
> 10.0 : 6 : 8
> <snip>
> 104.0 : 40 : 96
> 120.0 : 72 : 96
> >> 101 == 101.0
> 125.0 : 75 : 100
> 62
> 
> Tcl 8.5.2 (Activestate Dist.)
> ----------------------
> 5.0 : 3 : 4
> 10.0 : 6 : 8
> <snip>
> 104.0 : 40 : 96
> 120.0 : 72 : 96
> >> 101 == 100.99999999999999
> 125.0 : 75 : 100
> 62
> ###########################################
> 
> 
> TIA
> --
>                                                                      

-- 
Robert Heller             -- 978-544-6933
Deepwoods Software        -- Download the Model Railroad System
http://www.deepsoft.com/  -- Binaries for Linux and MS-Windows
heller@deepsoft.com       -- http://www.deepsoft.com/ModelRailroadSystem/
                                                  
0
heller (3031)
9/4/2009 7:03:51 PM
Bruce wrote:
> ....
> 
> why do you think 101 would have an exact representation in binary?

Because it is 1100101 ?

0
9/4/2009 8:29:09 PM
On Sep 4, 11:15 am, Bruce <Bruce_do_not_...@example.com> wrote:

> why do you think 101 would have an exact representation in binary?

Uhhh... because 101 = 1+4+32+64

In fact all integers up to the fp mantissa size are
exactly represented as fp numbers.

Throw in rational numbers whose fractions bear the pattern
of division by two, like 123.625 or even 48.828125 and you
have lots more numbers exactly representable.

Donald Arseneau


0
asnd (4601)
9/4/2009 8:37:19 PM
At Fri, 04 Sep 2009 13:29:09 -0700 Gerry Snyder <mesmerizerfan@gmail.com> wrote:

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

> 
>                                                            

-- 
Robert Heller             -- 978-544-6933
Deepwoods Software        -- Download the Model Railroad System
http://www.deepsoft.com/  -- Binaries for Linux and MS-Windows
heller@deepsoft.com       -- http://www.deepsoft.com/ModelRailroadSystem/
                                            
0
heller (3031)
9/5/2009 12:12:23 AM
On Sep 4, 5:12 pm, Robert Heller <hel...@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?

DA
0
asnd (4601)
9/5/2009 8:43:14 AM
On Sep 5, 6:46=A0am, solar <solaradmin2...@gmail.com> wrote:
>
> > By the way, you *do* realize that there are much better ways to
> > generate Pythagorean triples?
>
> Yes. I just found some code comparing the exectution time of triplets
> in php, erlang, ruby etc,. and rewrote the same in Tcl. It was only
> for a comparison.

Uh, so you're saying that the brilliant algorithm you're benchmarking
across languages uses the floating-point test

       A=3Dhypoth(B,C)
       round(A)=3D=3DA    (whatever '=3D=3D' means here)

instead of the accurate, integer-based test

       A=3Dround(hypoth(B,C))
       A*A=3D=3DB*B+C*C

as Bruce suggested ?

In a sense you're fortunate that Tcl helped you stumble on that OS-
specific hypot() inaccuracy...

-Alex
0
9/5/2009 8:51:39 AM
At Sat, 5 Sep 2009 01:43:14 -0700 (PDT) Donald Arseneau <asnd@triumf.ca> wrote:

> 
> On Sep 4, 5:12 pm, Robert Heller <hel...@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.  I.e. the match function
is not return 101, it is returning 100.9999999999999 or something like
that. Assuming that although *mathematically* the result should be an
integer result does not mean that the math library will in fact return
the exact integer result.  As soon as you start doing arithmetic on a
computer *using finite precision floating point*, you start picking up
errors and over time (number of calculations) you will pickup errors due
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.

> 
> DA
>                                   

-- 
Robert Heller             -- 978-544-6933
Deepwoods Software        -- Download the Model Railroad System
http://www.deepsoft.com/  -- Binaries for Linux and MS-Windows
heller@deepsoft.com       -- http://www.deepsoft.com/ModelRailroadSystem/
                                                                                                        
0
heller (3031)
9/5/2009 11:43:54 AM
On Sep 5, 1:51=A0pm, Alexandre Ferrieux <alexandre.ferri...@gmail.com>
wrote:
> On Sep 5, 6:46=A0am, solar <solaradmin2...@gmail.com> wrote:
>
>
>
> > > By the way, you *do* realize that there are much better ways to
> > > generate Pythagorean triples?
>
> > Yes. I just found some code comparing the exectution time of triplets
> > in php, erlang, ruby etc,. and rewrote the same in Tcl. It was only
> > for a comparison.
>
> Uh, so you're saying that the brilliant algorithm you're benchmarking
> across languages uses the floating-point test
>
> =A0 =A0 =A0 =A0A=3Dhypoth(B,C)
> =A0 =A0 =A0 =A0round(A)=3D=3DA =A0 =A0(whatever '=3D=3D' means here)
>
> instead of the accurate, integer-based test
>
> =A0 =A0 =A0 =A0A=3Dround(hypoth(B,C))
> =A0 =A0 =A0 =A0A*A=3D=3DB*B+C*C
>
> as Bruce suggested ?

Brilliant???. I don't know. I not qualified to comment on that.
See
http://tempe.st/2007/05/erlang-ruby-and-php-battle-it-out/
and
http://ocamlnews.blogspot.com/2007/05/ocaml-whups-lisp-on-pythagorean.html


>
> In a sense you're fortunate that Tcl helped you stumble on that OS-
> specific hypot() inaccuracy...

May be that's a good thing ;-)

Thanks,
--
0
9/5/2009 12:38:24 PM
Robert Heller wrote:
> At Sat, 5 Sep 2009 01:43:14 -0700 (PDT) Donald Arseneau <asnd@triumf.ca> wrote:
> 
>>
>> 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.  I.e. the match function
> is not return 101, it is returning 100.9999999999999 or something like
> that. Assuming that although *mathematically* the result should be an
> integer result does not mean that the math library will in fact return
> the exact integer result.  As soon as you start doing arithmetic on a
> computer *using finite precision floating point*, you start picking up
> errors and over time (number of calculations) you will pickup errors due
> 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.

Errr,...

A math library that says 15 / 3 = 4.999999 or sqrt(10201) = 100.9999999 
has an error that should be fixed.


Gerry
0
9/5/2009 3:11:48 PM
Gerry Snyder wrote:
> Robert Heller wrote:
>> At Sat, 5 Sep 2009 01:43:14 -0700 (PDT) Donald Arseneau 
>> <asnd@triumf.ca> wrote:
>>
>>>
>>> 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.  I.e. the match function
>> is not return 101, it is returning 100.9999999999999 or something like
>> that. Assuming that although *mathematically* the result should be an
>> integer result does not mean that the math library will in fact return
>> the exact integer result.  As soon as you start doing arithmetic on a
>> computer *using finite precision floating point*, you start picking up
>> errors and over time (number of calculations) you will pickup errors due
>> 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.
> 
> Errr,...
> 
> A math library that says 15 / 3 = 4.999999 or sqrt(10201) = 100.9999999 
> has an error that should be fixed.
> 
> 
> Gerry

15 / 3 = 5

but

15.0 / 3.0 may equal 4.99999...  since it is perfectly obvious that 4.99... 
and 5.0 the exact same number.

Note in fixed length floating point operations, round off errors do happen 
(normally at the hardware level but also at the software level). To expect 
anything less is to not live in the real world.

-- 
+------------------------------------------------------------------------+
| Gerald W. Lester                                                       |
|"The man who fights for his ideals is the man who is alive." - Cervantes|
+------------------------------------------------------------------------+
0
Gerald.Lester (2014)
9/5/2009 3:49:37 PM
* solar <solaradmin2009@gmail.com>
| --> if {[expr $a-int($a)] < 0.00001} { incr i }
| works for atleast N upto 15000, but its extremely slow going. :-(

Try

   if { ($a-int($a)) < 0.00001} { incr i }
   
should be a tick faster.

    proc foo {a} {set i 0 ; if {[expr $a-int($a)] < 0.00001} { incr i }}
    time {foo 1} 100000
    14.36562 microseconds per iteration

    proc foo {a} {set i 0 ; if {($a-int($a)) < 0.00001} { incr i }}
    time {foo 1} 100000
    1.88318 microseconds per iteration

R'
0
ralfixx (1283)
9/5/2009 6:39:27 PM
At Sat, 05 Sep 2009 08:11:48 -0700 Gerry Snyder <mesmerizerfan@gmail.com> wrote:

> 
> Robert Heller wrote:
> > At Sat, 5 Sep 2009 01:43:14 -0700 (PDT) Donald Arseneau <asnd@triumf.ca> wrote:
> > 
> >>
> >> 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.  I.e. the match function
> > is not return 101, it is returning 100.9999999999999 or something like
> > that. Assuming that although *mathematically* the result should be an
> > integer result does not mean that the math library will in fact return
> > the exact integer result.  As soon as you start doing arithmetic on a
> > computer *using finite precision floating point*, you start picking up
> > errors and over time (number of calculations) you will pickup errors due
> > 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.
> 
> 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

% set x 14.999997
14.999997
% set y 3
3
% expr {$x / $y}
4.999999
% format %.0f $x          
15
	

> 
> 
> Gerry
>                                                             

-- 
Robert Heller             -- 978-544-6933
Deepwoods Software        -- Download the Model Railroad System
http://www.deepsoft.com/  -- Binaries for Linux and MS-Windows
heller@deepsoft.com       -- http://www.deepsoft.com/ModelRailroadSystem/
                                                                                        
0
heller (3031)
9/5/2009 8:09:14 PM
Kevin Kenny wrote:
> 
> It appears that the MSVC library generates a result with an error
> of 1 unit in the least significant place for hypot(99.0, 20.0).
> There is very little that Tcl can do about bugs in the underlying
> C library.

And following up on this, it appears that the na�ve approach
of [expr {sqrt($a*$a + $b*$b)}] works better than [expr {hypot($a,$b)}]
for this case. (The chief advantage of hypot() is that it at least
supposedly does not overflow even if $a**2 or $b**2 does.)

-- 
73 de ke9tv/2, Kevin
0
kennykb (564)
9/6/2009 4:08:36 AM
On Sep 5, 11:39=A0pm, Ralf Fassel <ralf...@gmx.de> wrote:
> * solar <solaradmin2...@gmail.com>
> | --> if {[expr $a-int($a)] < 0.00001} { incr i }
> | works for atleast N upto 15000, but its extremely slow going. :-(
>
> Try
>
> =A0 =A0if { ($a-int($a)) < 0.00001} { incr i }
>
> should be a tick faster.
>
> =A0 =A0 proc foo {a} {set i 0 ; if {[expr $a-int($a)] < 0.00001} { incr i=
 }}
> =A0 =A0 time {foo 1} 100000
> =A0 =A0 14.36562 microseconds per iteration
>
> =A0 =A0 proc foo {a} {set i 0 ; if {($a-int($a)) < 0.00001} { incr i }}
> =A0 =A0 time {foo 1} 100000
> =A0 =A0 1.88318 microseconds per iteration
>
> R'

Thanks Ralf,

Speaking of improving speed, I thought that maybe using Tcl 8.5 would
help. But the "very original pythag" above takes about 2.3 times
longer
to execute on 8.5.2.
---------------------------------------
% info pa
8.4.18
% time {pythag 1000}
12705706 microseconds per iteration
% time {pythag 1000} 2
12709078.5 microseconds per iteration

% info pa
8.5.2
% time {pythag 1000}
29223915 microseconds per iteration
% time {pythag 1000} 2
29243033.5 microseconds per iteration
---------------------------------------

Now why would this be so? Even with your tip trying

---------------------------------------
proc foo1 {a} {set i 0 ; if {($a-int($a)) < 0.00001} { incr i }}
proc foo2 {a} {set i 0 ; if {[ expr $a-int($a)] < 0.00001} { incr i }}

% info pa
8.4.18
% time {foo1 1} 100000
1.08594 microseconds per iteration
% time {foo2 1} 100000
5.24368 microseconds per iteration
~~~~~~~~~~~
% info pa
8.5.2
% time {foo1 1} 100000
1.37046 microseconds per iteration
% time {foo2 1} 100000
6.48785 microseconds per iteration
---------------------------------------

Your tip would help speed it up a little bit though.

Thanks,
--
0
9/6/2009 4:35:21 AM
On Sep 6, 12:35 am, solar <solaradmin2...@gmail.com> wrote:

>
> % info pa
> 8.4.18
> % time {foo1 1} 100000
> 1.08594 microseconds per iteration
> % time {foo2 1} 100000
> 5.24368 microseconds per iteration
> ~~~~~~~~~~~
> % info pa
> 8.5.2
> % time {foo1 1} 100000
> 1.37046 microseconds per iteration
> % time {foo2 1} 100000
> 6.48785 microseconds per iteration
> ---------------------------------------

I tried the scripts above on Tcl 8.5.6 and 8.6b1.1 on SPARC Solaris8.
My times were significantly worse than what you saw:
 /var/tmp/Ac*6/bin/tclsh8.5 /var/tmp/timing.tcl
8.5.6
34.62197 microseconds per iteration
165.51024 microseconds per iteration
37.19967 microseconds per iteration
tclsh8.6 /var/tmp/timing.tcl
8.6b1.1
52.31166 microseconds per iteration
202.41822 microseconds per iteration
57.00032 microseconds per iteration

Note that the third timing is for:
proc foo3 {a} {set i 0 ; if {[ expr {$a-int($a) < 0.00001}]} { incr
i }}

0
lvirden (1938)
9/6/2009 6:35:27 AM
Robert Heller wrote:
> At Sat, 05 Sep 2009 08:11:48 -0700 Gerry Snyder <mesmerizerfan@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
0
9/6/2009 11:53:45 PM
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
0
kennykb (564)
9/7/2009 5:30:39 PM
Try http://docs.sun.com/source/806-3568/ncg_goldberg.html for more
than you ever wanted to know about floating point numbers.

Jackson
0
9/8/2009 8:03:24 AM
On Sep 5, 4:43 am, Robert Heller <hel...@deepsoft.com> wrote:
> At Sat, 5 Sep 2009 01:43:14 -0700 (PDT) Donald Arseneau <a...@triumf.ca> wrote:

> > On Sep 4, 5:12 pm, Robert Heller <hel...@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
0
asnd (4601)
9/9/2009 3:30:34 AM
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
0
kennykb (564)
9/9/2009 12:58:28 PM
On Sep 9, 5:58 am, Kevin Kenny <kenn...@acm.org> wrote:
> 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,

Then the fp hardware uses 80 bits precision internally, so the RTL
has
less excuse for the imprecise answer with 64 bits.

> As I said before, hypot(a,b) isn't a simple sqrt(a*a+b*b)

Yes, I know.  Yet even the naive implementation should give
an exact answer for 20,99 if it is using extra precision
internally.

And how am I getting sucked onto the other side?  I just dropped in
to dispute the implication that 101 was not an exact number in fp
representation.

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

Why the fuss about last-bit errors?  Certainly this example isn't
bad enough for Tcl to carry its own reimplementation!  People need
to have the habit of using fuzzy comparisons for equality as long
as there is some source of round-off errors, and there always will
be.

Donald Arseneau
0
asnd (4601)
9/10/2009 12:25:03 AM
Reply:

Similar Artilces:

Tcl + windows
Just quickly - When trying to run TCL programs under Windows (XP), if I leave the TCL_LIBRARY environment variable UNset, then it complains about not being able to find init.tcl, but if I set it - even if it is an incorrect path - it crashes (does the "send bug report?" thing). Thanks. gjb. "Gareth Bradley" <gb20@cs.waikato.ac.nz> wrote: >Just quickly - When trying to run TCL programs under Windows (XP), if I >leave the TCL_LIBRARY environment variable UNset, then it complains about >not being able to find init.tcl, but if I set it - even if it is an &g...

tcl on windows
Hi all, I just started learning tcl to do some small automation using expect and I am kind of struck with basic "Hello World" prog :). Can any one help me to resolve this? I am working on windows xp and when I use "puts" stmt, I wanted to write output to standard stdout ( command prompt). Instead what happens is it opens graphical window with name being the file name. #!/bin/sh # \ exec tclsh "$0" ${1+"$@"} puts stdout "Hello World" My task is to write a samll expect script that should return the o/p when called from another perl /python ...

running other tcl from a tcl
Is that possible? 1 - sequentially to call two tcl files, like: c:\> tclkit-8.4.7.exe first.tcl second.tcl 2 - from a tcl to call other one # first.tcl content below exec hello.tcl # end regards, mauro ps.: I'm not experience in Tcl. Mauro Silva wrote: > Is that possible? > > 1 - sequentially to call two tcl files, like: > > c:\> tclkit-8.4.7.exe first.tcl second.tcl > > 2 - from a tcl to call other one > > # first.tcl content below > > exec hello.tcl > > # end > > regards, > m...

tcl-gaul: Genetic Algorithms for Tcl. (Tcl package)
This is an announcement for a relatively new Tcl project: tcl-gaul Tcl-gaul is a Tcl extension for genetic/evolutionary algorithm processing.It relies on the GAUL library: http://gaul.sourceforge.net/ * A genetic algorithm (GA) is a search technique used in computing to find exact or approximate solutions to optimization and search problems. Genetic algorithms are categorized as global search heuristics. They are a particular class of evolutionary algorithms that use techniques inspired by evolutionary biology such as inheritance, mutation, selection, and crossover. For ...

TCL on windows?
hi what environment do you suggest me to start learning TCL on windows? thanks In article <mn.732e7d77c2d3be1c.77691@ginko.ginko>, ginko <ginko@ginko.ginko> wrote: >hi what environment do you suggest me to start learning TCL on windows? >thanks > > I'm not sure what you mean by your question. Possible answers: A. Read <URL: http://wiki.tcl.tk/15577 >; B. Install <URL: http://wiki.tcl.tk/1875 >; C. Read <URL: http://wiki.tcl.tk/1681 >; D. Purchase <URL: http://wiki.tcl.tk/3960 >; .... On Jul 14, 7:34 am, ginko <gi...@ginko.ginko&...

tcl-pam: PAM authentication for Tcl (Tcl package)
This is an announcement for a relatively new Tcl project: tcl-pam Tcl-pam is a Tcl interface to the PAM* service of Linux. It provides a Tcl package that allows Tcl scripts to use PAM to authenticate users and programs. It relies on linux-pam library: http://www.kernel.org/pub/linux/libs/pam/ * PAM (Pluggable Authentication Modules): A mechanism to integrate multiple low−level authentication schemes into a high−level application programming interface (API). This enables programs that rely on authentication to be written independently of the underlying authentication scheme. Platform:...

tcl-mq: POSIX Message Queues for Tcl. (Tcl package)
This is an announcement for a relatively new Tcl project: tcl-mp Tcl-mp is a Tcl interface to POSIX Message Queues*. It provides a Tcl package that allows scripts to create/open/close/unlink multiple parallel message queues, and to send/receive messages synchronously and asynchronously to/from them. * A POSIX message queue is an Inter-Process Communication mechanism available on Linux and some other POSIX-compliant operating systems. It allows to or more processes (or threads) to communicate under the same OS. The messages are buffered by the kernel, which gives them kernel persis...

tcl-mmap: A POSIX mmap interface for Tcl. (Tcl package)
This is an announcement for a relatively new Tcl project: tcl-mmap Tcl-mmap is a Tcl interface to the POSIX mmap* system call. It provides a Tcl package that allows Tcl scripts to: 1) Memory map files for improved access efficiency; 2) Share memory between related processes; 3) Easily implement cyclic persistent log files. * See the mmap(2) man page. Platform: Linux/Unix Home page: http://sourceforge.net/projects/tcl-mmap/ Man page: http://tcl-mmap.sourceforge.net/ Author: Alexandros Stergiakis On Sep 3, 11:48=A0am, Alexandros Stergiakis <alst...@gmail.com> wrote: > This is an ...

tcl-syslog: Unix system logging for Tcl (Tcl package)
This is an announcement for a relatively new Tcl project: tcl-syslog Tcl-syslog is a Tcl interface to the *nix syslog service. It provides a Tcl package that allows Tcl scripts to log messages to syslog. Platform: Linux/Unix Home page: http://sourceforge.net/projects/tcl-syslog/ Man page: http://tcl-syslog.sourceforge.net/ Author: Alexandros Stergiakis alsterg ...

Bug in Tcl fconfigure / TclHttpd (upload.tcl)
Hi there, had some trouble tracking down a segfault problem in tclhttpd. Has anybody experienced this before? Is there a bugfix floating around? (see below for details) Thanks Bernd Bug in Tcl fconfigure / TclHttpd (upload.tcl) --------------------------------------------- Author: bernd.worsch@@pawisda.de Date: 05-04-12 DESCRIPTION upfile.tcl from tclhttpd 3.5.1 segfaults in proc UploadDomain due to "fconfigure $sock -trans auto". This happens using Tcl 8.4.2 and the Problem remains when using Tcl 8.4.9. COMMENT I think this happens in thread mode only. SOLUTION Using "...

TCL/TK window with no window decoration, but with keyboard and mouse events
Hi all, I want to make a toplevel window that is shown without window decoration on a linux box, but I want it to receive all mouse and keyboard events. I am trying to build a GUI that looks a bit like an old DOS GUI: one full screen window without decoration, and with a menu on the top of the window that can be navigated both with the mouse and keyboard. I am trying the following piece of code, but if I type <Alt-F>, it does not open the file menu. Can anyone give me a hint on how to achieve this? Thanks a lot #!/usr/bin/wish # Hide the main window wm withdraw . # Make sure the main ...

windows tcl/tk bug does not appear under Wine
I posted about this four days ago without a response so I'm reposting in the hope that somebody understands what is going on. I sent a friend a Starpack so that he could test a program under MS Windows. When he ran it, he encountered the old bug in which assigning a hexadecimal number like 0xFF to a scale variable (not "scalar", but the variable named in the -variable option of a scale widget) fails because 0xFF is interpreted as non- numeric. I've encountered this before and fixed it (by using decimal), but evidently hadn't changed it here. The MS Windows runtime I used...

Bug in font handling under windows? [tcl 8.6]
Hi all, I have just tried the latest beta from ActiveState of 8.6, and I am having a peculiar behavior under vista. The fixed fonts in my application get larger each time I start my application :D The problem I have observed is what returned by font actual: I have created a font "TextFont" with size 10: > font actual TextFont -family Arial -size 10 -weight normal -slant roman -underline 0 -overstrike 0 > font configure TextFont -size 11 > font actual TextFont -family Arial -size 11 -weight normal -slant roman -underline 0 -overstrike 0 Everything seems normal her...

tcl-snmptools: SNMP v1/v2/v3 operations for Tcl. (Tcl package)
This is an announcement for a relatively new Tcl project: tcl-snmptools Tcl-snmptools is a Tcl interface to the Net-SNMP library which provides operations for the management of remote SNMP agents. It supports all the standard SNMP v1/v2/v3 operations: connect, close, get, set, getnext, walk, bulkget, bulkwalk, trap, translate and others. It is currently in a functional state, but more work and testing needs to be done. Home page: http://sourceforge.net/projects/tcl-snmptols/ Man page: http://tcl-snmptols.sourceforge.net/ Author: Alexandros Stergiakis alsterg ...

Web resources about - Possible bug in Tcl or Windows or Tcl on Windows - comp.lang.tcl

Price Drop: Rocket Video Cast for Chromecast: Best Browser to watch and stream movies to your TV
Rocket Video Cast for Chromecast: Best Browser to watch and stream movies to your TV 3.0 Device: iOS Universal Category: Productivity Price: ...

Trai rules in favour of net neutrality
Telecom regulator says no service provider can offer, charge discriminatory tariffs for data services on basis of content

Why older travelers are less likely to cancel travel plans over Zika virus
Travelers 60 years or older are less likely than young travelers to cancel their travel plans because of fear over the Zika virus, according ...

Hackers publish data on thousands of DHS and FBI workers
Following the publication of data belonging to about 10,000 Department of Homeland Security workers, hackers have now published the contact info ...

Newsonomics: The New York Times Restarts its New Product Model, in Spanish
If your image of Mexico comes only from election rhetoric, you might think of it as a land of drug-pushing, crime-committing illegals aiming ...

Slack, Honor, Gear VR, Stewart Butterfield, Bill Gurley win at the 9th Annual Crunchies
TechCrunch : Slack, Honor, Gear VR, Stewart Butterfield, Bill Gurley win at the 9th Annual Crunchies — And The Winners Of The 9th Annual Crunchies ...

Taiwan Continues its Search for Earthquake Victims on Lunar New Year
Taiwanese officials fear that more than 100 people remain trapped in the rubble of a collapsed high-rise building in the city of Tainan days ...

Judge again denies Texas' efforts to block Syrian refugees
Dallas Morning News Judge again denies Texas' efforts to block Syrian refugees Post-Bulletin AUSTIN, Texas (AP) — A federal judge Monday again ...

Cruise Ship Does A 180 After Getting Wrecked By Huge Storm
Over 4,500 passengers and 1,600 crew members were on their way from New York to Florida and the Bahamas when their ship got tossed around like ...

It’s official: Instagram for iOS now supports multiple account switching
Instagram has officially announced that it has begun rolling out support for switching between multiple accounts in its iOS app, without requiring ...

Resources last updated: 2/9/2016 7:56:22 AM