punto - que es coma flotante en informatica

¿Hay algún documento que describa cómo Clang maneja el exceso de precisión de punto flotante?

Question

Es casi imposible (*) proporcionar semántica estricta de IEEE 754 a un costo razonable cuando las únicas instrucciones de punto flotante que uno puede usar son las 387. Es particularmente difícil cuando uno desea mantener la FPU trabajando en el significado completo de 64 bits, de modo que el tipo long double esté disponible para una precisión extendida. La "solución" usual es hacer cálculos intermedios con la única precisión disponible, y convertir a la precisión más baja en ocasiones más o menos bien definidas.

Las versiones recientes de GCC manejan el exceso de precisión en cálculos intermedios de acuerdo con la interpretación presentada por Joseph S. Myers en una publicación de 2008 en la lista de correo de GCC . Esta descripción hace que un programa compilado con gcc -std=c99 -mno-sse2 -mfpmath=387 completamente predecible, hasta el último bit, por lo que yo entiendo. Y si por casualidad no lo hace, es un error y se solucionará: la intención declarada de Joseph S. Myers en su publicación es hacer que sea predecible.

¿Está documentado cómo Clang maneja el exceso de precisión (digamos cuando se usa la opción -mno-sse2 ) y dónde?

(*) EDITAR: esto es una exageración. Es ligeramente molesto pero no tan difícil de emular binary64 cuando se permite configurar la FPU x87 para usar un significado de 53 bits.

Siguiendo un comentario de R .. a continuación, aquí está el registro de una breve interacción mía con la versión más reciente de Clang que tengo:

Hexa:~ $ clang -v
Apple clang version 4.1 (tags/Apple/clang-421.11.66) (based on LLVM 3.1svn)
Target: x86_64-apple-darwin12.4.0
Thread model: posix
Hexa:~ $ cat fem.c
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <fenv.h>

double x;
double y = 2.0;
double z = 1.0;

int main(){
  x = y + z;
  printf("%d/n", (int) FLT_EVAL_METHOD);
}
Hexa:~ $ clang -std=c99 -mno-sse2 fem.c
Hexa:~ $ ./a.out 
0
Hexa:~ $ clang -std=c99 -mno-sse2 -S fem.c
Hexa:~ $ cat fem.s 
…
    movl    $0, %esi
    fldl    _y(%rip)
    fldl    _z(%rip)
    faddp   %st(1)
    movq    _x@GOTPCREL(%rip), %rax
    fstpl   (%rax)
…

Answer 1

Esto no responde a la pregunta planteada originalmente, pero si usted es un programador que trabaja con problemas similares, esta respuesta podría ayudarlo.

Realmente no veo dónde está la dificultad percibida. Proporcionar semántica estricta IEEE-754 binary64 mientras se limita a 80387 matemática de coma flotante y mantener el cálculo doble largo de 80 bits, parece seguir reglas de conversión C99 bien especificadas con GCC-4.6.3 y clang-3.0 (basado en LLVM 3.0).

Editado para agregar: Sin embargo, Pascal Cuoq está en lo cierto: ni gcc-4.6.3 ni clang-llvm-3.0 realmente imponen esas reglas correctamente para ''387 matemática de coma flotante. Dadas las opciones de compilación adecuadas, las reglas se aplican correctamente a las expresiones evaluadas en tiempo de compilación, pero no para las expresiones de tiempo de ejecución. Hay soluciones alternativas, enumeradas después del descanso a continuación.

Hago un código de simulación de dinámica molecular, y estoy muy familiarizado con los requisitos de repetibilidad / predictibilidad y también con el deseo de mantener la máxima precisión disponible cuando sea posible, así que afirmo que sé de lo que estoy hablando aquí. Esta respuesta debería mostrar que las herramientas existen y son fáciles de usar; los problemas surgen de no ser consciente de o no usar esas herramientas.

(Un ejemplo preferido que me gusta, es el algoritmo de suma de Kahan. Con C99 y fundición adecuada (agregando moldes a, por ejemplo, el código de ejemplo de Wikipedia), no se necesitan trucos ni variables temporales adicionales. La implementación funciona independientemente del nivel de optimización del compilador, incluso a -Ofast y -Ofast .)

