Better expm1

[tcaduser]

Better expm1

[j-fu]

So finally, thanks for this discussion, which lead to a more concise implementation.

[j-fu]

Thanks for the hint! Not sure why I was not aware of that...

But it is not that bad now either:

julia> x=1.0e-15
1.0e-15

julia> x/expm1(x)
0.9999999999999994

julia> VoronoiFVM.bernoulli_horner(x)
0.9999999999999994

So I need to see what I can gain. Will try it out.

[tcaduser]

I checked it in python, and the two formulas match to within 1e-16 at your threshold, which is at the limits of float64.

I guess the only difference would be the speed difference in the two implementations in julia. The implementation I use in my C++ solver is equivalent to:

y = expm1(x)
// this test may only be true for x == 0
if (x != y)
{
  b = x/y;
}
else
{
  // including x is probably overkill
  b = 1 / (1 + 0.5*x);
}

I am curious though about:

VoronoiFVM.jl/src/vfvm_functions.jl

Line 82 in 975d5b1

bm = x / (1.0 - 1.0 / expx)

and why not use:

bm = bp + x

as done here:

VoronoiFVM.jl/src/vfvm_functions.jl

Line 86 in 975d5b1

return y, x + y

[j-fu]

Well, I don't recall why...
I started a PR with an implentation using expm1 which indeed appears to be more robust.
However I indeed will keep the Horner scheme around zero. Bernoulli function with expm1 is good for small arguments, however the automatic derivative (which I use to create the Jacobians) fails:

julia> B(x)=x/expm1(x);

julia> B0(x)=1/(1+0.5x);

julia> B0(0.0)
1.0

julia> B0(nextfloat(0.0))
1.0

julia> ForwardDiff.derivative(B0,0)
-0.5

julia> ForwardDiff.derivative(B0,nextfloat(0.0))
-0.5

julia> ForwardDiff.derivative(B,nextfloat(0.0))
NaN

[j-fu]

As for the approximations which are now in the PR:

... thanks for nerd sniping me!

[tcaduser]

Very cool. I had to look up what Nerd Sniping was.

[tcaduser]

In my code, this is what I end up using for the derivative of B(x):

    const auto ex1 = expm1(x);

    if (x != ex1)
    {
      const auto ex2 = ex1 - (x * exp(x));
      ret = ex2;
      ret *= pow(ex1, -2);
    }
    else
    {
      DoubleType num = static_cast<DoubleType>(-0.5);
      DoubleType den = static_cast<DoubleType>( 1.0);
      num -= x / static_cast<DoubleType>(3.);
      den += x;
      ret = num / den;
    }

and so it works well even with extended precision. The full implementation is here:
https://github.com/devsim/devsim/blob/main/src/MathEval/Bernoulli.cc#L102

where I have to account for numerical issues I was having with the Android version of expm1 not working well with the x != expm1(x) test.