menos - La mejor aproximación de minimax polinómica optimizada para máquina a arcotangente en[-1,1]?

tangente a la menos 1

Question

Para la implementación simple y eficiente de las funciones matemáticas rápidas con una precisión razonable, las aproximaciones polinómicas minimax son a menudo el método de elección. Las aproximaciones mínimas se generan generalmente con una variante del algoritmo de Remez. Varias herramientas ampliamente disponibles como Maple y Mathematica tienen funcionalidad incorporada para esto. Los coeficientes generados generalmente se calculan usando aritmética de alta precisión. Es bien sabido que simplemente el redondeo de esos coeficientes a la precisión de la máquina conduce a una precisión subóptima en la implementación resultante.

En cambio, uno busca conjuntos de coeficientes estrechamente relacionados que son exactamente representables como números de máquina para generar una aproximación optimizada para la máquina. Dos documentos relevantes son:

Nicolas Brisebarre, Jean-Michel Muller y Arnaud Tisserand, "Aproximaciones polinomiales eficientes en el cálculo de la máquina", ACM Transactions on Mathematical Software, vol. 32, No. 2, junio de 2006, págs. 236-256.

Nicolas Brisebarre y Sylvain Chevillard, "Eficaces aproximaciones L∞ polinómicas", 18º Simposio IEEE sobre Aritmética Computacional (ARITH-18), Montpellier (Francia), junio de 2007, pp. 169-176.

Una implementación del algoritmo LLL del último documento está disponible como el fpminimax() de la herramienta Sollya . Tengo entendido que todos los algoritmos propuestos para la generación de aproximaciones optimizadas para la máquina se basan en la heurística, y que, por lo tanto, en general se desconoce qué precisión se puede lograr mediante una aproximación óptima. No me queda claro si la disponibilidad de FMA (fusionado multiplicar-agregar) para la evaluación de la aproximación influye en la respuesta a esa pregunta. Me parece ingenuamente que debería.

Actualmente estoy buscando una aproximación polinomial simple para arcotangente en [-1,1] que se evalúa en aritmética de precisión simple IEEE-754, usando el esquema de Horner y FMA. Ver la función atan_poly() en el código C99 a continuación. Debido a la falta de acceso a una máquina Linux en este momento, no utilicé Sollya para generar estos coeficientes, pero utilicé mi propia heurística que podría describirse como una mezcla del recocido más decente y simulado (para evitar quedar atrapado en los mínimos locales) . El error máximo de mi polinomio optimizado para máquina es muy cercano a 1 ulp, pero idealmente me gustaría que el error máximo de ulp sea inferior a 1 ulp.

Soy consciente de que podría cambiar mi cálculo para aumentar la precisión, por ejemplo, utilizando un coeficiente principal representado a más de precisión de precisión simple, pero me gustaría mantener el código exactamente como está (es decir, lo más simple posible) ajustando solo los coeficientes para entregar el resultado más preciso posible.

Un conjunto de coeficientes óptimos "probados" sería ideal; se agradecen los indicadores de la literatura pertinente. Realicé una búsqueda bibliográfica pero no pude encontrar ningún documento que avance significativamente el estado del arte más allá del fpminimax() de fpminimax() , y ninguno que examine el rol de las FMA (si corresponde) en este tema.

// max ulp err = 1.03143
float atan_poly (float a)
{
    float r, s;
    s = a * a;
    r =              0x1.7ed1ccp-9f;
    r = fmaf (r, s, -0x1.0c2c08p-6f);
    r = fmaf (r, s,  0x1.61fdd0p-5f);
    r = fmaf (r, s, -0x1.3556b2p-4f);
    r = fmaf (r, s,  0x1.b4e128p-4f);
    r = fmaf (r, s, -0x1.230ad2p-3f);
    r = fmaf (r, s,  0x1.9978ecp-3f);
    r = fmaf (r, s, -0x1.5554dcp-2f);
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}

// max ulp err = 1.52637
float my_atanf (float a)
{
    float r, t;
    t = fabsf (a);
    r = t;
    if (t > 1.0f) {
        r = 1.0f / r;
    }
    r = atan_poly (r);
    if (t > 1.0f) {
        r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r
    }
    r = copysignf (r, a);
    return r;
}

Answer 1

Esta no es una respuesta a la pregunta, pero es demasiado larga para caber en un comentario:

su pregunta es acerca de la elección óptima de los coeficientes C3, C5, ..., C17 en una aproximación polinómica a arcotangente donde ha fijado C1 a 1 y C2, C4, ..., C16 a 0.

El título de su pregunta dice que está buscando aproximaciones en [-1, 1], y una buena razón para anclar los coeficientes pares a 0 es que es suficiente y necesario para que la aproximación sea exactamente una función impar. El código en su pregunta "contradice" el título al aplicar la aproximación polinómica solo en [0, 1].

Si usas el algoritmo de Remez para buscar los coeficientes C2, C3, ..., C8 en una aproximación polinomial de arcotangente en [0, 1], puedes terminar con algo como los siguientes valores:

#include <stdio.h>
#include <math.h>

float atan_poly (float a)
{
  float r, s;
  s = a;
  //  s = a * a;

  r =             -3.3507930064626076153585890630056286726807491543578e-2;
  r = fmaf (r, s, 1.3859776280052980081098065189344699108643282883702e-1);
  r = fmaf (r, s, -1.8186361916440430105127602496688553126414578766147e-1);
  r = fmaf (r, s, -1.4583047494913656326643327729704639191810926020847e-2);
  r = fmaf (r, s, 2.1335202878219865228365738728594741358740655881373e-1);
  r = fmaf (r, s, -3.6801711826027841250774413728610805847547996647342e-3);
  r = fmaf (r, s, -3.3289852243978319173749528028057608377028846413080e-1);
  r = fmaf (r, s, -1.8631479933914856903459844359251948006605218562283e-5);
  r = fmaf (r, s, 1.2917291732886065585264586294461539492689296797761e-7);

  r = fmaf (r, a, a);
  return r;
}