C99 declara explícitamente (en por ejemplo, 5.4.2.2) que tanto la fundición como la asignación eliminan todo el alcance y la precisión adicionales. Esto significa que puede usar la aritmética long double definiendo las variables temporales utilizadas durante el cálculo como el long double , también fundiendo sus variables de entrada a ese tipo; cada vez que se necesita un binary64 IEEE-754, simplemente envía a un double .

En ''387, el elenco genera una asignación y una carga en ambos compiladores anteriores; esto redondea correctamente el valor de 80 bits a IEEE-754 binary64. Este costo es muy razonable en mi opinión. El tiempo exacto que toma depende de la arquitectura y el código que lo rodea; por lo general es y puede ser intercalado con otro código para reducir el costo a niveles insignificantes. Cuando MMX, SSE o AVX están disponibles, sus registros están separados de los registros 80387 de 80 bits, y el reparto generalmente se realiza moviendo el valor al registro MMX / SSE / AVX.

(Prefiero el código de producción para usar un tipo de coma flotante específico, por ejemplo tempdouble o tempdouble , para variables temporales, de modo que se pueda definir como double o long double dependiendo de la arquitectura y de las compensaciones de velocidad / precisión deseadas).

En una palabra:

No asuma que (expression) es de double precisión solo porque todas las variables y constantes literales son. Escríbelo como (double)(expression) si quiere el resultado con double precisión.

Esto también se aplica a las expresiones compuestas, y algunas veces puede llevar a expresiones poco manejables con muchos niveles de moldes.

Si tiene expr1 y expr2 que desea calcular con una precisión de 80 bits, pero también necesita primero el producto de cada uno redondeado a 64 bits, use

long double  expr1;
long double  expr2;
double       product = (double)(expr1) * (double)(expr2);

Tenga en cuenta que el product se calcula como un producto de dos valores de 64 bits; no calculado con precisión de 80 bits, luego redondeado hacia abajo. El cálculo del producto con una precisión de 80 bits, y luego el redondeo hacia abajo, sería

double       other = expr1 * expr2;

o, agregando moldes descriptivos que le dicen exactamente lo que está sucediendo,

double       other = (double)((long double)(expr1) * (long double)(expr2));

Debería ser obvio que el product y other menudo difieren.

Las reglas de conversión C99 son solo otra herramienta que debes aprender a utilizar, si trabajas con valores mixtos de punto flotante de 32 bits / 64 bits / 80 bits / 128 bits. En realidad, te encuentras con los mismos problemas exactos si mezclas flotadores binarios32 y binarios64 ( float y double en la mayoría de las arquitecturas).

Quizás reescribir el código de exploración de Pascal Cuoq, para aplicar correctamente las reglas de fundición, lo aclara.

#include <stdio.h>

#define TEST(eq) printf("%-56s%s/n", "" # eq ":", (eq) ? "true" : "false")

int main(void)
{
    double d = 1.0 / 10.0;
    long double ld = 1.0L / 10.0L;

    printf("sizeof (double) = %d/n", (int)sizeof (double));
    printf("sizeof (long double) == %d/n", (int)sizeof (long double));

    printf("/nExpect true:/n");
    TEST(d == (double)(0.1));
    TEST(ld == (long double)(0.1L));
    TEST(d == (double)(1.0 / 10.0));
    TEST(ld == (long double)(1.0L / 10.0L));
    TEST(d == (double)(ld));
    TEST((double)(1.0L/10.0L) == (double)(0.1));
    TEST((long double)(1.0L/10.0L) == (long double)(0.1L));

    printf("/nExpect false:/n");
    TEST(d == ld);
    TEST((long double)(d) == ld);
    TEST(d == 0.1L);
    TEST(ld == 0.1);
    TEST(d == (long double)(1.0L / 10.0L));
    TEST(ld == (double)(1.0L / 10.0));

    return 0;
}

La salida, tanto con GCC como con clang, es

