You are viewing a plain text version of this content. The canonical link for it is here.
Posted to issues@commons.apache.org by "Alex Herbert (Jira)" <ji...@apache.org> on 2021/11/29 08:58:00 UTC
[jira] [Comment Edited] (NUMBERS-177) Scaled complementary error function
[ 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)