int main() {
  for (float x = 0.0f; x < 1.0f; x+=0.1f)
    printf("x: %f/n%a/n%a/n/n", x, atan_poly(x), atan(x));
}

Esto tiene aproximadamente la misma complejidad que el código en su pregunta: el número de multiplicaciones es similar. En cuanto a este polinomio, no hay ninguna razón en particular para querer fijar cualquier coeficiente a 0. Si quisiéramos aproximar una función impar sobre [-1, 1] sin fijar los coeficientes pares, automáticamente aparecerían muy pequeños y sujetos a la absorción, y luego querríamos fijarlos a 0, pero para esta aproximación sobre [0, 1], no es así, así que no tenemos que fijarlos.

Pudo haber sido mejor o peor que un polinomio impar en tu pregunta. Resulta que es peor (ver abajo). Esta rápida y sucia aplicación de LolRemez 0.2 (código incluido en la parte inferior de la pregunta) parece ser, sin embargo, lo suficientemente buena como para plantear la cuestión de la elección de los coeficientes. En particular, me gustaría saber qué sucede si se someten los coeficientes de esta respuesta a la misma etapa de optimización de "mezcla de mayor pendiente y recocido simulado" que aplicó para obtener los coeficientes en su pregunta.

Entonces, para resumir esta observación publicada como respuesta, ¿está seguro de que está buscando los coeficientes óptimos C3, C5, ..., C17? Me parece que estás buscando la mejor secuencia de operaciones de punto flotante de precisión simple que produzca una aproximación fiel a arcotangente, y que esta aproximación no tiene que ser la forma de Horner de un polinomio impar de grado 17.

x: 0.000000 0x0p+0 0x0p+0 x: 0.100000 0x1.983e2cp-4 0x1.983e28938f9ecp-4 x: 0.200000 0x1.94442p-3 0x1.94441ff1e8882p-3 x: 0.300000 0x1.2a73a6p-2 0x1.2a73a71dcec16p-2 x: 0.400000 0x1.85a37ap-2 0x1.85a3770ebe7aep-2 x: 0.500000 0x1.dac67p-2 0x1.dac670561bb5p-2 x: 0.600000 0x1.14b1dcp-1 0x1.14b1ddf627649p-1 x: 0.700000 0x1.38b116p-1 0x1.38b113eaa384ep-1 x: 0.800000 0x1.5977a8p-1 0x1.5977a686e0ffbp-1 x: 0.900000 0x1.773388p-1 0x1.77338c44f8faep-1

Este es el código que vinculé a LolRemez 0.2 para optimizar la precisión relativa de una aproximación polinómica de grado 9 de arcotangente en [0, 1]:

#include "lol/math/real.h"
#include "lol/math/remez.h"

using lol::real;
using lol::RemezSolver;

real f(real const &y)
{
  return (atan(y) - y) / y;
}

real g(real const &y)
{
  return re (atan(y) / y);
}

int main(int argc, char **argv)
{
  RemezSolver<8, real> solver;
  solver.Run("1e-1000", 1.0, f, g, 50);
  return 0;
}

Answer 2

Esta no es una respuesta, sino un comentario extendido también.

Las CPU Intel recientes y algunas futuras CPU AMD tienen AVX2. En Linux, busque el indicador de avx2 en /proc/cpuinfo para ver si su CPU los admite.

AVX2 es una extensión que nos permite construir y calcular utilizando vectores de 256 bits, por ejemplo, ocho números de precisión simple o cuatro números de precisión doble, en lugar de solo escalares. Incluye compatibilidad con FMA3, es decir, fusionado y multiplicado para esos vectores. En pocas palabras, AVX2 nos permite evaluar ocho polinomios en paralelo, casi al mismo tiempo que evaluamos uno solo usando operaciones escalares.

La función error8() analiza un conjunto de coeficientes, usando valores predefinidos de x , comparando contra valores precalculados de atan(x) , y devuelve el error en ULP (por debajo y por encima del resultado deseado por separado), así como el número de resultados que coinciden exactamente con el valor de coma flotante deseado. Estos no son necesarios para simplemente probar si un conjunto de coeficientes es mejor que el conjunto actualmente mejor conocido, pero permiten diferentes estrategias sobre qué coeficientes probar. (Básicamente, el error máximo en ULP forma una superficie, y estamos tratando de encontrar el punto más bajo en esa superficie, conocer la "altura" de la superficie en cada punto nos permite hacer conjeturas sobre qué dirección tomar - - cómo cambiar los coeficientes.)

Hay cuatro tablas precalculadas utilizadas: known_x para los argumentos, known_f para los resultados de precisión simple redondeados correctamente, known_a para el valor "preciso" de doble precisión (solo espero que la biblioteca atan() sea ​​lo suficientemente precisa para esto - - ¡pero no se debe confiar en él sin verificarlo!), y known_m para escalar la diferencia de doble precisión con los ULP. Dado un rango de argumentos deseado, la función precalculate() precalculará estos utilizando la función de biblioteca atan() . (También se basa en que los formatos de coma flotante IEEE-754 y el orden de bytes flotante y entero son los mismos, pero esto es cierto en las CPU en las que se ejecuta este código).

Tenga en cuenta que las known_x , known_f , y known_a se pueden almacenar en archivos binarios; los contenidos de known_m se derivan trivialmente de known_a . Usar la biblioteca atan() sin verificar que no sea una buena idea, pero como los resultados de la mía coinciden con los de njuffa, no me molesté en buscar una referencia mejor atan() .

