# Is there an implementation of Trefethen's contour integrals for evaluating \[1 - exp(-x)\] / x as x-\>0?

**URL:** <https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680>\
**Category:** General Usage\
**Tags:** package\
**Created:** [March 22, 2021, 2:47am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680 "2021-03-22T02:47:32Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [March 22, 2021, 2:47am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/1 "2021-03-22T02:47:33Z")

</div>

This kind of evaluation frequently occurs in numerical analysis, and I believe there should be an implementation in Julia. If so, I am happy to use it rather than re-invent the wheel:) Thanks!

---

<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 22, 2021, 4:48am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/2 "2021-03-22T04:48:57Z")

</div>

Do you mean `(1-exp(x))/x`? Otherwise, this diverges near 0.

---

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [March 22, 2021, 5:18am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/3 "2021-03-22T05:18:48Z")

</div>

Yes. Sorry for the ambiguity in the original title. I have revised it accordingly. Thanks!

Currently, I use TaylorSeries.jl to evaluate the expression when x \< 0.01. It works great!

```julia
using TaylorSeries

h(x) = (1 - exp(-x)) / x

t = Taylor1(6)
taylorh = (1 - exp(-t)) / t

x = 0.001
taylorh(x)

```

---

<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 22, 2021, 5:34am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/4 "2021-03-22T05:34:22Z")

</div>

If you want something that will give good (but not great) accuracy over the entire domain, `-expm1(x)/x` should do what you want. Combining that with `TaylorSeries` should get you `<1.5` ULP over the whole domain (I think).

---

<div class="post-metadata">

**Author:** ![ettersi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ettersi/32/6829_2.png) [@ettersi](https://discourse.julialang.org/u/ettersi)\
**Post date:** [March 22, 2021, 6:16am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/5 "2021-03-22T06:16:26Z")

</div>

Out of curiosity, how did you get the 1.5 ULP?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 22, 2021, 6:46am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/6 "2021-03-22T06:46:52Z")

</div>

`ApproxFun.jl` and friends might be able to do that.

---

<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 22, 2021, 1:33pm UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/7 "2021-03-22T13:33:03Z")

</div>

I’ve done a bunch of testing of `expm1` for [faster expm1 by oscardssmith · Pull Request #37440 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/37440) so I know that the current implimentation is around `.75` ULP. One extra rounding caused by the division adds a bit more than `.5` ULP since the error isn’t quite linear.

---

<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 22, 2021, 5:40pm UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/8 "2021-03-22T17:40:55Z")

</div>

NB:  
[ULP](https://en.wikipedia.org/wiki/Unit_in_the_last_place)

---

<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 22, 2021, 7:03pm UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/9 "2021-03-22T19:03:46Z")

</div>

> [@Oscar\_Smith](#):
>
> If you want something that will give good (but not great) accuracy over the entire domain, `-expm1(x)/x` should do what you want. Combining that with `TaylorSeries` should get you `<1.5` ULP over the whole domain (I think).

I’m not sure why you would need `TaylorSeries` at all here; `expm1` is already accurate for small arguments.

```julia
f(x) = iszero(x) ? float(one(x)) : -expm1(-x)/x

x = exp10.(range(-300,log10(700),length=10^7)); x = [-x; x]
err(x) = let exact=f(big(x)); Float64(abs(f(x)-exact) / abs(exact)); end
maximum(err, x) / eps(Float64)

```

gives 1.02 (≈1ulp).

> [@liuyxpp](#):
>
> Currently, I use TaylorSeries.jl to evaluate the expression when x \< 0.01. It works great!

That will give a much larger error — your `h(x) = (1 - exp(-x)) / x` has an error of ≈24ulps when `x=0.01`.

---

<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 22, 2021, 7:05pm UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/10 "2021-03-22T19:05:12Z")

</div>

Good catch! I was a little worried that the 0/0 in the limit would cause issues, but apparently not.

---

<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 22, 2021, 7:16pm UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/11 "2021-03-22T19:16:45Z")

</div>

> [@Oscar\_Smith](#):
>
> I was a little worried that the 0/0 in the limit would cause issues, but apparently not.

(You do need an explicit `iszero(x)` check as in my code above, as otherwise you will get NaN when `x` is exactly zero.)

---

<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 22, 2021, 10:51pm UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/14 "2021-03-22T22:51:00Z")

</div>

@stevengj, sorry for asking a dumb question: as nextfloat(0.) is equal to 5.0e-324, why not starting the range after -324? 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 23, 2021, 12:22am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/15 "2021-03-23T00:22:43Z")

</div>

> [@rafael.guerra](#):
>
> as nextfloat(0.) is equal to 5.0e-324, why not starting the range after -324? Thank you.

For this particular `f(x)` you actually get `1.0` as soon as `x` is less than about `1e-16`, so there is not much point in worrying about the exact bottom end of the range.

In general calculations, however, if you work with [subnormal values](https://en.wikipedia.org/wiki/Denormal_number) you may get larger relative errors because the precision is reduced. The smallest non-subnormal value is `floatmin(Float64) ≈ 2.2e-308`. I usually just remember `1e-300` as a “round” approximation for this, which is why my quick test code started there.

---

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [March 23, 2021, 1:44am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/16 "2021-03-23T01:44:42Z")

</div>

Thank you! I only expect an error below 1e-8. So I am satisfied for now:) Thanks for all the in-depth discussion here. Although I use it to evaluate the Deby function and the h function, i.e.

```julia
g(x, f) = (2/x^2)(e^{-xf} + xf - 1) = (2/x)[f - h(x, f)]

h(x, f) = (1 - e^{-xf}) / x

```

These kind of expression is also very useful for computing the coefficient for the ETDRK4 solver.

---

<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 23, 2021, 4:52am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/17 "2021-03-23T04:52:47Z")

</div>

> [@liuyxpp](#):
>
> `(2/x^2)(e^{-xf} + xf - 1)`

To compute this function accurately you need more than just `expm1`, because the numerator cancels both the 0th- and 1st-order terms in `xf` (for small `xf`). I’m not sure if one of the [tricks for `expm1`](https://math.stackexchange.com/questions/185215/what-is-a-nice-way-to-compute-fx-x-expx-1) can be adapted here or if you just have to use a Taylor series for small `xf`.

If you do use a Taylor series and care about performance, I would hard-code the coefficients and use `evalpoly`, e.g. this gives at least 14 digits:

```julia
function g(x,f) 
    xf = x*f
    abs(xf) > 0.1 && return 2(exp(-xf) + xf - 1) / x^2
    return f^2 * evalpoly(xf, (1, -1/3, 1/12, -1/60, 1/360, -1/2520, 1/20160, -1/181440, 1/1814400))
end

```

---

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [March 23, 2021, 5:07am UTC](https://discourse.julialang.org/t/is-there-an-implementation-of-trefethens-contour-integrals-for-evaluating-1-exp-x-x-as-x-0/57680/18 "2021-03-23T05:07:16Z")

</div>

Thanks! I will use your function. It is exactly what I want.