tabla - Implementación vectorizable de la función de error complementario erfcf()

integral gaussiana ejemplos

Question

La función de error complementario, erfc , es una función especial estrechamente relacionada con la distribución normal estándar. Se utiliza con frecuencia en las estadísticas y las ciencias naturales (por ejemplo, problemas de difusión) donde las "colas" de esta distribución deben considerarse, y el uso de la función de error, erf , no es adecuado.

La función de error complementario estuvo disponible en la biblioteca matemática estándar ISO C99 como las funciones erfcf , erfc y erfcl ; estos fueron posteriormente adoptados en ISO C ++ también. Por lo tanto, el código fuente puede encontrarse fácilmente en implementaciones de código abierto de esa biblioteca, por ejemplo, en glibc .

Sin embargo, muchas implementaciones existentes son de naturaleza escalar, mientras que el hardware del procesador moderno está orientado a SIMD (explícitamente, como en la CPU x86, o implícitamente, como en las GPU). Por razones de rendimiento, una implementación vectorizable es, por lo tanto, altamente deseable. Esto significa que se deben evitar las ramas, excepto como parte de la asignación de selección. Del mismo modo, el uso extensivo de tablas no está indicado, ya que la búsqueda en paralelo es a menudo ineficiente.

¿Cómo se podría construir una implementación erfcf() eficiente de la función de precisión erfcf() ? La precisión, medida en ulp , debe ser aproximadamente la misma que la implementación escalar de glibc, que tiene un error máximo de 3.12575 ulps (determinado por pruebas exhaustivas). Se puede suponer la disponibilidad de Fused Multiplex-Add (Fused Multiplex-Add), ya que todas las arquitecturas de procesadores principales (CPU y GPU) lo ofrecen en este momento. Si bien se pueden ignorar los indicadores de estado de coma flotante y errno , los denormales, infinitos y NaN se deben manejar de acuerdo con los enlaces IEEE 754 para ISO C.

Answer 1

Después de examinar varios enfoques, el que parece más adecuado es el algoritmo propuesto en el siguiente documento:

MM Shepherd y JG Laframboise, "Aproximación de Chebyshev de (1 + 2 x ) exp ( x ² ) erfc x en 0 ≤ x <∞." Matemáticas de la computación , volumen 36, n.º 153, enero de 1981, págs. 249-253 ( copia en línea )

La idea básica del trabajo es crear una aproximación a (1 + 2 x ) exp ( x ² ) erfc ( x ), a partir de la cual podemos calcular erfcx ( x ) simplemente dividiendo por (1 + 2 x ), y erfc (x) multiplicándose luego con exp (- x ² ). El rango estrechamente delimitado de la función, con valores de función aproximadamente en [1, 1.3], y su "planeidad" general se presta bien a la aproximación polinómica. Las propiedades numéricas de este enfoque se mejoran aún más estrechando el intervalo de aproximación: el argumento original x se transforma por q = ( x - K) / ( x + K), donde K es una constante adecuadamente elegida, seguida de la computación p ( q ) , donde p es un polinomio.

Como erfc - x = 2 - erfc x , solo tenemos que considerar el intervalo [0, ∞] que está mapeado al intervalo [-1, 1] mediante esta transformación. Para la precisión simple IEEE-754, erfcf() desaparece (se convierte en cero) para x > 10.0546875, por lo que solo se debe considerar x ∈ [0, 10.0546875]. ¿Cuál es el valor "óptimo" de K para este rango? No conozco ningún análisis matemático que proporcione la respuesta, el documento sugiere K = 3.75 basado en experimentos.

Se puede establecer fácilmente que para el cálculo de precisión simple, una aproximación polinomial minimax de grado 9 es suficiente para varios valores de K en esa vecindad general. Generando sistemáticamente tales aproximaciones con el algoritmo Remez, con K variando entre 1.5 y 4 en pasos de 1/16, se observa un error de aproximación más bajo para K = {2, 2.625, 3.3125}. De estos, K = 2 es la opción más ventajosa, ya que se presta a un cálculo muy preciso de ( x - K) / ( x + K), como se muestra en esta pregunta .