Para simplificar, aquí está el código en forma de un programa de ejemplo:

#define _POSIX_C_SOURCE 200809L
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <immintrin.h>
#include <math.h>
#include <errno.h>

/** poly8() - Compute eight polynomials in parallel.
 * @x - the arguments
 * @c - the coefficients.
 *
 * The first coefficients are for degree 17, the second
 * for degree 15, and so on, down to degree 3.
 *
 * The compiler should vectorize the expression using vfmaddXXXps
 * given an AVX2-capable CPU; for example, Intel Haswell,
 * Broadwell, Haswell E, Broadwell E, Skylake, or Cannonlake;
 * or AMD Excavator CPUs. Tested on Intel Core i5-4200U.
 *
 * Using GCC-4.8.2 and
 *     gcc -O2 -march=core-avx2 -mtune=generic
 * this code produces assembly (AT&T syntax)
 *     vmulps       %ymm0, %ymm0, %ymm2
 *     vmovaps      (%rdi), %ymm1
 *     vmovaps      %ymm0, %ymm3
 *     vfmadd213ps  32(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  64(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  96(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  128(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  160(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  192(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  224(%rdi), %ymm2, %ymm1
 *     vmulps       %ymm2, %ymm1, %ymm0
 *     vfmadd132ps  %ymm3, %ymm3, %ymm0
 *     ret
 * if you omit the ''static inline''.
*/
static inline __v8sf poly8(const __v8sf x, const __v8sf *const c)
{
    const __v8sf xx = x * x;
    return (((((((c[0]*xx + c[1])*xx + c[2])*xx + c[3])*xx + c[4])*xx + c[5])*xx + c[6])*xx + c[7])*xx*x + x;
}

/** error8() - Calculate maximum error in ULPs
 * @x  - the arguments
 * @co - { C17, C15, C13, C11, C9, C7, C5, C3 }
 * @f  - the correctly rounded results in single precision
 * @a  - the expected results in double precision
 * @m  - 16777216.0 raised to the same power of two as @a normalized
 * @n  - number of vectors to test
 * @max_under - pointer to store the maximum underflow (negative, in ULPs) to
 * @max_over  - pointer to store the maximum overflow (positive, in ULPs) to
 * Returns the number of correctly rounded float results, 0..8*n.
*/
size_t error8(const __v8sf *const x, const float *const co,
              const __v8sf *const f, const __v4df *const a, const __v4df *const m,
              const size_t n,
              float *const max_under, float *const max_over)
{
    const __v8sf c[8] = { { co[0], co[0], co[0], co[0], co[0], co[0], co[0], co[0] },
                          { co[1], co[1], co[1], co[1], co[1], co[1], co[1], co[1] },
                          { co[2], co[2], co[2], co[2], co[2], co[2], co[2], co[2] },
                          { co[3], co[3], co[3], co[3], co[3], co[3], co[3], co[3] },
                          { co[4], co[4], co[4], co[4], co[4], co[4], co[4], co[4] },
                          { co[5], co[5], co[5], co[5], co[5], co[5], co[5], co[5] },
                          { co[6], co[6], co[6], co[6], co[6], co[6], co[6], co[6] },
                          { co[7], co[7], co[7], co[7], co[7], co[7], co[7], co[7] } };
    __v4df min = { 0.0, 0.0, 0.0, 0.0 };
    __v4df max = { 0.0, 0.0, 0.0, 0.0 };
    __v8si eqs = { 0, 0, 0, 0, 0, 0, 0, 0 };
    size_t i;

    for (i = 0; i < n; i++) {
        const __v8sf v = poly8(x[i], c);
        const __v4df d0 = { v[0], v[1], v[2], v[3] };
        const __v4df d1 = { v[4], v[5], v[6], v[7] };
        const __v4df err0 = (d0 - a[2*i+0]) * m[2*i+0];
        const __v4df err1 = (d1 - a[2*i+1]) * m[2*i+1];
        eqs -= (__v8si)_mm256_cmp_ps(v, f[i], _CMP_EQ_OQ);
        min = _mm256_min_pd(min, err0);
        max = _mm256_max_pd(max, err1);
        min = _mm256_min_pd(min, err1);
        max = _mm256_max_pd(max, err0);
    }

    if (max_under) {
        if (min[0] > min[1]) min[0] = min[1];
        if (min[0] > min[2]) min[0] = min[2];
        if (min[0] > min[3]) min[0] = min[3];
        *max_under = min[0];
    }

    if (max_over) {
        if (max[0] < max[1]) max[0] = max[1];
        if (max[0] < max[2]) max[0] = max[2];
        if (max[0] < max[3]) max[0] = max[3];
        *max_over = max[0];
    }

    return (size_t)((unsigned int)eqs[0])
         + (size_t)((unsigned int)eqs[1])
         + (size_t)((unsigned int)eqs[2])
         + (size_t)((unsigned int)eqs[3])
         + (size_t)((unsigned int)eqs[4])
         + (size_t)((unsigned int)eqs[5])
         + (size_t)((unsigned int)eqs[6])
         + (size_t)((unsigned int)eqs[7]);
}

