asserting the equality of double values

asserting the equality of double values

[Kamaraju S Kusumanchi]

I was browsing the source code of guile 1.8.2 and have couple of trivial
questions. I would be very grateful if you can answer the following
questions.

1) In test-suite/standalone/test-round.c (lines 90-96) we have


      /* 2^DBL_MANT_DIG-1
         In the past scm_c_round had incorrectly incremented this value, due
         to the way that x+0.5 would round upwards (in the usual default
         nearest-even mode on most systems).  */
      x = ldexp (1.0, DBL_MANT_DIG) - 1.0;
      assert (x == floor (x));      /* should be an integer already */

here ldexp and floor both return double values. Is it guaranteed that
asserting the equality of two double values will always work? When learning
C, I remember reading somewhere not to rely on such comparisons.

I also do not understand the comment /* should be an integer already */ at
the end of the line. AFAIK both ldexp, floor function return double values
and not integers. Can someone explain the comment?

I am interested in knowing this because guile 1.8.2 fails to build on a
Debian machine using alpha architecture. According to the build log of the
failure (
http://buildd.debian.org/fetch.cgi?&pkg=guile-1.8&ver=1.8.2%2B1-2&arch=alpha&stamp=1188100514&file=log )
it occurs around this part of the code.


2) The above code uses DBL_MANT_DIG macro. How can I easily find out where
this macro is defined? Is there any book or reference which discuss this
kind of macros. For now I did a grep in /usr/include and found that such a
macro is defined in /usr/include/c++/4.2/limits file. This brings up
another question. 

If you are using a macro defined in a standard c++ library header, shouldn't
the file be named test-round.cpp instead of test-round.c?

thanks
raju

-- 
Kamaraju S Kusumanchi
http://www.people.cornell.edu/pages/kk288/
http://malayamaarutham.blogspot.com/



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://lists.gnu.org/mailman/listinfo/guile-devel

Re: asserting the equality of double values

[Nelson H. F. Beebe]

Kamaraju S Kusumanchi asks today:

>> ...
>>       x = ldexp (1.0, DBL_MANT_DIG) - 1.0;
>>       assert (x == floor (x));      /* should be an integer already */
>>
>> here ldexp and floor both return double values. Is it guaranteed that
>> asserting the equality of two double values will always work? When learning
>> C, I remember reading somewhere not to rely on such comparisons.
>> ...

floor() returns a double that contains a whole number representing the
integer part of its argument.  Since a double is a 64-bit quantity on
most modern machines, and in the IEEE 754 floating-point system, the
64-bit format contains a 53-bit significand, the integer part of a
double is too big to store as a 32-bit int type in C.

The comparision x == floor(x) is perfectly okay.  Floating-point
comparisons are NOT fuzzy: they are exact.  If x contained a whole
number on entry, then floor(x) is bit-for-bit identical to it.

Regrettably, there is a lot of misinformation about floating-point
arithmetic promulgated by book authors and programs who lack
sufficient understanding of the subject.

If you do something like

	x = 1.0;
	y = 3.0;
	if (3.0 * (x/y) == 1.0) ...

you cannot expect the equality to be true everywhere, because 1/3 is
not exactly representable in a finite number of bits in a binary
arithmetic system.  This particular example evaluates to true (1) with
IEEE 754 arithmetic in default rounding, but not with other rounding
modes.

However, the point in the sample code at the start of this message is
quite different: it uses ldexp(u,v) to construct a number u * 2**v,
and since u = 1.0 here, the result is 2**v, an EXACTLY-REPRESENTABLE
number.  Since DBL_MANT_DIG = 53, the result is the number 2**53 =
9_007_199_254_740_992.  Subtracting one produces a significant of all
1-bits:  2**53 - 1 = +0x1.fffffffffffffp+52, and that is the
next-to-largest exactly representable whole number in this arithmetic
system.  Its floor is identical, so "x == floor()" is true.

On some historical floating-point designs, this code might not work
like it does in IEEE 754 arithmetic, because of different rounding
behavior.

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- University of Utah                    FAX: +1 801 581 4148                  -
- Department of Mathematics, 110 LCB    Internet e-mail: beebe@math.utah.edu  -
- 155 S 1400 E RM 233                       beebe@acm.org  beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://lists.gnu.org/mailman/listinfo/guile-devel