sizeof (double) = 8
sizeof (long double) == 12

Expect true:
d == (double)(0.1):                                     true
ld == (long double)(0.1L):                              true
d == (double)(1.0 / 10.0):                              true
ld == (long double)(1.0L / 10.0L):                      true
d == (double)(ld):                                      true
(double)(1.0L/10.0L) == (double)(0.1):                  true
(long double)(1.0L/10.0L) == (long double)(0.1L):       true

Expect false:
d == ld:                                                false
(long double)(d) == ld:                                 false
d == 0.1L:                                              false
ld == 0.1:                                              false
d == (long double)(1.0L / 10.0L):                       false
ld == (double)(1.0L / 10.0):                            false

excepto que las versiones recientes de GCC promueven el lado derecho de ld == 0.1 a largo doble primero (es decir, a ld == 0.1L ), dando true , y que con SSE / AVX, el long double es 128-bit.

Para las pruebas ''387 puras, utilicé

gcc -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test
clang -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test

con varias combinaciones de indicadores de optimización como ... , incluidos -fomit-frame-pointer , -O0 , -O1 , -O2 , -Os y -Os .

El uso de cualquier otro marcador o compilador C99 debería conducir a los mismos resultados, excepto para el long double tamaño long double (y ld == 1.0 para las versiones actuales de GCC). Si encuentra alguna diferencia, estaría muy agradecido de escuchar sobre ellos; Puede que necesite advertir a mis usuarios de dichos compiladores (versiones de compilador). Tenga en cuenta que Microsoft no es compatible con C99, por lo que no me interesan por completo.

Pascal Cuoq plantea un problema interesante en la cadena de comentarios a continuación, que no reconocí de inmediato.

Al evaluar una expresión, tanto GCC como clang con -mfpmath=387 especifican que todas las expresiones se evalúan utilizando una precisión de 80 bits. Esto lleva a, por ejemplo,

7491907632491941888 = 0x1.9fe2693112e14p+62 = 110011111111000100110100100110001000100101110000101000000000000
5698883734965350400 = 0x1.3c5a02407b71cp+62 = 100111100010110100000001001000000011110110111000111000000000000

7491907632491941888 * 5698883734965350400 = 42695510550671093541385598890357555200 = 100000000111101101101100110001101000010100100001011110111111111111110011000111000001011101010101100011000000000000000000000000

arrojando resultados incorrectos, porque esa cadena de unos en el medio del resultado binario está justo en la diferencia entre las mantisas de 53 y 64 bits (números de coma flotante de 64 y 80 bits, respectivamente). Entonces, mientras el resultado esperado es

42695510550671088819251326462451515392 = 0x1.00f6d98d0a42fp+125 = 100000000111101101101100110001101000010100100001011110000000000000000000000000000000000000000000000000000000000000000000000000

el resultado obtenido con just -std=c99 -m32 -mno-sse -mfpmath=387 es

42695510550671098263984292201741942784 = 0x1.00f6d98d0a43p+125 = 100000000111101101101100110001101000010100100001100000000000000000000000000000000000000000000000000000000000000000000000000000

En teoría, debería poder decirle a gcc y clang que apliquen las reglas correctas de redondeo C99 utilizando las opciones

-std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard

Sin embargo, esto solo afecta a las expresiones que el compilador optimiza, y no parece arreglar el manejo del 387 en absoluto. Si utiliza, por ejemplo, clang -O1 -std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard test.c -o test && ./test con test.c siendo el programa de ejemplo de Pascal Cuoq , obtendrá el resultado correcto según las reglas IEEE-754, pero solo porque el compilador optimiza la expresión, sin utilizar el 387 en absoluto.

En pocas palabras, en lugar de informática

(double)d1 * (double)d2

tanto gcc como clang en realidad le dicen al ''387 que calcule

(double)((long double)d1 * (long double)d2)