/** precalculate() - Allocate and precalculate tables for error8().
 * @x0   - First argument to precalculate
 * @x1   - Last argument to precalculate
 * @xptr - Pointer to a __v8sf pointer for the arguments
 * @fptr - Pointer to a __v8sf pointer for the correctly rounded results
 * @aptr - Pointer to a __v4df pointer for the comparison results
 * @mptr - Pointer to a __v4df pointer for the difference multipliers
 * Returns the vector count if successful,
 * 0 with errno set otherwise.
*/
size_t precalculate(const float x0, const float x1,
                 __v8sf **const xptr, __v8sf **const fptr,
                 __v4df **const aptr, __v4df **const mptr)
{
    const size_t align = 64;
    unsigned int i0, i1;
    size_t       n, i, sbytes, dbytes;
    __v8sf      *x = NULL;
    __v8sf      *f = NULL;
    __v4df      *a = NULL;
    __v4df      *m = NULL;

    if (!xptr || !fptr || !aptr || !mptr) {
        errno = EINVAL;
        return (size_t)0;
    }

    memcpy(&i0, &x0, sizeof i0);
    memcpy(&i1, &x1, sizeof i1);

    i0 ^= (i0 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U;
    i1 ^= (i1 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U;

    if (i1 > i0)
        n = (((size_t)i1 - (size_t)i0) | (size_t)7) + (size_t)1;
    else
    if (i0 > i1)
        n = (((size_t)i0 - (size_t)i1) | (size_t)7) + (size_t)1;
    else {
        errno = EINVAL;
        return (size_t)0;
    }

    sbytes = n * sizeof (float);
    if (sbytes % align)
        sbytes += align - (sbytes % align);

    dbytes = n * sizeof (double);
    if (dbytes % align)
        dbytes += align - (dbytes % align);

    if (posix_memalign((void **)&x, align, sbytes)) {
        errno = ENOMEM;
        return (size_t)0;
    }
    if (posix_memalign((void **)&f, align, sbytes)) {
        free(x);
        errno = ENOMEM;
        return (size_t)0;
    }
    if (posix_memalign((void **)&a, align, dbytes)) {
        free(f);
        free(x);
        errno = ENOMEM;
        return (size_t)0;
    }
    if (posix_memalign((void **)&m, align, dbytes)) {
        free(a);
        free(f);
        free(x);
        errno = ENOMEM;
        return (size_t)0;
    }

    if (x1 > x0) {
        float *const xp = (float *)x;
        float        curr = x0;

        for (i = 0; i < n; i++) {
            xp[i] = curr;
            curr = nextafterf(curr, HUGE_VALF);
        }

        i = n;
        while (i-->0 && xp[i] > x1)
            xp[i] = x1;
    } else {
        float *const xp = (float *)x;
        float        curr = x0;

        for (i = 0; i < n; i++) {
            xp[i] = curr;
            curr = nextafterf(curr, -HUGE_VALF);
        }

        i = n;
        while (i-->0 && xp[i] < x1)
            xp[i] = x1;
    }

    {
        const float *const xp = (const float *)x;
        float *const       fp = (float *)f;
        double *const      ap = (double *)a;
        double *const      mp = (double *)m;

        for (i = 0; i < n; i++) {
            const float curr = xp[i];
            int         temp;

            fp[i] = atanf(curr);
            ap[i] = atan((double)curr);

            (void)frexp(ap[i], &temp);
            mp[i] = ldexp(16777216.0, temp);
        }
    }

    *xptr = x;
    *fptr = f;
    *aptr = a;
    *mptr = m;

    errno = 0;
    return n/8;
}

static int parse_range(const char *const str, float *const range)
{
    float fmin, fmax;
    char  dummy;

    if (sscanf(str, " %f %f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %f:%f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %f,%f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %f/%f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff %ff %c", &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff:%ff %c", &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff,%ff %c", &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff/%ff %c", &fmin, &fmax, &dummy) == 2) {
        if (range) {
            range[0] = fmin;
            range[1] = fmax;
        }
        return 0;
    }

    if (sscanf(str, " %f %c",  &fmin, &dummy) == 1 ||
        sscanf(str, " %ff %c", &fmin, &dummy) == 1) {
        if (range) {
            range[0] = fmin;
            range[1] = fmin;
        }
        return 0;
    }

    return errno = ENOENT;
}

static int fix_range(float *const f)
{
    if (f && f[0] > f[1]) {
        const float tmp = f[0];
        f[0] = f[1];
        f[1] = tmp;
    }
    return f && isfinite(f[0]) && isfinite(f[1]) && (f[1] >= f[0]);
}

static const char *f2s(char *const buffer, const size_t size, const float value, const char *const invalid)
{
    char  format[32];
    float parsed;
    int   decimals, length;

    for (decimals = 0; decimals <= 16; decimals++) {
        length = snprintf(format, sizeof format, "%%.%df", decimals);
        if (length < 1 || length >= (int)sizeof format)
            break;

        length = snprintf(buffer, size, format, value);
        if (length < 1 || length >= (int)size)
            break;

        if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value)
            return buffer;

        decimals++;
    }

    for (decimals = 0; decimals <= 16; decimals++) {
        length = snprintf(format, sizeof format, "%%.%dg", decimals);
        if (length < 1 || length >= (int)sizeof format)
            break;

        length = snprintf(buffer, size, format, value);
        if (length < 1 || length >= (int)size)
            break;

        if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value)
            return buffer;

        decimals++;
    }

    length = snprintf(buffer, size, "%a", value);
    if (length < 1 || length >= (int)size)
        return invalid;

    if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value)
        return buffer;

    return invalid;
}

