# Numerical precision

**URL:** <https://discourse.julialang.org/t/numerical-precision/78749>\
**Category:** General Usage\
**Tags:** question, precision\
**Created:** [March 30, 2022, 2:43pm UTC](https://discourse.julialang.org/t/numerical-precision/78749 "2022-03-30T14:43:38Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![F-YF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/f-yf/32/17363_2.png) [@F-YF](https://discourse.julialang.org/u/F-YF)\
**Post date:** [March 30, 2022, 2:43pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/1 "2022-03-30T14:43:38Z")

</div>

Recently I had a small problem with Julia that confused me. The code is shown below:

```julia

Julia> x=2.0e-10
2.0e-10

Julia> -1.0*(-0.5*x^3 - 3.0*x^2 - 7.5*x - 7.5)*exp(-x)/x^6 + 1.0*(0.5*x^3 - 3.0*x^2 + 7.5*x - 7.5)*exp(x)/x^6
0.0

Julia> x=2.1e-10
2.1e-10

Julia> -1.0*(-0.5*x^3 - 3.0*x^2 - 7.5*x - 7.5)*exp(-x)/x^6 + 1.0*(0.5*x^3 - 3.0*x^2 + 7.5*x - 7.5)*exp(x)/x^6
-1.1150372599265312e43

```

I’m sorry if it’s a very simple question, but I don’t understand.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 30, 2022, 2:44pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/2 "2022-03-30T14:44:35Z")

</div>

[https://0.30000000000000004.com/](https://0.30000000000000004.com/)

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 30, 2022, 2:44pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/3 "2022-03-30T14:44:54Z")

</div>

> [@PSA: floating-point arithmetic](https://discourse.julialang.org/t/psa-floating-point-arithmetic/8678):
>
> Sometimes people are surprised by the results of floating-point calculations such as julia\> 5/6 0.8333333333333334 # shouldn't the last digit be 3? julia\> 2.6 - 0.7 - 1.9 2.220446049250313e-16 # shouldn't the answer be 0? These are not bugs in Julia. They’re consequences of the IEEE-standard 64-bit binary representation of floating-point numbers that is burned into computer hardware, which Julia and many other languages use by default. Brief explanation You can t…

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 30, 2022, 2:52pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/4 "2022-03-30T14:52:19Z")

</div>

In particular, you are seeing a [catastrophic cancellation](https://en.wikipedia.org/wiki/Catastrophic_cancellation).

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 30, 2022, 5:42pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/5 "2022-03-30T17:42:16Z")

</div>

Results will behave by using big floats.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 30, 2022, 5:43pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/6 "2022-03-30T17:43:11Z")

</div>

Only kind of. The errors will be smaller, but still non-zero.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 30, 2022, 6:06pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/7 "2022-03-30T18:06:09Z")

</div>