El valor K = 2 y el dominio de entrada para x sugerirían que es necesario usar la variante 4 de mi respuesta , sin embargo, una vez puede demostrarse experimentalmente que la variante menos costosa 5 logra la misma precisión aquí, lo que probablemente se debe a la muy superficial pendiente de la función aproximada para q > -0.5, lo que hace que cualquier error en el argumento q se reduzca en aproximadamente un factor de diez.

Dado que el cálculo de erfc() requiere pasos posteriores al procesamiento además de la aproximación inicial, es claro que la precisión de estos dos cálculos debe ser alta para lograr un resultado final suficientemente preciso. Se deben usar técnicas de corrección de errores.

Se observa que el coeficiente más significativo en la aproximación polinómica de (1 + 2 x ) exp ( x ² ) erfc ( x ) es de la forma (1 + s), donde s <0,5. Esto significa que podemos representar el coeficiente principal de forma más precisa dividiendo 1, y solo usando s en el polinomio. Entonces, en lugar de calcular un polinomio p (q), multiplicar por el recíproco r = 1 / (1 + 2 x ), es matemáticamente equivalente pero numéricamente ventajoso calcular la aproximación del núcleo como p ( q ) + 1, y usar p para calcular fma (p, r, r) .

La precisión de la división se puede mejorar calculando un cociente inicial q desde el recíproco r , calcule el residuo e = p +1 - q * (1 + 2 x ) con la ayuda de un FMA, luego use e para aplicar la corrección q = q + ( e * r ), nuevamente usando un FMA.

La exponenciación tiene propiedades de ampliación de error, por lo tanto, el cálculo de e ^{- x 2} debe realizarse con cuidado. La disponibilidad de FMA trivialmente permite el cálculo de - x ² como un doble float s _alto : s _bajo . e ^x es su propia derivada, por lo que uno puede calcular e ^{s _alto : s _bajo} como e ^{s _alto} + e ^{s _alto} * s _bajo . Este cálculo se puede combinar con la multiplicación del resultado intermedio previo r para dar r = r * e ^{s _alto} + r * e ^{s _alto} * s _bajo . Mediante el uso de FMA, uno se asegura de que el término más significativo r * e ^{s _alto} se calcule con la mayor precisión posible.

Combinando los pasos anteriores con algunas selecciones simples para manejar casos excepcionales y argumentos negativos, se llega al siguiente código C:

float my_expf (float);

/*  
 * Based on: M. M. Shepherd and J. G. Laframboise, "Chebyshev Approximation of 
 * (1+2x)exp(x^2)erfc x in 0 <= x < INF", Mathematics of Computation, Vol. 36,
 * No. 153, January 1981, pp. 249-253.  
 */  
float my_erfcf (float x)
{
    float a, d, e, m, p, q, r, s, t;

    a = fabsf (x); 

    /* Compute q = (a-2)/(a+2) accurately. [0, 10.0546875] -> [-1, 0.66818] */
    m = a - 2.0f;
    p = a + 2.0f;
    r = 1.0f / p;
    q = m * r;
    t = fmaf (q + 1.0f, -2.0f, a); 
    e = fmaf (q, -a, t); 
    q = fmaf (r, e, q); 

    /* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1, 0.66818] */
    p =             -0x1.a48024p-12f;  // -4.01020574e-4
    p = fmaf (p, q, -0x1.42a172p-10f); // -1.23073824e-3
    p = fmaf (p, q,  0x1.585784p-10f); //  1.31355994e-3
    p = fmaf (p, q,  0x1.1ade24p-07f); //  8.63243826e-3
    p = fmaf (p, q, -0x1.081b72p-07f); // -8.05991236e-3
    p = fmaf (p, q, -0x1.bc0b94p-05f); // -5.42047396e-2
    p = fmaf (p, q,  0x1.4ffc40p-03f); //  1.64055347e-1
    p = fmaf (p, q, -0x1.540840p-03f); // -1.66031361e-1
    p = fmaf (p, q, -0x1.7bf612p-04f); // -9.27639678e-2
    p = fmaf (p, q,  0x1.1ba03ap-02f); //  2.76978403e-1

    /* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
    t = a + a;
    d = t + 1.0f;
    r = 1.0f / d;
    q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
    e = (p - q) + fmaf (q, -t, 1.0f); // (p+1) - q*(1+2*a)
    r = fmaf (e, r, q);

    /* Multiply by exp(-a*a) ==> erfc(a) */
    s = a * a; 
    e = my_expf (-s);  
    t = fmaf (a, -a, s);
    r = fmaf (r, e, r * e * t);

    /* Handle NaN arguments to erfc() */
    if (!(a <= 0x1.fffffep127f)) r = x + x;

    /* Clamp result for large arguments */
    if (a > 10.0546875f) r = 0.0f;

    /* Handle negative arguments to erfc() */
    if (x < 0.0f) r = 2.0f - r; 

    return r;
}