int main(int argc, char *argv[])
{
    float xrange[2] = { 0.75f, 1.00f };
    float c17range[2], c15range[2], c13range[2], c11range[2];
    float c9range[2], c7range[2], c5range[2], c3range[2];
    float c[8];

    __v8sf *known_x;
    __v8sf *known_f;
    __v4df *known_a;
    __v4df *known_m;
    size_t  known_n;

    if (argc != 10 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
        fprintf(stderr, "/n");
        fprintf(stderr, "Usage: %s [ -h | --help ]/n", argv[0]);
        fprintf(stderr, "       %s C17 C15 C13 C11 C9 C7 C5 C3 x/n", argv[0]);
        fprintf(stderr, "/n");
        fprintf(stderr, "Each of the coefficients can be a constant or a range,/n");
        fprintf(stderr, "for example 0.25 or 0.75:1. x must be a non-empty range./n");
        fprintf(stderr, "/n");
        return EXIT_FAILURE;
    }

    if (parse_range(argv[1], c17range) || !fix_range(c17range)) {
        fprintf(stderr, "%s: Invalid C17 range or constant./n", argv[1]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[2], c15range) || !fix_range(c15range)) {
        fprintf(stderr, "%s: Invalid C15 range or constant./n", argv[2]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[3], c13range) || !fix_range(c13range)) {
        fprintf(stderr, "%s: Invalid C13 range or constant./n", argv[3]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[4], c11range) || !fix_range(c11range)) {
        fprintf(stderr, "%s: Invalid C11 range or constant./n", argv[4]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[5], c9range) || !fix_range(c9range)) {
        fprintf(stderr, "%s: Invalid C9 range or constant./n", argv[5]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[6], c7range) || !fix_range(c7range)) {
        fprintf(stderr, "%s: Invalid C7 range or constant./n", argv[6]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[7], c5range) || !fix_range(c5range)) {
        fprintf(stderr, "%s: Invalid C5 range or constant./n", argv[7]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[8], c3range) || !fix_range(c3range)) {
        fprintf(stderr, "%s: Invalid C3 range or constant./n", argv[8]);
        return EXIT_FAILURE;
    }

    if (parse_range(argv[9], xrange) || xrange[0] == xrange[1] ||
        !isfinite(xrange[0]) || !isfinite(xrange[1])) {
        fprintf(stderr, "%s: Invalid x range./n", argv[9]);
        return EXIT_FAILURE;
    }

    known_n = precalculate(xrange[0], xrange[1], &known_x, &known_f, &known_a, &known_m);
    if (!known_n) {
        if (errno == ENOMEM)
            fprintf(stderr, "Not enough memory for precalculated tables./n");
        else
            fprintf(stderr, "Invalid (empty) x range./n");
        return EXIT_FAILURE;
    }

    fprintf(stderr, "Precalculated %lu arctangents to compare to./n", 8UL * (unsigned long)known_n);
    fprintf(stderr, "/nC17 C15 C13 C11 C9 C7 C5 C3 max-ulps-under max-ulps-above correctly-rounded percentage cycles/n");
    fflush(stderr);

    {
        const double  percent = 12.5 / (double)known_n;
        size_t        rounded;
        char          c17buffer[64], c15buffer[64], c13buffer[64], c11buffer[64];
        char          c9buffer[64], c7buffer[64], c5buffer[64], c3buffer[64];
        char          minbuffer[64], maxbuffer[64];
        float         minulps, maxulps;
        unsigned long tsc_start, tsc_stop;

        for (c[0] = c17range[0]; c[0] <= c17range[1]; c[0] = nextafterf(c[0], HUGE_VALF))
        for (c[1] = c15range[0]; c[1] <= c15range[1]; c[1] = nextafterf(c[1], HUGE_VALF))
        for (c[2] = c13range[0]; c[2] <= c13range[1]; c[2] = nextafterf(c[2], HUGE_VALF))
        for (c[3] = c11range[0]; c[3] <= c11range[1]; c[3] = nextafterf(c[3], HUGE_VALF))
        for (c[4] = c9range[0]; c[4] <= c9range[1]; c[4] = nextafterf(c[4], HUGE_VALF))
        for (c[5] = c7range[0]; c[5] <= c7range[1]; c[5] = nextafterf(c[5], HUGE_VALF))
        for (c[6] = c5range[0]; c[6] <= c5range[1]; c[6] = nextafterf(c[6], HUGE_VALF))
        for (c[7] = c3range[0]; c[7] <= c3range[1]; c[7] = nextafterf(c[7], HUGE_VALF)) {
            tsc_start = __builtin_ia32_rdtsc();
            rounded = error8(known_x, c, known_f, known_a, known_m, known_n, &minulps, &maxulps);
            tsc_stop = __builtin_ia32_rdtsc();
            printf("%-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %lu %.3f %lu/n",
                   f2s(c17buffer, sizeof c17buffer, c[0], "?"),
                   f2s(c15buffer, sizeof c15buffer, c[1], "?"),
                   f2s(c13buffer, sizeof c13buffer, c[2], "?"),
                   f2s(c11buffer, sizeof c11buffer, c[3], "?"),
                   f2s(c9buffer, sizeof c9buffer, c[4], "?"),
                   f2s(c7buffer, sizeof c7buffer, c[5], "?"),
                   f2s(c5buffer, sizeof c5buffer, c[6], "?"),
                   f2s(c3buffer, sizeof c3buffer, c[7], "?"),
                   f2s(minbuffer, sizeof minbuffer, minulps, "?"),
                   f2s(maxbuffer, sizeof maxbuffer, maxulps, "?"),
                   rounded, (double)rounded * percent,
                   (unsigned long)(tsc_stop - tsc_start));
            fflush(stdout);

        }
    }

    return EXIT_SUCCESS;
}

El código se compila usando GCC-4.8.2 en Linux, pero podría tener que ser modificado para otros compiladores y / o sistemas operativos. (Sería feliz para incluir / aceptar las ediciones fijación de ellos, sin embargo. Simplemente no tienen ventanas o ICC, así que yo podía comprobar.)

Para compilar esto, recomiendo

gcc -Wall -O3 -fomit-frame-pointer -march=native -mtune=native example.c -lm -o example

Ejecutar sin argumentos para ver el uso; o

