Do not trust Floating Point

Introduction

The essence of this article in one sentence is "Do not blindly trust any result when computing with floating point numbers".

This may sound strange, but you have to be aware all the time, that computations involving floating point numbers are always in danger of being inexact and sometimes completely (absurdly) wrong. It is a good idea (i.e. highly recommended) that you check the expected error, precision and value ranges when doing floating point math.

An Example

Rump's Royal Pain¹⁾

This example shows absurdly wrong results from a seemingly innocent floating point computation. In addition to generating a completely wrong result in double precision IEEE floats, it also does so when using longer precision (80bit or 128 bit).

Try to compute the following expression (in your favorite programming language):

333.75y6 + x2(11x2y2 - y6 - 121y4 - 2) + 5.5y8 + x/(2y)

with:

 x = 77617 and y = 33096.

in Smalltalk, this could be written as:

x := 77617.
y := 33096.

(333.75 * (y ** 6)) 
+ ((x ** 2)
   * ((11 * (x ** 2) * (y ** 2))
      - (y ** 6)
      - (121 * (y ** 4))
      - 2))
+ (5.5 * (y ** 8))
+ (x / (2 * y))

or in JavaScript as:

x = 77617;
y = 33096;

(333.75 * (y ** 6)) 
+ ((x ** 2)
   * ((11 * (x ** 2) * (y ** 2))
      - (y ** 6)
      - (121 * (y ** 4))
      - 2))
+ (5.5 * (y ** 8))
+ (x / (2 * y))

when evaluated, the result will be 1.18059162071741e+21. Even when using higher precision IEEE floats (by using "x := 77617q" and "y := 33096q"), we get wrong results.
This is completely incorrect. The correct (approximated) result is: -0.8273960599..

So not even the sign is correct - and the IEEE value is off by 21 orders of magnitudes!

In Smalltalk/X, you can compute this by using large precision floats; simply write:

x := 77617QL.
y := 33096QL.
...

to get:

-0.827396059946821368141165095479816291999033115784384819

It could also be computed exactly, by using fractions instead:

|x y|

x := 77617.
y := 33096.

((333 + (3/4)) * (y ** 6)) 
+ ((x ** 2)
   * ((11 * (x ** 2) * (y ** 2))
      - (y ** 6)
      - (121 * (y ** 4))
      - 2))
+ ((11/2) * (y ** 8))
+ (x / (2 * y))

By the way, in expecco (and in many other programming languages), the above expression delivers the following results as per precision used:

Precision		Result
IEEE Single	x := 77617f. y := 33096f.	-1.18059162071741e+21	WRONG
IEEE Double	x := 77617. y := 33096.	1.18059162071741e+21	WRONG
IEEE Extended	x := 77617q. y := 33096q.	576460752303423489.2	WRONG
IEEE Quadruple	x := 77617Q. y := 33096Q.	1.17260394005317863185883490452	WRONG
IEEE Octuple	x := 77617QO. y := 33096QO.	-0.8273960599468213681411650...	OK
QDouble	x := 77617QD. y := 33096QD.	-0.8273960599468213681327577...	ALMOST
LargeFloat	x := 77617QL. y := 33096QL.	-0.8273960599468213681411650...	OK
Fractions	x := 77617. y := 33096.	(-54767/66192) which asLargeFloat is -0.8273960599468213681411650954...	EXACT

1) see Rump's Pain in books.google

---

Copyright © 2014-2024 eXept Software AG