Skip to content

Simplify cases for erf #4754

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
phannebohm wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Simplify cases for erf #4754

Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
fbc94ef
Simplify cases for `erf`
phannebohm Feb 28, 2026
c45d0e4
Fix linear approximation
phannebohm Mar 2, 2026
301d858
Fix indentation
phannebohm Mar 2, 2026
b5a3033
Use spaces after walrus operator
phannebohm Mar 2, 2026
3f99557
Use actual machine epsilon
phannebohm Mar 2, 2026
1ae4265
Use rounding epsilon, not interval epsilon
phannebohm Mar 2, 2026
8d298d0
Squash literals
phannebohm Mar 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
89 changes: 43 additions & 46 deletions Modelica/Math/Special.mo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,53 +7,50 @@ package Special "Library of special mathematical functions"
input Real u "Input argument";
output Real y "= 2/sqrt(pi)*Integral_0_u exp(-t^2)*dt";
protected
Boolean inv;
Real z;
Real zz;
constant Real y1 = 1.044948577880859375;
constant Real P[5] = {0.0834305892146531832907,
-0.338165134459360935041,
-0.0509990735146777432841,
-0.00772758345802133288487,
-0.000322780120964605683831};
constant Real Q[5] = {1,
0.455004033050794024546,
0.0875222600142252549554,
0.00858571925074406212772,
0.000370900071787748000569};
Boolean inv;
Real z;
Real zz;
constant Real y1 = 1.044948577880859375;
constant Real P[5] = {0.0834305892146531832907,
-0.338165134459360935041,
-0.0509990735146777432841,
-0.00772758345802133288487,
-0.000322780120964605683831};
constant Real Q[5] = {1,
0.455004033050794024546,
0.0875222600142252549554,
0.00858571925074406212772,
0.000370900071787748000569};
algorithm
if u >= 0 then
z :=u;
inv :=false;
else
z :=-u;
inv :=true;
end if;

if z < 0.5 then
if z <= 0 then
y := 0;
elseif z < 1.0e-10 then
y := z*1.125 + z*0.003379167095512573896158903121545171688;
else
// Maximum Deviation Found: 1.561e-17
// Expected Error Term: 1.561e-17
// Maximum Relative Change in Control Points: 1.155e-04
// Max Error found at double precision = 2.961182e-17
zz := z*z;
y := z*(y1 + Internal.polyEval(P, zz)/Internal.polyEval(Q, zz));
end if;

elseif z < 5.8 then
y :=1 - Internal.erfcUtil(z);

else
y :=1;
end if;

if inv then
y :=-y;
end if;
if u >= 0 then
z := u;
inv := false;
else
z := -u;
inv := true;
end if;

if z < sqrt(1.5*Modelica.Constants.eps) then
// Maclaurin series to first order
// alternating series estimated relative error < eps
y := z*2.0/sqrt(Modelica.Constants.pi);
elseif z < 0.5 then
// Maximum Deviation Found: 1.561e-17
// Expected Error Term: 1.561e-17
// Maximum Relative Change in Control Points: 1.155e-04
// Max Error found at double precision = 2.961182e-17
zz := z*z;
y := z*(y1 + Internal.polyEval(P, zz)/Internal.polyEval(Q, zz));
elseif z < 5.8 then
y := 1 - Internal.erfcUtil(z);

else
y := 1;
end if;

if inv then
y := -y;
end if;

annotation (Documentation(info="<html>
<h4>Syntax</h4>
Expand Down