./example 0x1.7ed24ap-9f -0x1.0c2c12p-6f  0x1.61fdd2p-5f -0x1.3556b0p-4f 0x1.b4e138p-4f -0x1.230ae2p-3f  0x1.9978eep-3f -0x1.5554dap-2f 0.75:1

para comprobar lo que informa para el grupo coeficiente de njuffa, en comparación contra biblioteca estándar C atan()función, con todos los posibles x en [0,75, 1] considerado.

En lugar de un coeficiente fijo, también se puede utilizar min:maxpara definir un rango para escanear (analiza todo único de precisión simple valores de punto flotante). Cada combinación posible de los coeficientes se prueba.

Porque prefiero notación decimal, pero necesito mantener los valores exactos, utilizo la f2s()función para mostrar los valores de punto flotante. Es una simple función de fuerza bruta ayudante, que utiliza el formato más corto que da el mismo valor cuando analizado de nuevo a flote.

Por ejemplo,

./example 0x1.7ed248p-9f:0x1.7ed24cp-9f -0x1.0c2c10p-6f:-0x1.0c2c14p-6f 0x1.61fdd0p-5f:0x1.61fdd4p-5f -0x1.3556aep-4f:-0x1.3556b2p-4f 0x1.b4e136p-4f:0x1.b4e13ap-4f -0x1.230ae0p-3f:-0x1.230ae4p-3f 0x1.9978ecp-3f:0x1.9978f0p-3f -0x1.5554d8p-2f:-0x1.5554dcp-2f 0.75:1

calcula todos los 6561 (3 ⁸ ) combinaciones de coeficientes de ± 1 ULP alrededor conjunto de njuffa para x en [0,75, 1]. (De hecho, muestra que la disminución de C ₁₇ por 1 ULP a 0x1.7ed248p-9f produce los mismos resultados exactos.)

(Que se ejecutan tardó 90 segundos en Core i5-4200U a 2,6 GHz - bastante en línea en mi estimación de 30 conjuntos de coeficientes por segundo por GHz por núcleo Mientras que este código no está enhebrado, las funciones de las teclas son seguros para subprocesos, por lo. enhebrado no debería ser demasiado difícil. Esta Core i5-4200U es un ordenador portátil, y hace bastante calor incluso cuando haciendo hincapié en un solo núcleo, por lo que no se molestó.)

(Considero que el código anterior para estar en dominio público, o A0 con licencia, donde la dedicación de dominio público no es posible. De hecho, no estoy seguro si es lo suficientemente creativos para ser protegible en absoluto. De todos modos, no dude en usarlo en cualquier lugar, en cualquier forma que desee, siempre y cuando no me culpes si se rompe.)

¿Preguntas?Mejoras? Edita para fijar Linux / GCC''isms son bienvenidos!

Answer 3

La siguiente función es una implementación fielmente redondeada de arctan en [0, 1] :

float atan_poly (float a) {
  float s = a * a, u = fmaf(a, -a, 0x1.fde90cp-1f);
  float r1 =               0x1.74dfb6p-9f;
  float r2 = fmaf (r1, u,  0x1.3a1c7cp-8f);
  float r3 = fmaf (r2, s, -0x1.7f24b6p-7f);
  float r4 = fmaf (r3, u, -0x1.eb3900p-7f);
  float r5 = fmaf (r4, s,  0x1.1ab95ap-5f);
  float r6 = fmaf (r5, u,  0x1.80e87cp-5f);
  float r7 = fmaf (r6, s, -0x1.e71aa4p-4f);
  float r8 = fmaf (r7, u, -0x1.b81b44p-3f);
  float r9 = r8 * s;
  float r10 = fmaf (r9, a, a);
  return r10;
}

El arnés de prueba siguiente se cancelará si la función atan_poly no se redondea fielmente en [1e-16, 1] e imprime "éxito" en caso contrario:

int checkit(float f) {
  double d = atan(f);
  float d1 = d, d2 = d;
  if (d1 < d) d2 = nextafterf(d1, 1.0/0.0);
  else d1 = nextafterf(d1, -1.0/0.0);
  float p = atan_poly(f);
  if (p != d1 && p != d2) return 0;
  return 1;
}

int main() {
  for (float f = 1; f > 1e-16; f = nextafterf(f, -1.0/0.0)) {
    if (!checkit(f)) abort();
  }
  printf("success/n");
  exit(0);
}

El problema con el uso de s en cada multiplicación es que los coeficientes del polinomio no decaen rápidamente. Las entradas cercanas a 1 dan como resultado muchas y muchas cancelaciones de números casi iguales, lo que significa que está tratando de encontrar un conjunto de coeficientes para que el redondeo acumulado al final del cálculo se aproxime mucho al residual de arctan .

La constante 0x1.fde90cp-1f es un número cercano a 1 para el cual (arctan(sqrt(x)) - x) / x^3 está muy cerca del flotador más cercano. Es decir, es una constante que entra en el cálculo de u para que el coeficiente cúbico esté casi completamente determinado. (Para este programa, el coeficiente cúbico debe ser -0x1.b81b44p-3f o -0x1.b81b42p-3f ).

La alternancia de multiplicaciones por s u tiene el efecto de reducir el efecto del error de redondeo en ri sobre r{i+2} en un factor de como máximo 1/4, ya que s*u < 1/4 sea ​​lo que sea a. Esto le da un margen considerable para elegir los coeficientes de quinto orden y más allá.

Encontré los coeficientes con la ayuda de dos programas:

