# How to Break Math
## Numerical Computing with Floating-Point Numbers

Real numbers may have an infinite number of decimal digits: for example, $\frac{1}{3} = 0.333\dots$, and $\pi = 3.14159\dots$. However, computers do not have an infinite amount of storage, and for efficiency, most computer programs represent real numbers using floating-point values.

Floating point is similar to scientific notation. For example, the speed of light can be written simply as $299792458\,ms^{-1}$ or in scientific notation, $2.99792458 \times 10^8\,ms^{-1}$. A number in scientific notation has three parts: the *significand* or *mantissa*, the *base*, and the *exponent*:
\begin{equation}
299792458 = \underbrace{2.99792458}_{\text{significand}} \times {\underbrace{10}_{\text{base}}}^{\hspace{-0.6em}\overbrace{8}^{\text{exponent}}}
\end{equation}

By convention, the significand has one non-zero digit before the point (unless the number is exactly zero).

In binary floating point, the base is always 2, and the signficand is made up entirely of the digits 0 and 1. For example, the number 42 can be represented as a binary integer as 101010, or in binary floating point:
\begin{equation}
42 = 1.01010 \times 2^5
\end{equation}

Computers use the [IEEE-754 standard floating point format](https://en.wikipedia.org/wiki/IEEE_754), which defines standard "single" and "double" precision floating-point numbers. The bits which represent an IEEE-754 value are divided into three fields: the sign bit, the exponent, and the significand.




A sign bit of `0` represents a positive number, while a sign bit of `1` represents a negative number. 

The leading '1' in the significand is not stored, so the value `0b0100...` represents the signficand '1.0100...'. 

The exponent is *biased* - it is read as an integer value, from which the value of the bias is subtracted. For example, in single precision the bias is 127, so the exponent `0b10000100` (decimal 132) is read as (132-127) = 5, and the exponent `0b01111100` (decimal 124) is read as (124-127) = -3.


```julia
function print_bitstring(a::Float32)
    s = bitstring(a)
    printstyled(s[1], color=:blue)
    print(" ")
    printstyled(s[2:9], color=:green)
    print(" ")
    printstyled(s[10:32], color=:red)
    println()
end

function print_bitstring(a::Float64)
    # TODO complete - print the bits of the double-precision floating-point number a, separated into fields
end

function exponent_ieee754(a::Float32)
    s = bitstring(a)
    e = parse(Int32,s[2:9],base=2)
    return e - (2^7-1)
end

function exponent_ieee754(a::Float64)
    # TODO complete - extract exponent field from the double-precision floating point number a
end
    
x = Float32(42)
print(x," = ")
print_bitstring(x)
println("exponent = ", exponent_ieee754(x))

y = Float32(0.15625)
print(y," = ")
print_bitstring(y)
println("exponent = ", exponent_ieee754(y))

```

### Rounding Error

Rounding error occurs whenever an input value or the result of a calculation cannot be represented exactly in the chosen floating-point format, for example:


```julia
print_bitstring(0.1)
print_bitstring(Float64(pi))
print_bitstring(1.0/3.0)
```

### Cancellation

Cancellation occurs when two close-together values are subtracted, leading to a large relative error in the result. For example:


```julia
epsilon = 1e-8
x = (1+epsilon)
print("x ")
print_bitstring(x)
y = (1.0)
print("y ")
print_bitstring(y)
difference = epsilon
z = x-y
print("z ")
print_bitstring(z)
println("Correct result is ", difference, " but got ", z, " (relative error ", (difference-z)/difference, ")")
```

## Rump's Example

Siegfried Rump created a [nasty example](http://www.ti3.tu-harburg.de/paper/rump/Ru88a.pdf) of floating point error. He asks us to compute the value of the function:

$333.75 b^6 + a^2 (11 a^2 b^2 – b^6 – 121 b^4 – 2) + 5.5 b^8 + \frac{a}{(2b)}$

for $a = 77617.0$ and $b = 33096.0$ . We can define a Julia function to compute the value using floating-point numbers of type `T`:


```julia
function rump(T)
  a = T(77617.0)
  b = T(33096.0)
  a2 = a*a
  b2 = b*b
  b4 = b2*b2
  b6 = b4*b2
  b8 = b4*b4
  return T(333.75)*b6 + a2 * (11*a2*b2 - b6 - 121*b4 - 2) + T(5.5)*b8 + (a/(2*b))
end
```

Let's compute that function in single-precision floating point:


```julia
rump(Float32)
```

To check the correctness, let's compute it again in double-precision:


```julia
rump(Float64)
```

Oh dear! Our results aren't even remotely the same - they're out by more than eight orders of magnitude. 

Let's try increasing the precision still further. Julia supports arbitrary-precision arithmetic using [GNU MPFR](https://www.mpfr.org/) with a special type `BigFloat`. We can use this to create a function that will compute Rump's example using a chosen number of bits of precision in the significand:


```julia
function rump(numbits::Int)
  Base.MPFR.setprecision(numbits)
  return rump(BigFloat)
end
```

First, let's replicate our previous results using BigFloat. First single-precision, which in [IEEE-754](https://en.wikipedia.org/wiki/IEEE_754) has 24 bits in the significand:


```julia
rump(24)
```

Next, double precision, which has 53 bits in the significand:


```julia
rump(53)
```

Now let's try computing in quadruple precision, which in IEEE-754 has 113 bits in the significand:


```julia
rump(113)
```

This looks even worse - now even the sign of the result has changed!


```julia
for i = 113:130
    println("Precision ", i, " result is ", rump(i))
end
```

Rump gives the exact answer as $\frac{-54767}{66192} = -0.8273960599\dots$

Why do we need such high precision to correctly compute Rump's example? Use the `bitstring_ieee754(Float64)` function you created earlier to print the representation of some of the subexpressions that are computed in the function, to see if you can understand why the error is so great.


```julia
function rump_print(T)
  a = T(77617.0)
  b = T(33096.0)
  a2 = a*a
  b2 = b*b
  b4 = b2*b2
  b6 = b4*b2
  b8 = b4*b4
  # TODO print some of the subexpressions computed below
  intermediate_value1 = # ?
  intermediate_value2 = # ?
  print_bitstring(intermediate_value1)
  print_bitstring(intermediate_value2)
  return T(333.75)*b6 + a2 * (11*a2*b2 - b6 - 121*b4 - 2) + T(5.5)*b8 + (a/(2*b))
end

result = rump_print(Float64)
print_bitstring(result)
```

# Patriot vs. Scud

The [Patriot missile battery](http://www-users.math.umn.edu/~arnold//disasters/patriot.html) counted time using a clock which ticked ten times each second. To compute the expected location of an incoming target, the number of ticks was divided by 10 to produce a fixed-point value representing the number of seconds since the system started up. We can simulate this using a single-precision floating point value.

Modify the code below to compute the actual time in seconds given the current value of the clock. Why is `current_time_seconds` different? Hint: try printing out the bitstring representation of `tick_length`.


```julia
using Printf
const TICKS_PER_SECOND = 10
const SCUD_MISSILE_SPEED_METRES_PER_SECOND = 1676

function on_the_hour(ticks)
  return (ticks % (TICKS_PER_SECOND * 60 * 60) == 0)
end

Base.MPFR.setprecision(20)
tick_internal = BigFloat(1 / TICKS_PER_SECOND)
tick_length = Float32(tick_internal)

current_ticks = 0
current_time_seconds = 0.0
one_hundred_hours = 100 * 60 * 60 * TICKS_PER_SECOND
for i=0:one_hundred_hours
    current_ticks += 1
    current_time_seconds += tick_length
    if (on_the_hour(current_ticks))
        actual_time_seconds = # TODO compute correct time in seconds
        h = actual_time_seconds / 3600
        error_time_seconds = actual_time_seconds - current_time_seconds
        @printf("hour = %3d; error (time) = %.2e s; error (target position) %.2e m\n", h, 
            error_time_seconds, error_time_seconds*SCUD_MISSILE_SPEED_METRES_PER_SECOND)
    end
end
```

9c9a07f2ffd629c01005356be2ff12b6ff94e711
ipynb
Jupyter Notebook
numerics.ipynb
milthorpe/numerics
e8c93da721a205be584c32b15da45a2c3424a208
CC-BY-4.0
2020-03-27T08:01:31.000Z
2020-03-27T08:01:31.000Z
numerics.ipynb
milthorpe/numerics
e8c93da721a205be584c32b15da45a2c3424a208
CC-BY-4.0
numerics.ipynb
milthorpe/numerics
e8c93da721a205be584c32b15da45a2c3424a208
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES