I'm a little confused about how sbcl tests float equivalence.
And looking into the HS I did not see whether it is specified what
floating point equivalence means.
Are two floating point numbers equal if their binary bits are
identical?
Or if when printed in base 10, the digits are the same?
Or is equivalence allowed to ignore some guard digits?
( as in http://docs.sun.com/source/806-3568/ncg_goldberg.html )
The example below seems to indicate that they are equal if when
printed in base 10 the digits are the same.
I was expecting the following number not to be exactly zero.
(- (* 2 (/ 1 1000.0)) (/ 2 1000.0))
but it seems in fact to be.
(zerop (- (* 2 (/ 1 1000.0)) (/ 2 1000.0)))
returns T.
a*(1/1000.0) - (a/1000.0)
It is a bit more perplexing to me when i try this calculation with
numbers from 1 to 29. Somtimes the similar caluclation comes up
with zero, sometimes it does not. I'm supprised that it ever does
when (/ a b) has a non-terminating binary representation.
In the example where, zero occurs on 1, 2, 3, 4, 6, 7, 8, 12,13,14,
16,17,
24,25,26,27,and 28, but not with 5, 9, 10 ,11, 16,18, 19,20, 21,22,29
and 29.
FILL> (loop for denom from 1 below 30 do
(let ((b (/ denom 1000.0))
(a (/ 1 1000.0)))
(let ((c (* a denom)))
(format t "denom=~A a=~A b=~A c=~A c-b=~A 0? ~A~%"
denom a b c (- c b) (zerop (- c b))))))
denom=1 a=0.001 b=0.001 c=0.001 c-b=0.0 0? T
denom=2 a=0.001 b=0.002 c=0.002 c-b=0.0 0? T
denom=3 a=0.001 b=0.003 c=0.003 c-b=0.0 0? T
denom=4 a=0.001 b=0.004 c=0.004 c-b=0.0 0? T
denom=5 a=0.001 b=0.005 c=0.0050000004 c-b=4.656613e-10 0?
NIL
denom=6 a=0.001 b=0.006 c=0.006 c-b=0.0 0? T
denom=7 a=0.001 b=0.007 c=0.007 c-b=0.0 0? T
denom=8 a=0.001 b=0.008 c=0.008 c-b=0.0 0? T
denom=9 a=0.001 b=0.009 c=0.009000001 c-b=9.313226e-10 0? NIL
denom=10 a=0.001 b=0.01 c=0.010000001 c-b=9.313226e-10 0? NIL
denom=11 a=0.001 b=0.011 c=0.011000001 c-b=9.313226e-10 0?
NIL
denom=12 a=0.001 b=0.012 c=0.012 c-b=0.0 0? T
denom=13 a=0.001 b=0.013 c=0.013 c-b=0.0 0? T
denom=14 a=0.001 b=0.014 c=0.014 c-b=0.0 0? T
denom=15 a=0.001 b=0.015 c=0.015000001 c-b=9.313226e-10 0?
NIL
denom=16 a=0.001 b=0.016 c=0.016 c-b=0.0 0? T
denom=17 a=0.001 b=0.017 c=0.017 c-b=0.0 0? T
denom=18 a=0.001 b=0.018 c=0.018000001 c-b=1.8626451e-9 0?
NIL
denom=19 a=0.001 b=0.019 c=0.019000001 c-b=1.8626451e-9 0?
NIL
denom=20 a=0.001 b=0.02 c=0.020000001 c-b=1.8626451e-9 0? NIL
denom=21 a=0.001 b=0.021 c=0.021000002 c-b=1.8626451e-9 0?
NIL
denom=22 a=0.001 b=0.022 c=0.022000002 c-b=1.8626451e-9 0?
NIL
denom=23 a=0.001 b=0.023 c=0.023000002 c-b=1.8626451e-9 0?
NIL
denom=24 a=0.001 b=0.024 c=0.024 c-b=0.0 0? T
denom=25 a=0.001 b=0.025 c=0.025 c-b=0.0 0? T
denom=26 a=0.001 b=0.026 c=0.026 c-b=0.0 0? T
denom=27 a=0.001 b=0.027 c=0.027 c-b=0.0 0? T
denom=28 a=0.001 b=0.028 c=0.028 c-b=0.0 0? T
denom=29 a=0.001 b=0.029 c=0.029000001 c-b=1.8626451e-9 0?
NIL
NIL
FILL>
Continuing the experiment, I wanted to print the floating point numbers
in base 2. But after setting *print-base* to 2, I notice that
*print-base* only applies to rationals, not to floats :-(
Does anyone know how to print floats in base 2?
If anyone is interested, here is a function that will print a
non-negative floating
point number according to *print-base*. It happens to do something
quite useful
for rationals as well; it prints the equivalent radix representation
with the number of
digits of percision requested.
(defun print-float (float percision stream)
(if (null stream)
(let ((fstr (make-array '(0) :element-type 'character
:fill-pointer 0 :adjustable t)))
(with-output-to-string (s fstr)
(print-float float percision s))
fstr)
(multiple-value-bind (right left) (truncate float)
(format stream "~A." right)
(dotimes (ign percision)
(setf left (* *print-base* left))
(multiple-value-setq (right left) (truncate left))
(format stream (if (minusp right)
"[~A]" "~A")
right))
stream)))
Jimka wrote:
> Continuing the experiment, I wanted to print the floating point numbers
> in base 2. But after setting *print-base* to 2, I notice that
> *print-base* only applies to rationals, not to floats :-(
>
> Does anyone know how to print floats in base 2?
On 2006-08-26 14:09:34 +0100, "Jimka" <·····@rdrop.com> said:
> I think i see why these floating point results are happening.
> It seems simply to be a coincidence that all the floating point
> round offs cancel each other out in some cases but not in others.
> I've written a program that prints the numbers in binary, and when
> i look at the digits i see that sometimes the result is zero because
> of the accident of two different round-offs returning the same result.
There are a bunch of functions in CL which let you look at the bits of
floating point numbers. Look at DECODE-FLOAT & friends.
--tim
From: Alexander Schmolck
Subject: Re: floating point calculations, what does equal mean?
Date:
Message-ID: <yfswt8vodrj.fsf@oc.ex.ac.uk>
"Jimka" <·····@rdrop.com> writes:
> I'm a little confused about how sbcl tests float equivalence.
> And looking into the HS I did not see whether it is specified what
> floating point equivalence means.
>
> Are two floating point numbers equal if their binary bits are
> identical?
(= nan nan) tends to be false one ieee754 conforming systems, but other than
that I'd hope so.
> Or if when printed in base 10, the digits are the same?
Not necessarily. Printing floating point numbers tends to be underspecified
(CL and other languages don't mandate the best possible decimal approximation,
for example). It's also very unlikely that the base 10 representation of an
ieee float is internally used as an equality criterion in any sane programming
language.
> Or is equivalence allowed to ignore some guard digits?
> ( as in http://docs.sun.com/source/806-3568/ncg_goldberg.html )
> The example below seems to indicate that they are equal if when
> printed in base 10 the digits are the same.
No, I think your reasoning is confused here.
> I was expecting the following number not to be exactly zero.
> (- (* 2 (/ 1 1000.0)) (/ 2 1000.0))
Why not? (/ 1 1000.0) and (/ 2 1000.0) have the same significand -- generally
multiplying or dividing by 2 is unlikely to affect the significand.
Try playing with this -- it tries to print the machine representation of a
given float (careful: it's not properly tested and implementing single-floats
as well as denormalized floats is left as an excercise to the reader):
(defun float-to-hex (f)
(check-type f double-float)
(cond
;; NaNs -- must come first; loose info which one
((/= f f)
"fff8000000000000")
((> f most-positive-double-float)
"7ff0000000000000")
((< f most-negative-double-float)
"fff0000000000000")
(:else
(multiple-value-bind (f e s) (integer-decode-float f)
(format nil "~x" (logior (ash (* 1/2 (- 1 s)) 63)
(ash (+ e 52 1023) 52)
f))))))
'as
p.s.: If any of the slime people is reading this -- I think it would be rather
nice addition to show rationals also as floats and floats also as hex (as
above) in the same contexts in which ints are displayed in multiple bases.
From: Pascal Bourguignon
Subject: Re: floating point calculations, what does equal mean?
Date:
Message-ID: <87d5any3av.fsf@informatimago.com>
"Jimka" <·····@rdrop.com> writes:
> I'm a little confused about how sbcl tests float equivalence.
> And looking into the HS I did not see whether it is specified what
> floating point equivalence means.
>
> Are two floating point numbers equal if their binary bits are
> identical?
> Or if when printed in base 10, the digits are the same?
> Or is equivalence allowed to ignore some guard digits?
> ( as in http://docs.sun.com/source/806-3568/ncg_goldberg.html )
> The example below seems to indicate that they are equal if when
> printed in base 10 the digits are the same.
>
> I was expecting the following number not to be exactly zero.
> (- (* 2 (/ 1 1000.0)) (/ 2 1000.0))
> but it seems in fact to be.
If your computer uses IEEE-754 floating point numbers, then
multiplicating or dividing by 2 only change the exponent, not the
mantisa.
> (zerop (- (* 2 (/ 1 1000.0)) (/ 2 1000.0)))
> returns T.
I'm not surprized.
> a*(1/1000.0) - (a/1000.0)
>
> It is a bit more perplexing to me when i try this calculation with
> numbers from 1 to 29. Somtimes the similar caluclation comes up
> with zero, sometimes it does not. I'm supprised that it ever does
> when (/ a b) has a non-terminating binary representation.
>
> In the example where, zero occurs on 1, 2, 3, 4, 6, 7, 8, 12,13,14,
> 16,17,
> 24,25,26,27,and 28, but not with 5, 9, 10 ,11, 16,18, 19,20, 21,22,29
> and 29.
Sometimes you have numbers that have an exact representation in
IEEE-754, sometimes you don't.
What Every Computer Scientist Should Know About Floating-Point Arithmetic
http://docs-pdf.sun.com/800-7895/800-7895.pdf
http://portal.acm.org/citation.cfm?id=103163
http://focus.hut.fi/docs/WorkShop/common/ug/goldberg1.doc.html
--
__Pascal Bourguignon__
Pascal Bourguignon wrote:
> > In the example where, zero occurs on 1, 2, 3, 4, 6, 7, 8, 12,13,14,
> > 16,17,
> > 24,25,26,27,and 28, but not with 5, 9, 10 ,11, 16,18, 19,20, 21,22,29
> > and 29.
>
> Sometimes you have numbers that have an exact representation in
> IEEE-754, sometimes you don't.
>
Hi Pascal, i do not think the situation is only a result of exact
arithmetic. i do not think that 13/1000.0 has an exact representation;
however (zerop (- (/ 13 1000.0) (* 13 (/ 1 1000.0))))
evaluates to T.
I think it is just a coincidence that with the number 13 all the
floating point errors offset each other, but with 19 they dont.
(zerop (- (/ 19 1000.0) (* 19 (/ 1 1000.0))))
evaluates to NIL.
On 26 Aug 2006 08:27:00 -0700, "Jimka" <·····@rdrop.com>:
>I think it is just a coincidence that with the number 13 all the
>floating point errors offset each other, but with 19 they dont.
"Offset each other" is not exact IMHO.
With 13, the result of "-" is so tiny that it rounds down to 0.
With 19, it doesn't.
In a decimal 4-decimal floating-point system, you'd get:
1/3 + 2/3 = 0.3333 + 0.6667 = 1
1/3 + 1/3 + 1/3 = 0.3333 + 0.3333 + 0.3333 = 0.9999
You could consider floating-point calculation errors as a (rather bad)
random number generator:
Start with n = 3.
Compute (if (zerop (- (/ n 1000.0) (* n (/ 1 1000.0)))) 1 0)
Increment n by 2 to get the next random number.