special_functions.cpp File Reference

#include <cmath>
#include "./special_functions.h"
#include "../_c_arithmetics/c_arithmetics.h"

Include dependency graph for special_functions.cpp:

Go to the source code of this file.

Functions
double	_sign (const double x)
	sign function.
double	erf_inv (const double x)
	Inverse error function.

Function Documentation

_sign()

double _sign

(

const double

x

)

sign function.

Definition at line 35 of file special_functions.cpp.

{
    return (x > 0) - (x < 0);
}

Referenced by erf_inv().

Here is the caller graph for this function:

erf_inv()

double erf_inv

(

const double

x

)

Inverse error function.

The function inverse is found based on Newton method using the evaluation of the error function erf from standard math library and its derivative. The Newton method here uses two refinements.

For further details on the algorithm, refer to: http://www.mimirgames.com/articles/programming/approximations- of-the-inverse-error-function/

The accuracy of this method for the whole input interval of [-1, 1] is in the order of 1e-15 compared to scipy.special.erfinv function.

Parameters

[in]

x

Input value, a float number between -1 to 1.

Returns

The inverse error function ranging from -INFINITY to INFINITY.

Definition at line 65 of file special_functions.cpp.

{
    // Output
    double r;
 
    // Check extreme values
    if (c_arithmetics::is_equal(x, 1.0) || c_arithmetics::is_equal(x, -1.0))
    {
        r = _sign(x) * INFINITY;
        return r;
    }
 
    double a[4] = {0.886226899, -1.645349621, 0.914624893, -0.140543331};
    double b[5] = {1.0, -2.118377725, 1.442710462, -0.329097515, 0.012229801};
    double c[4] = {-1.970840454, -1.62490649, 3.429567803, 1.641345311};
    double d[3] = {1.0, 3.543889200, 1.637067800};
 
    double z = _sign(x) * x;
 
    if (z <= 0.7)
    {
        double x2 = z * z;
        r = z * (((a[3] * x2 + a[2]) * x2 + a[1]) * x2 + a[0]);
        r /= (((b[4] * x2 + b[3]) * x2 + b[2]) * x2 + b[1]) * x2 + b[0];
    }
    else
    {
        double y = sqrt(-log((1.0 - z) / 2.0));
        r = (((c[3] * y + c[2]) * y + c[1]) * y + c[0]);
        r /= ((d[2] * y + d[1]) * y + d[0]);
    }
 
    r = r * _sign(x);
    z = z * _sign(x);
 
    // These two lines below are identical and repeated for double refinement.
    // Comment one line below for a single refinement of the Newton method.
    r -= (erf(r) - z) / ((2.0 / sqrt(M_PI)) * exp(-r * r));
    r -= (erf(r) - z) / ((2.0 / sqrt(M_PI)) * exp(-r * r));
 
    return r;
}

References _sign(), and c_arithmetics::is_equal().

Referenced by ConvergenceTools< DataType >::average_estimates(), and ConvergenceTools< DataType >::check_convergence().

Here is the call graph for this function:

Here is the caller graph for this function: