float vs fraction (just playing around)

[mike.f...@gmail.com]

I thought it neat that clojure supports fractions (just like
Smalltalk).
user=> (/ 3)
1/3

I was sort of surprised by this.

user=> (+ (float (* (/ 2) (/ 3))) (float (* (/ 2) (/ 3))) )
0.33333334

(= (+ (float (* (/ 2) (/ 3))) (float (* (/ 2) (/ 3))) ) (/ 3) )
yields true

How is tracking these number internally?

And, why does it present the float version of 1/3,
as rounded up?

[Mark H.]

On Mar 5, 10:00 pm, "mike.farn...@gmail.com" <mike.farn...@gmail.com>
wrote:

> I was sort of surprised by this.
>
> user=> (+ (float (* (/ 2) (/ 3))) (float (* (/ 2) (/ 3))) )
> 0.33333334
>
>  (=  (+ (float (* (/ 2) (/ 3))) (float (* (/ 2) (/ 3))) ) (/ 3) )
> yields true
>
> How is tracking these number internally?
>
> And, why does it present the float version of 1/3,
> as rounded up?

This is a question about floating-point arithmetic, not about
Clojure. See

http://en.wikipedia.org/wiki/IEEE_754

and

http://docs.sun.com/source/806-3568/ncg_goldberg.html

In particular, note that Java uses binary representation for floating-
point numbers, which probably explains the "rounding up" observation.

mfh

[David Sletten]

On Mar 5, 2009, at 8:00 PM, mike.f...@gmail.com wrote:

>
> I thought it neat that clojure supports fractions (just like
> Smalltalk).
> user=> (/ 3)
> 1/3
>
> I was sort of surprised by this.
>
> user=> (+ (float (* (/ 2) (/ 3))) (float (* (/ 2) (/ 3))) )
> 0.33333334

It comes as no surprise that certain numbers cannot be represented by
a finite string of decimal digits. We all realize that 1/3 =
0.3333..., and the "..." part is critical. Take it away and the
equality goes away too. In other words, if that string of 3's stops,
then we don't have 1/3 anymore.

This situation is true inside the computer too. The floating-point
representation of 1/3 is binary 0.01010101010101... Whereas in
decimal we have 1/3 = 3*1/10 + 3*1/100 + 3*1/1000 + ..., in binary we
have 1/3 = 0*1/2 + 1*1/4 + 0*1/8 + 1*1/16 + ... If this sum stops we
once again lose equality. However, the computer only has finite
memory in which to represent a floating-point number, so the string
of "01"s does in fact stop. The computer cannot represent 1/3 exactly
in floating-point. Suppose for the sake of argument that the computer
only stores the equivalent of 4 decimal digits. So we get 1/3 ->
0.3333. But the true value is 0.3333 < 1/3 < 0.3334. In terms of the
binary representation, depending on how many bits are used the
decimal representation may be slightly closer to 0.3334 rather than
0.3333. You saw that by using the "float" function, which coerces the
Ratio 1/3 to a single-precision float:
(float 1/3) => 0.33333334

If you had coerced to a double-precision float you would have seen this:
(double 1/3) => 0.3333333333333333

Neither one of these is more "correct" per se (the double is more
precise of course). They simply reflect the limitations of different
finite representations of 1/3.

You can explore on your own by using Clojure's "rationalize" function
(which as far as I can tell mirrors Common Lisp's RATIONAL function
rather than its RATIONALIZE function):
(rationalize (float 1/3)) => 416666679084301/1250000000000000
(rationalize (double 1/3)) => 3333333333333333/10000000000000000

How close are these to the true values?
(- (rationalize (float 1/3)) 1/3) => 37252903/3750000000000000
(- 1/3 (rationalize (double 1/3))) => 1/30000000000000000

So the single-precision representation is slightly more than 1/3,
which is why you see the terminal 4.

Here's a little program to help you look inside 1/3. Just type
(dec2bin/display-one-third <PRECISION>) where <PRECISION> is either
float or double (you will need to widen your terminal window):
(ns dec2bin
(:use clojure.contrib.math))

(defn integer-to-binary [i]
(if (zero? i)
"0"
(loop [i i
result ""]
(if (zero? i)
(apply str result)
(recur (quot i 2)
(if (zero? (rem i 2))
(cons \0 result)
(cons \1 result)))) )))

