Algoritmo de boost :: math :: erf

Question

¿hay detalles disponibles para el algoritmo detrás de la función erf de boost? La documentación del módulo no es muy precisa. Todo lo que descubrí es que varios métodos son mixtos. Para mí, parece variaciones de Abramowitz y Stegun.

• ¿Qué métodos se mezclan?
• ¿Cómo se mezclan los métodos?
• ¿Cuál es la complejidad de la función erf (tiempo constante)?

Sebastian

Answer 1

la documentación para Boost Math Toolkit tiene una larga lista de references, entre los que Abramowitz y Stegun. La interfaz erf-function contiene un parámetro de plantilla policy que se puede usar para controlar la precisión numérica (y, por lo tanto, su complejidad en tiempo de ejecución).

#include <boost/math/special_functions/erf.hpp> 
namespace boost{ namespace math{ 

template <class T> 
calculated-result-type erf(T z); 

template <class T, class Policy> 
calculated-result-type erf(T z, const Policy&); 

template <class T> 
calculated-result-type erfc(T z); 

template <class T, class Policy> 
calculated-result-type erfc(T z, const Policy&); 

}} // namespaces 

ACTUALIZACIÓN:

A continuación una copia literal de la sección "Aplicación" de la anterior referencia proporcionada a la ERF-función:

Implementación

Todas las versiones de estas funciones usan primero las fórmulas de reflexión habituales para hacer que sus argumentos sean positivos:

erf(-z) = 1 - erf(z); 

erfc(-z) = 2 - erfc(z); // preferred when -z < -0.5 

erfc(-z) = 1 + erf(z); // preferred when -0.5 <= -z < 0 

Las versiones genéricas de estas funciones se implementan en términos de la función gamma incompleta.

Cuando se reconoce el tamaño significante (mantisa) (actualmente para 53, 64 y 113 bits reales, más precisión única de 24 bits manejada mediante promoción a doble), se utiliza una serie de aproximaciones racionales diseñadas por JM.

Para z < = 0,5 entonces se utiliza una aproximación racional a erf, basado en la observación de que erf es una función impar y por lo tanto erf se calcula usando:

erf(z) = z * (C + R(z*z)); 

donde la aproximación racional R (z * z) está optimizado para el error absoluto: siempre que su error absoluto es lo suficientemente pequeño en comparación con la constante C, entonces cualquier error de redondeo incurrido durante el cálculo de R (z * z) desaparecerá efectivamente del resultado. Como resultado, el error para erf y erfc en esta región es muy bajo: el último bit es incorrecto en muy pocos casos.

Para z> 0.5 se observa que más de un pequeño intervalo [a, b) a continuación:

erfc(z) * exp(z*z) * z ~ c 

para alguna constante c.

Por lo tanto para z> 0,5 calculamos erfc usando:

erfc(z) = exp(-z*z) * (C + R(z - B))/z; 

De nuevo R (z - B) está optimizado para el error absoluto, y la constante C es la media de erfc (z) * exp (z * z) * z tomado en los puntos finales del rango.Una vez más, siempre que el error absoluto en R (z - B) sea pequeño en comparación con c, entonces c + R (z - B) se redondeará correctamente, y el error en el resultado dependerá únicamente de la precisión de la exp función. En la práctica, en todo menos en un número muy pequeño de casos, el error se limita al último bit del resultado. La constante B se elige de manera que el extremo izquierdo de la gama de la aproximación racional es 0.

Para grandes z en un intervalo [a, + ∞] a la aproximación anterior se modifica a:

erfc(z) = exp(-z*z) * (C + R(1/z))/z; 

Las aproximaciones racionales se explican en excruciating detail. Si necesita más detalles, siempre puede mirar el source code.