Simplify cases for erf

[phannebohm]

The local variable z is by construction non-negative and the linear case gives y = 0 for z = 0.

This makes the symbolic derivative of erf continuous around zero. (No longer relevant as soon as #4753 is merged.)

This simplifies the algorithm and clarifies the constant values.

[CLAassistant]

All committers have signed the CLA.

[AHaumer]

@phannebohm did you sign and re-check the CLA? Thanks!

[MartinOtter]

I asked Claude 4.6 Sonnet (free version; clause.ai), and Claude gave a very detailed answer and especially the Modelica package erf_with_derivative.mo (attached as erf_with_derivative.txt)

• Provides the analytic derivative of erf(u)
• Implemented it as a Modelica function with derivative annotation
• Provided a test model: der(y) = erf(time)

I simulated in OpenModelica and got the following plot which looks good.

der(y) = er(sin(time)) from OpenModelica/OpenModelica#15118 (comment) works also without problems.

I asked Claude to provide a better implementation of Modelica.Math.Special.erf and after a while, Claude provided file Erfimpl.mo (attached as ErfImpl.txt). Claude states:

• The function value is computed via Padé rational approximations derived
  with 80-digit arithmetic (mpmath) from the Taylor series of erf(u)/u and
  erfc(x)*exp(x^2), and verified to achieve a maximum relative error below
  7e-16 (within 4 ULP of IEEE 754 double precision) over the whole real line.

I quickly checked with OpenModelica and this looks good (but more operations as Modelica.Math.Special.erf, not clear whether this is needed).

I do not have time to provide a proper pull request, with documentation and test model included in ModelicaTest. Hopefully, someone else can do it (the quickest solution is to just add erf_derivative(..) from file erf_derivative.txt to MSL, include the derivative annotation to erf(..) and add ErfDemo to ModelicaTest and adapt the documentations).
In a similar way, also the derivative of erfc(..) could be derived from Claude and added to MSL.

[AHaumer]

@MartinOtter I think your comment belongs rather to #4753 - here I provided a full PR (without Claude).
You (resp. Claude) didn't change the implementation of erf itself, but @phannebohm aims at an enhancement of erf.
The derivative of erf is trivial.

[HansOlsson]

@MartinOtter I think your comment belongs rather to #4753 - here I provided a full PR (without Claude). You (resp. Claude) didn't change the implementation of erf itself, but @phannebohm aims at an enhancement of erf. The derivative of erf is trivial.

I agree, but to me the function could be simplified further.

The proposed improvement has:

   z < 0.5 then
      if z < 1.0e-10 then
         y := z*1.125 + z*0.003379167095512573896158903121545171688;
      else

I can understand that 1.125 and 0.003379167095512573896158903121545171688 (from the original) are computed differently - but adding them to generate

   z < 0.5 then
      if z < 1.0e-10 then
         y := z*(1.125+0.003379167095512573896158903121545171688);
      else

would be straightforward, and avoids having to check if there was any error in the code.

[phannebohm]

I fixed the constant in the linear approximation, as suggested by @HansOlsson, and I computed a better error bound for that case as well, the relative error in the Maclaurin series is actually below machine precision until z = 2.5e-8 way beyond 1e-10. If we added the cubic term we could get that to $\sqrt[4]{10 \epsilon_M} \approx$ 2.1e-4.

I also cleaned up the formatting and simplified the nested cases for z.

[HansOlsson]

Oh, so z*1.125 + z*0.003379167095512573896158903121545171688; was just z*2/sqrt(pi) split into two parts, really good to have that changed.

[casella]

@phannebohm I'm not sure I get the point of this PR. The original motivation was to avoid getting a bad symbolic derivative close to zero.

As I understand, erf is the integral of a gaussian function. As such, it is not computable in closed form, so we resort to smart approximation functions. But it doesn't really make sense to rely on symbolic differentiation in this case, since the derivative is known analytically and is very simple.

Bottom line, I think what is really needed to get a good derivative is to add a derivative annotation.

Do you agree?

[phannebohm]

Bottom line, I think what is really needed to get a good derivative is to add a derivative annotation.

Do you agree?

Absolutely! As you say the better solution for that problem is #4753.

So the remaining purpose for this PR is to simplify, and potentially clarify the smart approximation for erf.
I figured, since I already looked at the function anyway, I might as well do the PR...

[HansOlsson]

Looks good

[HansOlsson]

Bottom line, I think what is really needed to get a good derivative is to add a derivative annotation.

Yes, and I assume we all agree on this, but just wanted it stated:

• functions like erf may have a derivative annotation.
• analytic functions in closed form like the derivative of erf should just have smoothOrder-annotation, allowing tools to safely differentiate them.

Trying to differentiate a function like erf without smoothOrder annotation is error-prone; as automatic differentiation for the original variant would have given incorrect results at 0, and here you have a nice closed form for the derivative - which can just be used.

[AHaumer]

@phannebohm pls. wait until #4753 is merged and pull those changes to your PR:
annotations(derivative) in erf and erfc and new functions erfDer and erfcDer

[phannebohm]

@phannebohm pls. wait until #4753 is merged and pull those changes to your PR: annotations(derivative) in erf and erfc and new functions erfDer and erfcDer

I think github is able to tell if the changes overlap and in that case it does not allow merging before the conflicts are resolved.