(defn fraction-to-binary [x digits]
(loop [n digits
fraction x
result ""]
(if (zero? n)
(apply str (reverse result))
(let [whole (quot (* 2 fraction) 1)
fraction (rem (* 2 fraction) 1)]
(recur (dec n) fraction (cons (char (+ whole (int \0)))
result)))) ))

(defn decimal-to-binary
([x] (decimal-to-binary x 10))
([x digits]
(let [whole (integer-to-binary (quot x 1))
fraction (fraction-to-binary (rem x 1) digits)]
(format "%s.%s" whole fraction))))

(defn display-one-third [f]
(prn (format "%5s %-20s %-40s %s" "bits" "floating-point"
"fraction" "binary"))
(let [precision (if (= f float) 16 32)]
(dotimes [i precision]
(let [j (inc i)
sum (reduce + (map #(/ (expt 2 (* 2 %))) (range 1 (inc
j))))]
(prn (format "%5d %-20s %-40s %s" (* 2 j) (str (f sum)) (str
sum) (decimal-to-binary sum (* 2 j)))) ))))

>
> (= (+ (float (* (/ 2) (/ 3))) (float (* (/ 2) (/ 3))) ) (/ 3) )
> yields true
>
>

You have to be very careful comparing floating-point numbers. While
the above discussion of 1/3 should not be surprising, many people do
get surprised when they look at other numbers in floating-point. A
perfect example is 0.1. In decimal, this looks like a "nice" fraction
that terminates after 1 decimal place. However, binary fractions are
not represented using powers of 10: 1/10, 1/100, 1/1000, etc. We have
to represent 0.1 in terms of: 1/2, 1/4, 1/8, etc. The binary
representation of 0.1 is:
0.00011001100110011001100110011001100110011001100110011...
Again this is an infinite string represented in a finite amount of
memory, so 0.1 is not exact. That's why you get this behavior:
(loop [i 1
x 0.1]
(prn (format "%d %s %s" i (str x) (str (* i 0.1))))
(when-not (= i 10)
(recur (inc i) (+ x 0.1))))

"1 0.1 0.1"
"2 0.2 0.2"
"3 0.30000000000000004 0.30000000000000004"
"4 0.4 0.4"
"5 0.5 0.5"
"6 0.6 0.6000000000000001"
"7 0.7 0.7000000000000001"
"8 0.7999999999999999 0.8"
"9 0.8999999999999999 0.9"
"10 0.9999999999999999 1.0"

Apparently adding 0.1 6 times is different from 6 * 0.1:
(rationalize 0.6) => 3/5
(rationalize (* 6 0.1)) => 6000000000000001/10000000000000000
(rationalize (+ 0.1 0.1 0.1 0.1 0.1 0.1)) => 3/5

Even worse is the potential for infinite loops:
(defn trap [x]
(prn x)
(cond (= x 1) "How nice. 10 * 0.1 = 1"
(> x 1) "Whoops. Good thing I had an escape hatch."
:else (trap (+ x 0.1))))

(trap 0.1)
0.1
0.2
0.30000000000000004
0.4
0.5
0.6
0.7
0.7999999999999999
0.8999999999999999
0.9999999999999999
1.0999999999999999
"Whoops. Good thing I had an escape hatch."

If you've read this far you'll want to take a look at the links Mark
provided earlier.

Aloha,
David Sletten

[mike.f...@gmail.com]

Thanks for the very detailed explanation and all of the code examples.
I was a bit twiddler a long time ago.
(I enjoyed reading the opcodes in the mainframe dumps)

Lisp is a new language. And, I think it counts as the Language of
theYear,
since it is so different than Java.

I was playing with java.util.Date and range, by converting to number
with .getTime.
I couldn't figure out how to create a Date. I finally re-read the Date
constructor
and realized I had to coerce the number to a long. The compiler/
interpreter was
telling me the truth. There wasn't a Date constructor for number.

I was also playing with "with-precision" and "rounding".




On Mar 6, 9:06 pm, David Sletten <da...@bosatsu.net> wrote:

[Mark H.]

On Mar 6, 8:06 pm, David Sletten <da...@bosatsu.net> wrote:
> It comes as no surprise that certain numbers cannot be represented by  
> a finite string of decimal digits. We all realize that 1/3 =  
> 0.3333..., and the "..." part is critical. Take it away and the  
> equality goes away too. In other words, if that string of 3's stops,  
> then we don't have 1/3 anymore.

> ...

That was about the kindest, friendliest introduction to floating-point
arithmetic I've ever seen -- great job!

mfh