---
engine: julia
---

```{julia}
#| error: false
#| echo: false
#| output: false
using InteractiveUtils
```

```{julia}
#| error: false
#| echo: false
#| output: false
flush(stdout)
```

# Example of Floating-Point Arithmetic Stability



## Calculation of $\pi$ according to Archimedes

A lower bound for $2\pi$, the circumference of the unit circle, is obtained by the
sum of the side lengths of a regular $n$-gon inscribed in the unit circle.
The figure on the left shows how one can iteratively double the number of vertices starting from a square with side length $s_4=\sqrt{2}$.

:::{.narrow}

| Fig. 1 | Fig. 2 |
| :-:    | :-:   |
| ![](../images/pi1.png) | ![](../images/pi2.png)  |
: {tbl-colwidths="[57,43]"}

:::


The second figure shows the geometry of the vertex doubling.

With
$|\overline{AC}|= s_{n},\quad |\overline{AB}|= |\overline{BC}|= s_{2n},\quad |\overline{MN}| =a, \quad |\overline{NB}| =1-a,$ the Pythagorean theorem applied to triangles $MNA$ and
$NBA$ respectively yields
$$
  a^2 + \left(\frac{s_{n}}{2}\right)^2 = 1\qquad\text{resp.} \qquad
  (1-a)^2 + \left(\frac{s_{n}}{2}\right)^2 = s_{2n}^2
$$
Elimination of $a$ gives the recursion
$$
   s_{2n} = \sqrt{2-\sqrt{4-s_n^2}} \qquad\text{with initial value}\qquad
 s_4=\sqrt{2}
$$
for the length $s_n$ **of one** side of the inscribed regular
$n$-gon.


The sequence $(n\cdot s_n)$
converges monotonically from below to the
limit $2\pi$:
$$
    n\, s_n \rightarrow 2\pi % \qquad\text{and} \qquad {n c_n}\downarrow 2\pi
$$
The relative error of the approximation of 2π by an $n$-gon is:
$$
     \epsilon_n = \left| \frac{n\, s_n-2 \pi}{2\pi} \right|
$$


## Two Iteration Formulas^[by: Christoph Überhuber, "Computer-Numerik" Vol. 1, Springer 1995, Chap. 2.3]
The equation
$$
   s_{2n} = \sqrt{2-\sqrt{4-s_n^2}}\qquad \qquad \text{(Iteration A)}
$$
is mathematically equivalent to
$$
    s_{2n} = \frac{s_n}{\sqrt{2+\sqrt{4-s_n^2}}} \qquad \qquad \text{(Iteration B)}
$$

(Please verify!)

However, Iteration A is ill-conditioned and numerically unstable, as the following code demonstrates. The output shows the respective approximation for π.

```{julia}
using Printf

iterationA(s) =  sqrt(2 - sqrt(4 - s^2))
iterationB(s) = s / sqrt(2 + sqrt(4 - s^2))

s_B = s_A = sqrt(2)    # initial value 

ϵ(x) = abs(x - 2π)/2π  # rel. error   

ϵ_A = Float64[]   # vectors for the plot
ϵ_B = Float64[]
is  = Float64[]

@printf """                        approx. value of π
          n-gon       Iteration A            Iteration B
"""

for i in 3:35
    push!(is, i)
    s_A = iterationA(s_A)         
    s_B = iterationB(s_B)         
    doublePi_A = 2^i * s_A
    doublePi_B = 2^i * s_B
    push!(ϵ_A, ϵ(doublePi_A))
    push!(ϵ_B, ϵ(doublePi_B))
    @printf "%14i %20.15f  %20.15f\n" 2^i doublePi_A/2 doublePi_B/2 
end
```

While Iteration B stabilizes at a value correct within machine precision, Iteration A quickly becomes unstable. A plot of the relative errors $\epsilon_i$ confirms this:


```{julia}
using PlotlyJS

layout = Layout(xaxis_title="Iteration steps", yaxis_title="rel. error",
            yaxis_type="log", yaxis_exponentformat="power", 
            xaxis_tick0=2, xaxis_dtick=2)

plot([scatter(x=is, y=ϵ_A, mode="markers+lines", name="Iteration A", yscale=:log10),   
      scatter(x=is, y=ϵ_B, mode="markers+lines", name="Iteration B", yscale=:log10)],
      layout) 
 ```
 
## Stability and Cancellation

At $i=26$, Iteration B reaches a relative error on the order of machine epsilon:


```{julia}
ϵ_B[22:28]
```

Further iterations do not improve the result. It stabilizes at a relative error of approximately 2.5 machine epsilons:


```{julia}
ϵ_B[end]/eps(Float64)
```


Iteration A is unstable. Already at $i=16$ the relative error begins to grow again.

The cause is a typical cancellation effect. The side lengths $s_n$ become very small very quickly. Thus
$a_n=\sqrt{4-s_n^2}$ is only slightly smaller than 2, and when computing $s_{2n}=\sqrt{2-a_n}$ a typical cancellation effect occurs.


```{julia}
setprecision(80) # precision for BigFloat

s = sqrt(BigFloat(2))

@printf "     a = √(4-s^2) as BigFloat       and as Float64\n\n"

for i = 3:44
    s = iterationA(s)
    x = sqrt(4-s^2)
    if i > 20
        @printf "%i %30.26f %20.16f \n" i  x  Float64(x)
    end 
end


```

One sees the decrease in the number of significant digits. One also sees that using `BigFloat` with a mantissa length of 80 bits here only slightly delays the onset of the cancellation effect.

**Countermeasures against unstable algorithms usually require better, stable algorithms, not more precise machine numbers!**

:::{.content-hidden unless-format="xxx"}

Clearly, the computation of $2-a_n$ already shows a decrease in the number of significant digits (cancellation) relatively early,
before the computation of $a_n=\sqrt{4-s_n^2}$ itself leads to an unusable result due to cancellation.

:::