[
https://issues.apache.org/jira/browse/NUMBERS-177?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17450271#comment-17450271
]
Alex Herbert edited comment on NUMBERS-177 at 11/29/21, 8:57 AM:
-----------------------------------------------------------------
I have created a working implementation of this and evaluated it compared to
the Matlab implementation. Results are within 1 or 2 ULP for the the valid
range from -26.6 to 7e7. Note for negative x:
{noformat}
erfcx(x) = 2*exp(x*x) - erfcx(|x|)
{noformat}
The switch to the use of the continued fraction requires some further
investigation. The current implementation for erfc( x ) uses a rational
function approximation for 4.5 < x < 28. This is scaled down by exp(-x*x). At x
~ 26.54 the result is subnormal. At ~27.2 the result is 0:
{noformat}
erfc(x) = exp(-x*x) * P / Q
erfcx(x) = P / Q
{noformat}
Given that P / Q has been optimised to work when x < 28 and in the knowledge
that x > 26.54 will be sub-normal (and thus requires less precision) the value
of P / Q should be checked that it is better than using the rational function
approximation.
When x > 50 then the continued fraction converges very fast. It is possible to
evaluate the nth convergent directly without using an iterative method, e.g.
using the first 5 terms of the continued fraction. The region 26.5 < x < 50 may
require iteration.
was (Author: alexherbert):
I have created a working implementation of this and evaluated it compared to
the Matlab implementation. Results are within 1 or 2 ULP for the the valid
range from -26.6 to 7e7. Note for negative x:
{noformat}
erfcx(x) = 2*exp(x*x) - erfcx(|x|)
{noformat}
The switch to the use of the continued fraction requires some further
investigation. The current implementation for erfc(x) uses a rational function
approximation for 4.5 < x < 28. This is scaled down by exp(-x*x). At x ~ 26.54
the result is subnormal. At ~27.2 the result is 0:
{noformat}
erfc(x) = exp(-x*x) * P / Q
erfcx(x) = P / Q
{noformat}
Given that P / Q has been optimised to work when x < 28 and in the knowledge
that x > 26.54 will be sub-normal (and thus requires less precision) the value
of P / Q should be checked that it is better than using the rational function
approximation.
When x > 50 then the continued fraction converges very fast. It is possible to
evaluate the nth convergent directly without using an iterative method, e.g.
using the first 5 terms of the continued fraction. The region 26.5 < x < 50 may
require iteration.
> 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
>
>
> 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)