Esto es de hecho Creo que este es un error del compilador que afecta tanto a gcc-4.6.3 como a clang-llvm-3.0, y es fácilmente reproducible. (Pascal Cuoq señala que FLT_EVAL_METHOD=2 significa que las operaciones con argumentos de doble precisión siempre se realizan con mayor precisión, pero no veo ningún motivo razonable, aparte de tener que reescribir partes de la libm en ''387, para hacer eso en C99 y considerando que las reglas IEEE-754 son alcanzables por el hardware. Después de todo, el compilador puede realizar fácilmente la operación correcta , modificando la palabra de control ''387 para que coincida con la precisión de la expresión. Y, dadas las opciones del compilador que deberían forzar esto comportamiento - -std=c99 -ffloat-store -fexcess-precision=standard - no tiene sentido si realmente se desea el comportamiento FLT_EVAL_METHOD=2 , tampoco hay problemas de compatibilidad con versiones anteriores.) Es importante tener en cuenta que dado el correcto indicadores del compilador, las expresiones evaluadas en tiempo de compilación se evalúan correctamente y solo las expresiones evaluadas en tiempo de ejecución obtienen resultados incorrectos.

La solución más simple, y la portátil, es usar fesetround(FE_TOWARDZERO) (de fenv.h ) para redondear todos los resultados hacia cero.

En algunos casos, redondear hacia cero puede ayudar con la predictibilidad y los casos patológicos. En particular, para intervalos como x = [0,1) , el redondeo hacia cero significa que el límite superior nunca se alcanza mediante el redondeo; importante si evalúa, por ejemplo, splines por partes.

Para los otros modos de redondeo, debe controlar el hardware 387 directamente.

Puede usar __FPU_SETCW() desde #include <fpu_control.h> , o abrirlo por código. Por ejemplo, precision.c :

#include <stdlib.h>
#include <stdio.h>
#include <limits.h>

#define FP387_NEAREST   0x0000
#define FP387_ZERO      0x0C00
#define FP387_UP        0x0800
#define FP387_DOWN      0x0400

#define FP387_SINGLE    0x0000
#define FP387_DOUBLE    0x0200
#define FP387_EXTENDED  0x0300

static inline void fp387(const unsigned short control)
{
    unsigned short cw = (control & 0x0F00) | 0x007f;
    __asm__ volatile ("fldcw %0" : : "m" (*&cw));
}

const char *bits(const double value)
{
    const unsigned char *const data = (const unsigned char *)&value;
    static char buffer[CHAR_BIT * sizeof value + 1];
    char       *p = buffer;
    size_t      i = CHAR_BIT * sizeof value;
    while (i-->0)
        *(p++) = ''0'' + !!(data[i / CHAR_BIT] & (1U << (i % CHAR_BIT)));
    *p = ''/0'';
    return (const char *)buffer;
}


int main(int argc, char *argv[])
{
    double  d1, d2;
    char    dummy;

    if (argc != 3) {
        fprintf(stderr, "/nUsage: %s 7491907632491941888 5698883734965350400/n/n", argv[0]);
        return EXIT_FAILURE;
    }

    if (sscanf(argv[1], " %lf %c", &d1, &dummy) != 1) {
        fprintf(stderr, "%s: Not a number./n", argv[1]);
        return EXIT_FAILURE;
    }
    if (sscanf(argv[2], " %lf %c", &d2, &dummy) != 1) {
        fprintf(stderr, "%s: Not a number./n", argv[2]);
        return EXIT_FAILURE;
    }

    printf("%s:/td1 = %.0f/n/t    %s in binary/n", argv[1], d1, bits(d1));
    printf("%s:/td2 = %.0f/n/t    %s in binary/n", argv[2], d2, bits(d2));

    printf("/nDefaults:/n");
    printf("Product = %.0f/n/t  %s in binary/n", d1 * d2, bits(d1 * d2));

    printf("/nExtended precision, rounding to nearest integer:/n");
    fp387(FP387_EXTENDED | FP387_NEAREST);
    printf("Product = %.0f/n/t  %s in binary/n", d1 * d2, bits(d1 * d2));

    printf("/nDouble precision, rounding to nearest integer:/n");
    fp387(FP387_DOUBLE | FP387_NEAREST);
    printf("Product = %.0f/n/t  %s in binary/n", d1 * d2, bits(d1 * d2));

    printf("/nExtended precision, rounding to zero:/n");
    fp387(FP387_EXTENDED | FP387_ZERO);
    printf("Product = %.0f/n/t  %s in binary/n", d1 * d2, bits(d1 * d2));

    printf("/nDouble precision, rounding to zero:/n");
    fp387(FP387_DOUBLE | FP387_ZERO);
    printf("Product = %.0f/n/t  %s in binary/n", d1 * d2, bits(d1 * d2));

    return 0;
}