• Un programa conecta un grupo de puntos de prueba, escribe un sistema de desigualdades lineales y calcula los límites en los coeficientes de ese sistema de desigualdades. Observe que, dado a , uno puede calcular el rango de r8 que conduce a un resultado fielmente redondeado. Para obtener desigualdades lineales , fingí que r8 se computaría como un polinomio en el float s s y u en aritmética de números reales ; las desigualdades lineales limitaban este número real r8 a estar en algún intervalo. Usé la biblioteca de Polyhedra de Parma para manejar estos sistemas de restricción.
• Otro programa probó aleatoriamente conjuntos de coeficientes en ciertos rangos, conectando primero un conjunto de puntos de prueba y luego todos los float de 1 a 1e-8 en orden descendente y comprobando que atan_poly produce un redondeo fiel de atan((double)x) . Si alguna x falla, imprime esa x y por qué falló.

Para obtener coeficientes, pirateé este primer programa para corregir c3 , calcular límites en r7 para cada punto de prueba, y luego obtener límites en los coeficientes de orden superior. Luego lo pirateé para arreglar c3 y c5 y obtener límites en los coeficientes de orden superior. Hice esto hasta que tuve todos menos los tres coeficientes de orden más alto, c13 , c15 y c17 .

Crecí el conjunto de puntos de prueba en el segundo programa hasta que dejó de imprimir algo o imprimió "éxito". Necesitaba sorprendentemente pocos puntos de prueba para rechazar casi todos los polinomios incorrectos. Cuento 85 puntos de prueba en el programa.

Aquí muestro algo de mi trabajo seleccionando los coeficientes. Para obtener un arctan fielmente redondeado para mi conjunto inicial de puntos de prueba suponiendo que r1 a r8 se evalúan en aritmética real (y redondeados de alguna manera desagradable pero de una manera que no recuerdo) pero r9 y r10 se evalúan en aritmética float , Necesito:

-0x1.b81b456625f15p-3 <= c3 <= -0x1.b81b416e22329p-3
-0x1.e71d48d9c2ca4p-4 <= c5 <= -0x1.e71783472f5d1p-4
0x1.80e063cb210f9p-5 <= c7 <= 0x1.80ed6efa0a369p-5
0x1.1a3925ea0c5a9p-5 <= c9 <= 0x1.1b3783f148ed8p-5
-0x1.ec6032f293143p-7 <= c11 <= -0x1.e928025d508p-7
-0x1.8c06e851e2255p-7 <= c13 <= -0x1.732b2d4677028p-7
0x1.2aff33d629371p-8 <= c15 <= 0x1.41e9bc01ae472p-8
0x1.1e22f3192fd1dp-9 <= c17 <= 0x1.d851520a087c2p-9

Tomando c3 = -0x1.b81b44p-3, suponiendo que r8 también se evalúa en aritmética float :

-0x1.e71df05b5ad56p-4 <= c5 <= -0x1.e7175823ce2a4p-4
0x1.80df529dd8b18p-5 <= c7 <= 0x1.80f00e8da7f58p-5
0x1.1a283503e1a97p-5 <= c9 <= 0x1.1b5ca5beeeefep-5
-0x1.ed2c7cd87f889p-7 <= c11 <= -0x1.e8c17789776cdp-7
-0x1.90759e6defc62p-7 <= c13 <= -0x1.7045e66924732p-7
0x1.27eb51edf324p-8 <= c15 <= 0x1.47cda0bb1f365p-8
0x1.f6c6b51c50b54p-10 <= c17 <= 0x1.003a00ace9a79p-8

Tomando c5 = -0x1.e71aa4p-4, suponiendo que r7 se hace en aritmética float :

0x1.80e3dcc972cb3p-5 <= c7 <= 0x1.80ed1cf56977fp-5
0x1.1aa005ff6a6f4p-5 <= c9 <= 0x1.1afce9904742p-5
-0x1.ec7cf2464a893p-7 <= c11 <= -0x1.e9d6f7039db61p-7
-0x1.8a2304daefa26p-7 <= c13 <= -0x1.7a2456ddec8b2p-7
0x1.2e7b48f595544p-8 <= c15 <= 0x1.44437896b7049p-8
0x1.396f76c06de2ep-9 <= c17 <= 0x1.e3bedf4ed606dp-9

Tomando c7 = 0x1.80e87cp-5, suponiendo que r6 se hace en aritmética float :

0x1.1aa86d25bb64fp-5 <= c9 <= 0x1.1aca48cd5caabp-5
-0x1.eb6311f6c29dcp-7 <= c11 <= -0x1.eaedb032dfc0cp-7
-0x1.81438f115cbbp-7 <= c13 <= -0x1.7c9a106629f06p-7
0x1.36d433f81a012p-8 <= c15 <= 0x1.3babb57bb55bap-8
0x1.5cb14e1d4247dp-9 <= c17 <= 0x1.84f1151303aedp-9

Tomando c9 = 0x1.1ab95ap-5, suponiendo que r5 se hace en aritmética float :

-0x1.eb51a3b03781dp-7 <= c11 <= -0x1.eb21431536e0dp-7
-0x1.7fcd84700f7cfp-7 <= c13 <= -0x1.7ee38ee4beb65p-7
0x1.390fa00abaaabp-8 <= c15 <= 0x1.3b100a7f5d3cep-8
0x1.6ff147e1fdeb4p-9 <= c17 <= 0x1.7ebfed3ab5f9bp-9

Escogí un punto cerca del medio del rango para c11 y elegí aleatoriamente c13 , c15 y c17 .

EDITAR: ahora he automatizado este procedimiento. La siguiente función también es una implementación fielmente redondeada de arctan en [0, 1] :

float c5 = 0x1.997a72p-3;
float c7 = -0x1.23176cp-3;
float c9 = 0x1.b523c8p-4;
float c11 = -0x1.358ff8p-4;
float c13 = 0x1.61c5c2p-5;
float c15 = -0x1.0b16e2p-6;
float c17 = 0x1.7b422p-9;

