# Difference between f(1.0e-7) & f(1.0e-8)

**URL:** <https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [December 29, 2019, 10:09am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779 "2019-12-29T10:09:15Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Stephen](https://avatars.discourse-cdn.com/v4/letter/s/96bed5/32.png) [@Stephen](https://discourse.julialang.org/u/Stephen)\
**Post date:** [December 29, 2019, 10:09am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/1 "2019-12-29T10:09:15Z")

</div>

the limit of the following function at 0 is 1,

![QQ截图20191229180612](https://global.discourse-cdn.com/julialang/original/3X/d/2/d2d3dda3247943edcaed4c167d4dab9145c09ce3.png)

yet f(1.0e-8) gives 0 while f(1.0e-7) gives the correct answer, why?

```julia
julia> f(x)=2*(1-x*exp(-x)-exp(-x))/(1-exp(-x))^2
f (generic function with 1 method)

julia> f(1.0e-7)
0.9992008230547272

julia> f(1.0e-8)
0.0

```

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [December 29, 2019, 11:36am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/2 "2019-12-29T11:36:29Z")

</div>

The problem is [catastrophic cancellation](https://en.wikipedia.org/wiki/Loss_of_significance) due to limited accuracy of double precision floating point arithmetic:

```julia
julia> 1 - (1e-7) * exp(-1e-7)
0.99999990000001

julia> exp(-1e-7)
0.999999900000005

julia> 1 - (1e-7) * exp(-1e-7) - exp(-1e-7)
4.9960036108132044e-15

julia> 1 - (1e-8) * exp(-1e-8)
0.9999999900000001

julia> exp(-1e-8)
0.9999999900000001

julia> 1 - (1e-8) * exp(-1e-8) - exp(-1e-8)
0.0

```

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [December 29, 2019, 11:45am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/3 "2019-12-29T11:45:46Z")

</div>

Note also that your `f` is probably quite wrong also for `x = 1e-7`:

```julia
julia> f(big"1e-7")
1.000000033333333333333322222222222222226190476190476189153439153575262308010033

```

Using textbook formulas involving very small numbers is often wrong in numerical computing.

---

<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:** [December 29, 2019, 12:42pm UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/4 "2019-12-29T12:42:49Z")

</div>

The following implementation computes your g(\delta) accurately in double precision, by switching over to a Taylor series for small \delta and exploiting the `expm1` function for larger \delta.

```julia
function g(δ::Union{Float64,ComplexF64})
    if abs(δ) < 0.2
        return 1 + δ*@evalpoly(δ^2, 1/3,-1/90,1/2520,-1/75600,1/2395008,-691/54486432000)
    else
        y = expm1(-δ)
        return -2(y + δ*exp(-δ)) / y^2
    end
end

```

(It’s possible that there is some clever algebraic transformation that lets you avoid the explicit Taylor series, but I’m not seeing it at the moment.)

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [December 29, 2019, 2:39pm UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/5 "2019-12-29T14:39:34Z")

</div>

you can also take the fact of the existence of `expm1` and expanding each term:

```julia
f1(x)= 2*(1-x*exp(-x)-exp(-x))/(1-exp(-x))^2 #for comparison
function f2(x)
    denom= -expm1(-x)
    return 2/denom+ 2*x/denom- 2*x/(denom*denom)
end

```

Testing for equality on other number

```julia
julia> f1(4)
1.885271061051406

julia> f2(4)
1.8852710610514052

```

near zero:

```julia
julia> f1(1e-7)
0.9992008230547272

julia> f2(1e-7)
1.0000000335276127

julia> f2(1e-8)
1.0

julia> f1(1e-8)
0.0

```

to be fair it also fails, but around `1e-11`

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [December 29, 2019, 3:17pm UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/6 "2019-12-29T15:17:35Z")

</div>

Often times for numerical computing, in the field I was in, we would either rigorously determine, or crudely determine limits for significance/guard digits. Rules of thumb for floats were usually 1e-5, for doubles usually 1e-10 (because we used ffts and things a lot).

Main take away is, there’s error in every calculation you do basically. There’s also a finite resolution to what a digital object can store/handle. When you do lots of calculations that error can propagate, or do unexpected things.

---

<div class="post-metadata">

**Author:** ![StevenSiew](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevensiew/32/218393_2.png) [@StevenSiew](https://discourse.julialang.org/u/StevenSiew)\
**Post date:** [December 29, 2019, 11:20pm UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/7 "2019-12-29T23:20:48Z")

</div>

> [@anon92994695](#):
>
> Often times for numerical computing, in the field I was in, we would either rigorously determine, or crudely determine limits for significance/guard digits. Rules of thumb for floats were usually 1e-5, for doubles usually 1e-10 (because we used ffts and things a lot).

Does this means that in any serious calculations where human lives are at stake, you need to do the calculations twice? Once in Float64 and once in Float32. Then compare the results to make sure they come up with roughly the same value in case there are any numerical instabilities?

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [December 29, 2019, 11:25pm UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/8 "2019-12-29T23:25:47Z")

</div>

Well I’m not so sure about the answer to that. Another thing that can be done, instead, would be to plot a range of solutions near the answer, and test results of known similar conditions.

I mostly did work like this. When something goes awry you’ll be able to tell most of the time. Complex real and imag values become underdetermined you’ll see steppy swaps. Really though you can either measure or find the floating point error for most calculations. Julia is actually amazing for this… I need to write a paper about it - but have negative time to do so. I have tons of results and code for it though…

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [December 29, 2019, 11:39pm UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/9 "2019-12-29T23:39:33Z")

</div>

When someone has done a careful analysis, this is not necessary—in some cases one can prove that the result has a certain accuracy.

---

<div class="post-metadata">

**Author:** ![Stephen](https://avatars.discourse-cdn.com/v4/letter/s/96bed5/32.png) [@Stephen](https://discourse.julialang.org/u/Stephen)\
**Post date:** [December 29, 2019, 11:54pm UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/10 "2019-12-29T23:54:14Z")

</div>

thanks all of you for such detailed explanation and diverse solutions

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [December 30, 2019, 12:43am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/11 "2019-12-30T00:43:08Z")

</div>

> [@StevenSiew](#):
>
> Does this means that in any serious calculations where human lives are at stake, you need to do the calculations twice? Once in Float64 and once in Float32. Then compare the results to make sure they come up with roughly the same value in case there are any numerical instabilities?

Although it is not entirely foolproof, this method is likely to work well in practice. For a rather extensive discussion about such techniques, see for example this paper by Kahan:

> **[Mindless.pdf](https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf)**
>
> 375.12 KB

The discussion about possible round-off error analysis schemes and Kahan’s opinion about them starts at page 13.  
The first of the five schemes listed in this paper is the one you mention: running the same computation several times with different FP precisions. The second scheme (Kahan’s favorite) is very similar: it consists in running the same computation several times with the same FP precision, but with different rounding modes.

---

<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:** [December 30, 2019, 1:00am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/12 "2019-12-30T01:00:22Z")

</div>

Sometimes numerical work may inform decisions and influence actions relevant to another’s health, wealth, or contentment. I have found it better to compute more accurately than may obtain using `Float32` or `Float64`, then converting the likely more accurate result to `Float32` or `Float64`. The choice of computational type and the choice of presentational type are made with an understanding of the specifics of the work and its context.

There are some publicly available Julia packages designed for this work. For your purposes, I suggest becoming familiar with [DoubleFloats.jl](https://github.com/JuliaMath/DoubleFloats.jl). Where algorithmic choices, data uncertainties, or natural inexactness suggest characterizing the results as spans, The types `ArbReal` and `ArbComplex` from [ArbNumerics.jl](https://github.com/JeffreySarnoff/ArbNumerics.jl) do just that. Also see packages from [JuliaIntervals](https://github.com/JuliaIntervals).

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [December 30, 2019, 2:09am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/13 "2019-12-30T02:09:42Z")

</div>

I saw your talk on this!!!

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [December 30, 2019, 7:52am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/14 "2019-12-30T07:52:33Z")

</div>

You can also use `TaylorSeries.jl` to do Taylor series manipulations:

```julia
julia> using TaylorSeries

julia> n = 16
16

julia> t = Taylor1(Rational{Int}, n)
 1//1 t + 𝒪(t¹⁷)

julia> ee = sum( (-t)^i / factorial(i) for i in 0:n )
 1//1 - 1//1 t + 1//2 t² - 1//6 t³ + 1//24 t⁴ - 1//120 t⁵ + 1//720 t⁶ - 1//5040 t⁷ + 1//40320 t⁸ - 1//362880 t⁹ + 1//3628800 t¹⁰ - 1//39916800 t¹¹ + 1//479001600 t¹² - 1//6227020800 t¹³ + 1//87178291200 t¹⁴ - 1//1307674368000 t¹⁵ + 1//20922789888000 t¹⁶ + 𝒪(t¹⁷)

julia> 2 * (1 - t * ee - ee) / ((1 - ee)^2)
 1//1 + 1//3 t - 1//90 t³ + 1//2520 t⁵ - 1//75600 t⁷ + 1//2395008 t⁹ - 691//54486432000 t¹¹ + 1//2668723200 t¹³ - 59513//156920924160000 t¹⁵ - 42397//94152554496000 t¹⁶ + 𝒪(t¹⁷)

```

(Note that although `exp(t)` is defined for a Taylor series, it currently converts the arguments to floats.)

You can also bound the remainder using `TaylorModels.jl`.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [December 30, 2019, 8:03am UTC](https://discourse.julialang.org/t/difference-between-f-1-0e-7-f-1-0e-8/32779/15 "2019-12-30T08:03:47Z")

</div>

That’s a neat writeup of very relevant problems and methods from a giant in the field, but it is surprising how jaded and bitter the tone of the article is.