/* Compute exponential base e. Maximum ulp error = 0.87161 */
float my_expf (float a)
{
    float c, f, r;
    int i;

    // exp(a) = exp(i + f); i = rint (a / log(2)) 
    c = 0x1.800000p+23f; // 1.25829120e+7
    r = fmaf (0x1.715476p+0f, a, c) - c; // 1.44269502e+0
    f = fmaf (r, -0x1.62e400p-01f, a); // -6.93145752e-1 // log_2_hi 
    f = fmaf (r, -0x1.7f7d1cp-20f, f); // -1.42860677e-6 // log_2_lo
    i = (int)r;
    // approximate r = exp(f) on interval [-log(2)/2,+log(2)/2]
    r =             0x1.6a98dap-10f;  // 1.38319808e-3
    r = fmaf (r, f, 0x1.1272cap-07f); // 8.37550033e-3
    r = fmaf (r, f, 0x1.555a20p-05f); // 4.16689515e-2
    r = fmaf (r, f, 0x1.55542ep-03f); // 1.66664466e-1
    r = fmaf (r, f, 0x1.fffff6p-02f); // 4.99999851e-1
    r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
    r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
    // exp(a) = 2**i * exp(f);
    r = ldexpf (r, i);
    // handle special cases
    if (!(fabsf (a) < 104.0f)) {
        r = a + a; // handle NaNs
        if (a < 0.0f) r = 0.0f;
        if (a > 0.0f) r = 1e38f * 1e38f; // + INF
    }
    return r;
}

expf() mi propia implementación de expf() en el código anterior para aislar mi trabajo de las diferencias en las implementaciones de expf() en diferentes plataformas de cómputo. Pero cualquier implementación de expf() cuyo error máximo sea cercano a 0.5 ulp debería funcionar bien. Como se muestra arriba, es decir, cuando se usa my_expf() , my_erfcf() tiene un error máximo de 2.65712 ulps. Siempre que haya disponibilidad de un expf() vectorizable, el código anterior debería vectorizarse sin problema. Hice una comprobación rápida con el compilador de Intel 13.1.3.198. Puse una llamada a my_erfcf() en un bucle, agregué #include <mathimf.h> , reemplacé la llamada a my_expf() con una llamada a expf() , luego compilé usando estos modificadores de línea de comando:

/Qstd=c99 /O3 /QxCORE-AVX2 /fp:precise /Qfma /Qimf-precision:high:expf /Qvec_report=2

El compilador de Intel informó que el bucle había sido vectorizado, lo que comprobé dos veces mediante la inspección del código binario desensamblado.

Como my_erfcf() solo usa recíprocos en lugar de divisiones completas, es susceptible de uso de implementaciones recíprocas rápidas, siempre que entreguen resultados casi correctamente redondeados. Para los procesadores que proporcionan una aproximación recíproca rápida de precisión simple en hardware, esto puede lograrse fácilmente acoplando esto con una iteración de Halley con convergencia cúbica. Un ejemplo (escalar) de este enfoque para los procesadores x86 es:

/* Compute 1.0f / a almost correctly rounded. Halley iteration with cubic convergence */
float fast_recipf (float a)
{
    __m128 t;
    float e, r;
    t = _mm_set_ss (a);
    t = _mm_rcp_ss (t);
    _mm_store_ss (&r, t);
    e = fmaf (r, -a, 1.0f);
    e = fmaf (e, e, e);
    r = fmaf (e, r, r);
    return r;
}