> [@rafael.guerra](#):
>
> Results will behave by using big floats.

@F-YF posted an ill-conditioned sum whose terms (in exact arithmetic) should _exactly_ cancel to give zero. _Any_ roundoff error can give a nonzero answer, which is an _infinite_ relative error compared to `0.0` (_no_ digits are correct), no matter how much precision you have.

(From another perspective, compared to the magnitude of the _summand terms_ \sim 1/x^6 \approx 10^{58}, an error of 10^{43} is quite small: 10^{43} / 10^{58} = 10^{-15}.)

This is exactly the same as the issue discussed in another thread here: ["sum([1.0, 1e100, 1.0, -1e100])" fails - #10 by stevengj](https://discourse.julialang.org/t/sum-1-0-1e100-1-0-1e100-fails/76502/10)

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [March 30, 2022, 7:01pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/8 "2022-03-30T19:01:25Z")

</div>

The general guidance is to mathematically rewrite your expressions to avoid catastrophic cancellation. In this case Wolfram Alpha simplifies this to

(x(x^2 + 15)\cosh(x) - 3(2x^2 + 5)\sinh(x))/x^6

This doesn’t exhibit the same cancellation issue:

```julia
julia> f(x) = (x*(x^2 + 15)*cosh(x) - 3(2x^2 + 5)*sinh(x))/x^6
f (generic function with 1 method)

julia> f(2e-10)
0.0

julia> f(2.1e-10)
0.0

```

It could exhibit catastrophic cancellation elsewhere, however. But this is the general way you need to deal with such problems.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [March 30, 2022, 7:34pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/9 "2022-03-30T19:34:55Z")

</div>

I tried [herbie](https://herbie.uwplse.org/demo/780219c7723b885dd2f9c6cf2547b4876ca48160.ea9754c765d6f2e3443236047de835c7bf56f0de/graph.html) for fun, and it suggested:  
1.2025012025012024 \cdot 10^{-5} \cdot {x}^{5} + \left(0.0005291005291005291 \cdot {x}^{3} + x \cdot 0.009523809523809525\right)

Not exact, but the error is small

```julia
julia> function f(x)
         (1.2025012025012024e-5 * ^(x, 5.0)) + ((0.0005291005291005291 * ^(x, 3.0)) + (x * 0.009523809523809525));
       end
f (generic function with 1 method)

julia> f(2.1e-10)
2.0000000000000004e-12

julia> f(2.0e-10)
1.904761904761905e-12

```

It used a Taylor expansion to get this. However, their plot suggests it should be good for a wide range of inputs.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 30, 2022, 8:36pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/10 "2022-03-30T20:36:36Z")

</div>

> [@Oscar\_Smith](#):
>
> Only kind of. The errors will be smaller, but still non-zero.

Humble question: using big floats the error is close to zero (`~2e-12`) and a 5% perturbation of the input doesn’t catastrophically explode the output from near 0 to `~ -1e43`. This gain in stability should be important too, right?

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [March 30, 2022, 8:59pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/11 "2022-03-30T20:59:01Z")

</div>

right

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 30, 2022, 10:30pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/12 "2022-03-30T22:30:09Z")

</div>

> [@rafael.guerra](#):
>
> Humble question: using big floats the error is close to zero ( `~2e-12` ) and a 5% perturbation of the input doesn’t catastrophically explode the output from near 0 to `~ -1e43` .

Compared to what? Both of `2e-12` and `1e43` are infinitely bigger than `0.0`. Both are a catastrophic explosion of the forward error (infinite relative error, no correct significant digits).

> [@rafael.guerra](#):
>
> This gain in stability should be important too, right?

“Stability” has a technical meaning in numerical analysis, and whether an algorithm is “stable” or “unstable” is _independent of the precision_. It’s a category error to talk about “gaining stability” by increasing the precision.

This algorithm for f(x) is numerically unstable (regardless of the precision), I think (though I haven’t performed a formal analysis). You can only make the computation stable by changing the algorithm (e.g. by refactoring as @StefanKarpinski suggests).

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [March 30, 2022, 11:36pm UTC](https://discourse.julialang.org/t/numerical-precision/78749/13 "2022-03-30T23:36:50Z")

</div>

> [@stevengj](#):
>
> Compared to what? Both of `2e-12` and `1e43` are infinitely bigger than `0.0` .

Using big floats and plotting around `2e-12` the `cosh()...` equivalent expression provided above by @StefanKarpinski , we get:

 ![Plot_big_floats](https://global.discourse-cdn.com/julialang/original/3X/7/8/783de72e236dd90781a63f72b6cfad03f248d679.png)

So, the formula does not evaluate to zero at `x0 = 2.0e-12`.  
In this case, if we perturb the input by 5% from `x0` to `x0+dx = 2.1e-12`, we would not expect `dy` to vary by `1e48` (in absolute terms) when evaluating the function between x0 and x0+dx.

And please do not take this as a challenge, as I am light-years away from your knowledge. It is just my limited understanding. Thank you.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 31, 2022, 12:02am UTC](https://discourse.julialang.org/t/numerical-precision/78749/14 "2022-03-31T00:02:16Z")

</div>

> [@rafael.guerra](#):
>
> So, the formula does not evaluate to zero at `x0 = 2.0e-12` .

You’re quite right, the exact value is not zero, so the forward relative error is finite, not infinite; I mis-read the formula.

The basic numerical instability, however, is that no matter what precision you use, for sufficiently small input `x` you will get a catastrophic cancellation in the original formula.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 31, 2022, 12:09am UTC](https://discourse.julialang.org/t/numerical-precision/78749/15 "2022-03-31T00:09:53Z")

</div>

> [@StefanKarpinski](#):
>
> This doesn’t exhibit the same cancellation issue:
> 
> ```julia
> julia> f(x) = (x*(x^2 + 15)*cosh(x) - 3(2x^2 + 5)*sinh(x))/x^6
> 
> ```

Unfortunately, it’s still susceptible to spurious underflow:

```julia
julia> f(1e-100)
NaN

julia> setprecision(2^16)
65536

julia> Float64(f(big"1e-100"))
9.523809523809523e-103

```

and it fact it looks like there is still a catastrophic cancellation in the numerator for small `x`, where `x * cosh(x) ≈ sinh(x)`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 31, 2022, 12:38am UTC](https://discourse.julialang.org/t/numerical-precision/78749/16 "2022-03-31T00:38:00Z")

</div>

I don’t see a clever trick to implement this function accurately except by Taylor expanding around x=0 and switching between this series and the exact formula at some precision-dependent threshold:

```julia
function f(x::Union{Float64,ComplexF64})
    if abs(x) < 1.0 # empirical threshold for double precision
        return x * evalpoly(x^2, (1/105, 1/1890, 1/83160, 1/6486480, 1/778377600, 1/132324192000, 1/30169915776000))
    else
        return (x*(x^2 + 15)*cosh(x) - 3(2x^2 + 5)*sinh(x))/x^6
    end
end

```