A note for beginners - Coding with floating point numbers

[Dilawar Singh]

Some people trust computer to a degree that they will accept its computation as
correct when algorithm is right and program has no bug. Some mistrust it to an
extent that they consider all of its calculation wrong. Between these two
extreme there are million others whom computer is a platform for various
activities which it does quite OK.

Issues involved in computation are rather serious. Not only one can not
represent 10/3 in accurately as finite decimal representation, but one can't
also represent some finite decimal representation such as 0.1 correctly in binary
format [1]. See [2] for a detailed commentary on issues involved with
computation with real numbers.

These issues have led many to explore alternatives of computation [3]. Till
people come up with a nice model to overcome these fundamental limitations, we
can find ways to deal with them in our day to day activities.

On Windows machine, try to compute $\sqrt[2]{4} - 2$ or $2- \sqrt[2]{4}$ in its
inbuilt calculator. If you find its answer surprising and think of it as bug
like [this guy](http://www.youtube.com/watch?v=oBNyBhzu8zI) you should read [1].

Reference [1] explains why and when $x == x + 1.0$ can be true and when and why $x
== x$ is false when dealing with floating point numbers. Its a great reference
for electrical engineering students who wants to implement a processor with
floating point units in some HDL such as verilog, VHDL, SystemC or Bluespec.

When adding two floating point numbers, their relative values make a huge
difference and one should be prepared for surprises. One can start his research
into this by this [wiki
article](http://en.wikipedia.org/wiki/Loss_of_significance). Ordering of
addition matters in such cases. Adding float a into float b can be very
different from adding b into a. This particular case is highlighted using
Haskell.

Following program computes the function $f(x) = (2- x) - (2^3 -
\frac{x^3}{12})$. One can translate it into other language also. This also
serve as a beauty of Haskell syntax.

```haskell
-- This is type signature of function. It says that function calc takes an
-- argument of type Float and returns a Float
calc :: Float -> Float

-- This is the definition of function. Quite readable.
calc x = (x - (x ^ 3) / 12)

-- Haskell can infer the type by itself. So we can choose not to write a type
-- declaration. This function subtract f(x) from f(2).
calc2 x = (calc 2) - (calc x)

-- This function takes a list of Float and apply function calc2 on all of its
-- elements and add them up. In python, equivalent is
--
-- def calcList1(l):
-- newlist = []
-- for a in l:
-- newlist.append(calc2(a))
-- return (0.0 + sum(newlist))
--
calcList1 :: [Float] -> Float
calcList1 l = foldl (+) 0.0 (map calc2 l)

-- This is like the previous function but sum is done in reverse order. This is
-- not entirely correct. we are using foldr (+) and foldl (+) function to sum a
-- list. While the foldr adds first to the last, foldl adds from last to first.
-- Again this is not entirely correct.

calcList2 :: [Float] -> Float
calcList2 l = foldr (+) 0.0 (map calc2 l)

-- Here we are applying the same function on the element of list but adding them
-- in different order. The input step construct a list l where each element is
-- step size away from the previos one.
-- e.g. [0.0,0.5..2.0] is [0.0, 0.5, 1.0, 1.5, 2.0]

test1 :: Float -> Float
test1 step = (calcList1 l) - (calcList2 l)
where
l = [0.0,step..2.0]

```

Let's see now what happens when we run this script `filename.hs` in `ghci
filename.hs`


*Main> test1 0.1
9.536743e-7
*Main> test1 0.01
2.2888184e-5
*Main> test1 0.001
2.4414063e-4
*Main> test1 0.0001
-3.7109375e-2
*Main>

One expects 0.0 but this is not what is happening here. Life with floating point
numbers is not easy, and it can be terribly confusing. If one is writing
simulator, it pays to be aware of these limitation.

References:

1. [What every computer scientist should know about floating-point
arithmatic](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html)
2. [How are real number specified in computation](http://cstheory.stackexchange.com/questions/2058/how-are-real-numbers-specified-in-computation/2070#2070)
3. [Continuum Computing – The Next Step in the Theory of Computing](http://www.ee.iitb.ac.in/web/schedule/seminar/Talk_by_NarendraKarmarkar)

[Manish Goregaokar]

The nice thing about Mathematica  (besides being good with the functional paradigm and list-based coding) is that it does away with floating point and stores every real number symbolically. 35242252425234/352514135413134 is exactly that fraction, no more, no less. The same goes for `Sin[1]`, it is `Sin[1]` and no t0.8414709848078965 . There are functions that let one control when something gets evaluated to a numerical value -- otherwise it tries to keep everything in terms of integers and functions of integers and known entities like Pi/E. This is a common newbie trip-up, for example typing `Sin[1]+Cos[1]` will just spit out the same string, unless you pass it through `N[]`. But once you get used to it, it makes the software very powerful in you hands as you have greater  knowledge of exactly where the floating point errors will occur, and a better degree of control over their effect. Of course, numerical methods like `NIntegrate[]`, `NDSolve[]` return floats by default

This actually extends to expressions as well, you can give it (x-1)(x+1) and tell it to hold its form for a period of computation, and then release the hold later on whereupon it will simplify and substitute.

Pretty neat.

-Manish Goregaokar

--
--
The website for the club is http://wncc-iitb.org/
To post to this group, send email to wncc...@googlegroups.com
--- You received this message because you are subscribed to the Google Groups "Web and Coding Club IIT Bombay" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wncc_iitb+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

[Saket Choudhary]

On 8 February 2014 04:13, Manish Goregaokar <manis...@gmail.com> wrote:
> The nice thing about Mathematica (besides being good with the functional
> paradigm and list-based coding) is that it does away with floating point and
> stores every real number symbolically. 35242252425234/352514135413134 is
> exactly that fraction, no more, no less. The same goes for `Sin[1]`, it is
> `Sin[1]` and no t0.8414709848078965 . There are functions that let one
> control when something gets evaluated to a numerical value -- otherwise it
> tries to keep everything in terms of integers and functions of integers and
> known entities like Pi/E. This is a common newbie trip-up, for example
> typing `Sin[1]+Cos[1]` will just spit out the same string, unless you pass
> it through `N[]`. But once you get used to it, it makes the software very
> powerful in you hands as you have greater knowledge of exactly where the
> floating point errors will occur, and a better degree of control over their
> effect. Of course, numerical methods like `NIntegrate[]`, `NDSolve[]` return
> floats by default
>
>
> This actually extends to expressions as well, you can give it (x-1)(x+1) and
> tell it to hold its form for a period of computation, and then release the
> hold later on whereupon it will simplify and substitute.
>
> Pretty neat.
>

Have a look at Sympy[1] or SAGE[2] for these things done the 'Python way'.

[1] http://sympy.org/en/index.html
[2] http://www.sagenb.org/

>> 3. [Continuum Computing - The Next Step in the Theory of

>> Computing](http://www.ee.iitb.ac.in/web/schedule/seminar/Talk_by_NarendraKarmarkar)
>>
>> --
>> --
>> The website for the club is http://wncc-iitb.org/
>> To post to this group, send email to wncc...@googlegroups.com
>> --- You received this message because you are subscribed to the Google
>> Groups "Web and Coding Club IIT Bombay" group.
>> To unsubscribe from this group and stop receiving emails from it, send an

>> email to wncc_iitb+...@googlegroups.com.

>> For more options, visit https://groups.google.com/groups/opt_out.
>
>

> --
> --
> The website for the club is http://wncc-iitb.org/
> To post to this group, send email to wncc...@googlegroups.com
> ---
> You received this message because you are subscribed to the Google Groups
> "Web and Coding Club IIT Bombay" group.
> To unsubscribe from this group and stop receiving emails from it, send an

> email to wncc_iitb+...@googlegroups.com.

[Dilawar Singh]

[Dilawar Singh]

[Saket Choudhary]

On 8 February 2014 14:02, Dilawar Singh <dilawar....@gmail.com> wrote:
>
> --
> --
> The website for the club is http://wncc-iitb.org/
> To post to this group, send email to wncc...@googlegroups.com
> ---
> You received this message because you are subscribed to the Google Groups
> "Web and Coding Club IIT Bombay" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to wncc_iitb+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>

> Python documentation warns its user,
>
> It's easy to forget that the stored value is an approximation to the
> original decimal fraction, because of the way that floats are displayed at
> the interpreter prompt. Python only prints a decimal approximation to the
> true decimal value of the binary approximation stored by the machine. If
> Python were to print the true decimal value of the binary approximation
> stored for 0.1, it would have to display
>
>>>> 0.1
> 0.1000000000000000055511151231257827021181583404541015625
>
> For instance, hand calculation and sympy calculation returned by following
> script which computes output at one resistor of potential divider differ a
> lot. Sympy returns (with or without using pi as symbol) 1.15789473684211 and
> hand calculation returns (using pi as symbol, cancellation occurs) 1.148....
>

I haven't looked at the calculations, but this looks interesting!

Saket

> Script, problem statement and a wrong solution to confuse students
>
> Two disaster caused by floating point errors
>
> PS: Most of my hand-calculation tends to be wrong. But anyway, I am counting
> on this one. -- Dilawar

[Dilawar Singh]

Another example (occurred today morning while simulation a model a rat's Neuron axon) is following. This equation  in model is computed by WollframAlfa (click me) but it could not be solved by GiNaC library.

[100%] Built target Axon
And the equation is (3.5014087480216975145E-4-(0.12095775674984045018)*r)*(0.015915494309189535082+(3.1830988618379066628)*r)^(-1)
terminate called after throwing an instance of 'GiNaC::pole_error'
  what():  power::eval(): division by zero
./run.sh: line 2: 24336 Aborted                 (core dumped) ./Axon

A `divide by zero` error occurs while substituting root returned by WolframAplha gives a near to zero result. A small perturbation here leads to very different values. This library has been used by octave to do symbolic computation and it is not as mature as Maple or Mathematica but it is wonderful example for newcomer to realize that CAS, despite what they claim, can go horribly wrong on certain kind of equations.

Here is the program. Change the variable `wolfram` at line 34 in axon.cpp file and run the run.sh file. You need to install cmake and libginac and clang++. You can use g++ also. Just change the compiler in CMakeLists.txt file.

--
Dilawar

On Saturday, February 8, 2014 7:38:21 PM UTC+5:30, Dilawar Singh wrote:

Neo: Why am I here?

Architect: Your life is the sum of a remainder of an unbalanced equation inherent to the programming of the matrix. you are the eventuallity of an anomaly which despite my sincerest efforts I have been unable to eliminate from what is otherwise a harmony of mathematical precision. While it remains a burden asciduously avoided it is not unexpected, and thus not beyond a measure of control. Which has led you inexorably….here

• From the movie 'Matrix' or something like that