float juffa_poly (float a) {
  float s = a * a;
  float r1 =              c17;
  float r2 = fmaf (r1, s, c15);
  float r3 = fmaf (r2, s, c13);
  float r4 = fmaf (r3, s, c11);
  float r5 = fmaf (r4, s, c9);
  float r6 = fmaf (r5, s, c7);
  float r7 = fmaf (r6, s, c5);
  float r8 = fmaf (r7, s, -0x1.5554dap-2f);
  float r9 = r8 * s;
  float r10 = fmaf (r9, a, a);
  return r10;
}

Me parece sorprendente que este código incluso exista. Para los coeficientes cercanos a estos, puedes obtener un límite en la distancia entre r10 y el valor del polinomio evaluado en aritmética real del orden de unos pocos ulps gracias a la lenta convergencia de este polinomio cuando s está cerca de 1 . Esperaba que el error de redondeo se comportara de una manera que era fundamentalmente "indomable" simplemente mediante el ajuste de coeficientes.

Answer 4

Reflexioné sobre las diversas ideas que recibí en los comentarios y también realicé algunos experimentos basados ​​en esa retroalimentación. Al final, decidí que una búsqueda heurística refinada era la mejor manera de avanzar. Ahora he logrado reducir el error máximo para atanf_poly() a 1.01036 ulps, con solo tres argumentos que exceden mi objetivo declarado de un límite de error de 1 ulp:

ulp = -1.00829 @ |a| =  9.80738342e-001 0x1.f62356p-1 (3f7b11ab)
ulp = -1.01036 @ |a| =  9.87551928e-001 0x1.f9a068p-1 (3f7cd034)
ulp =  1.00050 @ |a| =  9.99375939e-001 0x1.ffae34p-1 (3f7fd71a)

De acuerdo con la manera de generar la aproximación mejorada, no hay garantía de que sea la mejor aproximación; ningún avance científico aquí. Como el error ulp de la solución actual aún no está perfectamente equilibrado, y dado que continuar la búsqueda continúa entregando mejores aproximaciones (aunque a intervalos de tiempo exponencialmente crecientes), creo que se puede lograr un límite de error de 1 ulp, pero al mismo tiempo parecen estar muy cerca de la mejor aproximación optimizada para la máquina.

La mejor calidad de la nueva aproximación es el resultado de un proceso de búsqueda refinado. Observé que todos los errores ulp más grandes en el polinomio ocurren cerca de la unidad, por ejemplo en [0.75.1.0] para ser conservador. Esto permite un escaneo rápido para conjuntos de coeficientes interesantes cuyo error máximo es menor que algunos límites, digamos 1.08 ulps. Luego puedo probar en detalle y exhaustivamente todos los conjuntos de coeficientes dentro de un hipercono elegido heurísticamente y anclado en ese punto. Este segundo paso busca el error mínimo de ulp como el objetivo principal y el porcentaje máximo de resultados redondeados correctamente como un objetivo secundario. Al utilizar este proceso de dos pasos en los cuatro núcleos de mi CPU, pude acelerar significativamente el proceso de búsqueda: hasta ahora he podido verificar alrededor de 2 ²¹ conjuntos de coeficientes.

Basado en el rango de cada coeficiente en todas las soluciones "cercanas", ahora estimo que el espacio de búsqueda útil total para este problema de aproximación es> = 2 ²⁴ conjuntos de coeficientes en lugar del número más optimista de 2 ²⁰ que eché antes. Esto parece ser un problema factible de resolver para alguien que sea muy paciente o que tenga mucha potencia de cálculo.

Mi código actualizado es el siguiente:

// max ulp err = 1.01036
float atanf_poly (float a)
{
    float r, s;
    s = a * a;
    r =              0x1.7ed22cp-9f;
    r = fmaf (r, s, -0x1.0c2c2ep-6f);
    r = fmaf (r, s,  0x1.61fdf6p-5f);
    r = fmaf (r, s, -0x1.3556b4p-4f);
    r = fmaf (r, s,  0x1.b4e12ep-4f);
    r = fmaf (r, s, -0x1.230ae0p-3f);
    r = fmaf (r, s,  0x1.9978eep-3f);
    r = fmaf (r, s, -0x1.5554dap-2f);
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}

// max ulp err = 1.51871
float my_atanf (float a)
{
    float r, t;
    t = fabsf (a);
    r = t;
    if (t > 1.0f) {
        r = 1.0f / r;
    }
    r = atanf_poly (r);
    if (t > 1.0f) {
        r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r
    }
    r = copysignf (r, a);
    return r;
}

Actualización (después de revisar el problema dos años y medio después)

Usando el borrador de publicación de T. Myklebust como punto de partida, encontré la aproximación arcotangente en [-1,1] que tiene el error más pequeño para tener un error máximo de 0.94528 ulp.

/* Based on: Tor Myklebust, "Computing accurate Horner form approximations 
   to special functions in finite precision arithmetic", arXiv:1508.03211,
   August 2015. maximum ulp err = 0.94528
*/
float atanf_poly (float a)
{
    float r, s;
    s = a * a;                        
    r =              0x1.6d2086p-9f;  //  2.78569828e-3
    r = fmaf (r, s, -0x1.03f2ecp-6f); // -1.58660226e-2
    r = fmaf (r, s,  0x1.5beebap-5f); //  4.24722321e-2
    r = fmaf (r, s, -0x1.33194ep-4f); // -7.49753043e-2
    r = fmaf (r, s,  0x1.b403a8p-4f); //  1.06448799e-1
    r = fmaf (r, s, -0x1.22f5c2p-3f); // -1.42070308e-1
    r = fmaf (r, s,  0x1.997748p-3f); //  1.99934542e-1
    r = fmaf (r, s, -0x1.5554d8p-2f); // -3.33331466e-1
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}