Usando clang-llvm-3.0 para compilar y ejecutar, obtengo los resultados correctos,

clang -std=c99 -m32 -mno-sse -mfpmath=387 -O3 -W -Wall precision.c -o precision
./precision 7491907632491941888 5698883734965350400

7491907632491941888:    d1 = 7491907632491941888
        0100001111011001111111100010011010010011000100010010111000010100 in binary
5698883734965350400:    d2 = 5698883734965350400
        0100001111010011110001011010000000100100000001111011011100011100 in binary

Defaults:
Product = 42695510550671098263984292201741942784
          0100011111000000000011110110110110011000110100001010010000110000 in binary

Extended precision, rounding to nearest integer:
Product = 42695510550671098263984292201741942784
          0100011111000000000011110110110110011000110100001010010000110000 in binary

Double precision, rounding to nearest integer:
Product = 42695510550671088819251326462451515392
          0100011111000000000011110110110110011000110100001010010000101111 in binary

Extended precision, rounding to zero:
Product = 42695510550671088819251326462451515392
          0100011111000000000011110110110110011000110100001010010000101111 in binary

Double precision, rounding to zero:
Product = 42695510550671088819251326462451515392
          0100011111000000000011110110110110011000110100001010010000101111 in binary

En otras palabras, puede solucionar los problemas del compilador utilizando fp387() para establecer la precisión y el modo de redondeo.

La desventaja es que algunas librerías matemáticas ( libm.a , libm.so ) pueden escribirse con la suposición de que los resultados intermedios siempre se computan con una precisión de 80 bits. Al menos la biblioteca C de GNU fpu_control.h en x86_64 tiene el comentario "la libm requiere una precisión extendida" . Afortunadamente, puede tomar las implementaciones ''387 de, por ejemplo, la biblioteca GNU C, e implementarlas en un archivo de cabecera o escribir una libm trabajo libm , si necesita la funcionalidad math.h ; de hecho, creo que podría ayudar allí.

Answer 2

Para el registro, a continuación se encuentra lo que encontré por experimentación. El siguiente programa muestra varios comportamientos cuando se compila con Clang:

#include <stdio.h>

int r1, r2, r3, r4, r5, r6, r7;

double ten = 10.0;

int main(int c, char **v)
{
  r1 = 0.1 == (1.0 / ten);
  r2 = 0.1 == (1.0 / 10.0);
  r3 = 0.1 == (double) (1.0 / ten);
  r4 = 0.1 == (double) (1.0 / 10.0);
  ten = 10.0;
  r5 = 0.1 == (1.0 / ten);
  r6 = 0.1 == (double) (1.0 / ten);
  r7 = ((double) 0.1) == (1.0 / 10.0);
  printf("r1=%d r2=%d r3=%d r4=%d r5=%d r6=%d r7=%d/n", r1, r2, r3, r4, r5, r6, r7);
}

Los resultados varían con el nivel de optimización:

$ clang -v
Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn)
$ clang -mno-sse2 -std=c99  t.c && ./a.out
r1=0 r2=1 r3=0 r4=1 r5=1 r6=0 r7=1
$ clang -mno-sse2 -std=c99 -O2  t.c && ./a.out
r1=0 r2=1 r3=0 r4=1 r5=1 r6=1 r7=1

El elenco (double) que diferencia r5 y r6 a -O2 no tiene efecto en -O0 y para las variables r3 y r4 . El resultado r1 es diferente de r5 en todos los niveles de optimización, mientras que r6 solo difiere de r3 en -O2 .