Re: asserting the equality of double values

[Stephen Compall]

[-- Attachment #1.1: Type: text/plain, Size: 1069 bytes --]

On Sat, 2007-09-15 at 13:53 -0400, Kamaraju S Kusumanchi wrote:
> 2) The above code uses DBL_MANT_DIG macro. How can I easily find out where
> this macro is defined? Is there any book or reference which discuss this
> kind of macros. For now I did a grep in /usr/include and found that such a
> macro is defined in /usr/include/c++/4.2/limits file. This brings up
> another question. 
> 
> If you are using a macro defined in a standard c++ library header, shouldn't
> the file be named test-round.cpp instead of test-round.c?

The GNU C Library Reference Manual is a great reference for the standard
C library.  According to its text, the macro is defined to be included
by float.h.  For details, see Info node (libc)Floating Point Parameters.

-- 
;;; Stephen Compall ** http://scompall.nocandysw.com/blog **
"Peta" is Greek for fifth; a petabyte is 10 to the fifth power, as
well as fifth in line after kilo, mega, giga, and tera.
  -- Lee Gomes, performing every Wednesday in his tech column
     "Portals" on page B1 of The Wall Street Journal

[-- Attachment #1.2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 189 bytes --]

[-- Attachment #2: Type: text/plain, Size: 143 bytes --]

_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://lists.gnu.org/mailman/listinfo/guile-devel

Re: asserting the equality of double values

[Kevin Ryde]

"Nelson H. F. Beebe" <beebe@math.utah.edu> writes:
>
> Subtracting one produces a significant of all
> 1-bits:  2**53 - 1 = +0x1.fffffffffffffp+52,

Yes.

> and that is the
> next-to-largest exactly representable whole number in this arithmetic
> system.  Its floor is identical, so "x == floor()" is true.

Well, it's supposed to be.  I think I put in that test.  If it's failing
on alpha someone will want to dig around and see how it went wrong.  It
ought to be pretty simple, all the values are integers, so no rounding
at all (and it hasn't even reached the actual scm_c_round).


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://lists.gnu.org/mailman/listinfo/guile-devel

Re: asserting the equality of double values

[Kamaraju S Kusumanchi]

Kevin Ryde wrote:

> "Nelson H. F. Beebe" <beebe@math.utah.edu> writes:
>>
>> Subtracting one produces a significant of all
>> 1-bits:  2**53 - 1 = +0x1.fffffffffffffp+52,
> 
> Yes.
> 
>> and that is the
>> next-to-largest exactly representable whole number in this arithmetic
>> system.  Its floor is identical, so "x == floor()" is true.
> 
> Well, it's supposed to be.  I think I put in that test.  If it's failing
> on alpha someone will want to dig around and see how it went wrong.  It
> ought to be pretty simple, all the values are integers, so no rounding
> at all (and it hasn't even reached the actual scm_c_round).
> 

For the record, I have reported this bug to the Debian maintainers of glibc 
http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=442568 .

raju

-- 
Kamaraju S Kusumanchi
http://www.people.cornell.edu/pages/kk288/
http://malayamaarutham.blogspot.com/



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://lists.gnu.org/mailman/listinfo/guile-devel

Re: asserting the equality of double values

[Kamaraju S Kusumanchi]

Stephen Compall wrote:

> On Sat, 2007-09-15 at 13:53 -0400, Kamaraju S Kusumanchi wrote:
>> 2) The above code uses DBL_MANT_DIG macro. How can I easily find out
>> where this macro is defined? Is there any book or reference which discuss
>> this kind of macros. For now I did a grep in /usr/include and found that
>> such a macro is defined in /usr/include/c++/4.2/limits file. This brings
>> up another question.
>> 
>> If you are using a macro defined in a standard c++ library header,
>> shouldn't the file be named test-round.cpp instead of test-round.c?
> 
> The GNU C Library Reference Manual is a great reference for the standard
> C library.  According to its text, the macro is defined to be included
> by float.h.  For details, see Info node (libc)Floating Point Parameters.
> 

Yes, It is indeed part of float.h. Thank you.

raju

-- 
Kamaraju S Kusumanchi
http://www.people.cornell.edu/pages/kk288/
http://malayamaarutham.blogspot.com/



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://lists.gnu.org/mailman/listinfo/guile-devel