[
https://issues.apache.org/jira/browse/NUMBERS-177?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17455249#comment-17455249
]
Alex Herbert commented on NUMBERS-177:
--------------------------------------
h2. Erfcx(z < 0.5)
The previous post did not show performance for erfcx(z < 0.5). When |z| < 0.5
the erfc function is computed using 1 - erf. The erfcx result is:
{noformat}
(1 - erf(z)) * exp(z*z)
exp(z*z) - erf(z) * exp(z*z)
{noformat}
At small values of z the exp(z*z) result approaches 1. We can rearrange the
formula to use expm1:
{noformat}
exp(z*z) - erf(z) * exp(z*z)
(exp(z*z) - 1) - erf(z) * (exp(z*z) - 1) - erf(z) + 1
expm1(z*z) - erf(z) * expm1(z*z) - erf(z) + 1
{noformat}
When z is negative (0 > z > -0.5) all terms are positive and can be summed in
order of magnitude. When z is positive (0 < z < 0.5) then the term -erf(z) is
negative. This can create cancellation. However the magnitude of each term is
larger than the round-off error from adding the previous term to the sum. At
the limit:
{noformat}
a + b + c + d
a = -erf(0.5) * expm1(0.5^2) = -0.147
b = expm1(0.5^2) = 0.284
c = -erf(0.5) = -0.520
d = = 1.0
-0.147 +
0.284 = 0.136 +
-0.520 = -0.384 +
1.0 = 0.616 = erfcx(0.5)
{noformat}
Thus cancellation is not significant. I computed the result of erfcx(z < 0.5)
for 1e9 random values between 2^-27 and 0.5 across 10000 bins. Note that below
2^-27 the value exp(z*z) is 1 and the result erfcx(z) = 1 - erf(z). The terms
were added in standard precision or added using 128-bit precision (xp). The
following is a plot of the max error and the difference between the two results:
!erfcx.small.png!
Summing in extended precision is better. But it cannot overcome the main source
of error which is error in the term erf(z). This is largest term in magnitude
that is computed (since the term 1 is exact). The difference in the max error
is very low. Thus summing the terms in extended precision is not a useful
addition.
Note that the step in the max error at large z occurs around z = 0.4736. This
corresponds to changes in the precision of expm1(z*z) and erf(z):
{noformat}
sqrt(log1p(0.25)) ~ 0.4723
erf inverse(0.5) ~ 0.4769
{noformat}
> Scaled complementary error function
> -----------------------------------
>
> Key: NUMBERS-177
> URL: https://issues.apache.org/jira/browse/NUMBERS-177
> Project: Commons Numbers
> Issue Type: New Feature
> Components: gamma
> Affects Versions: 1.1
> Reporter: Alex Herbert
> Priority: Minor
> Fix For: 1.1
>
> Attachments: erfc.png, erfcx.large.png, erfcx.medium.png,
> erfcx.small.png, expn.png
>
>
> Add a scaled complementary error function:
> {noformat}
> erfcx(x) = exp(x^2) * erfc(x)
> {noformat}
> For small x this can use the current rational function approximations for erf
> or erfc and remove the scaling down by exp(-x^2). For large x this can use
> the [continued fraction expansion of
> erfc|https://en.wikipedia.org/wiki/Error_function#Continued_fraction_expansion]
> and remove the scaling by e^-z:
> {noformat}
> z 1
> erfc z = -------- e^-z^2 -----------------------------
> sqrt(pi) z^2 + 1/2
> -----------------------
> 1 + 1
> -------------------
> z^2 + 3/2
> -------------
> 1 + 2
> ---------
> z^2 + ...
> {noformat}
> At very large x the rational function cannot be evaluated as {{z^2 + 0.5 ~
> z^2}} and the following approximation can be used:
> {noformat}
> 1 1
> erfcx z = -------- * -
> sqrt(pi) z
> {noformat}
> This occurs at 2^26 or approximately 6.71e7.
--
This message was sent by Atlassian Jira
(v8.20.1#820001)
