Life hack: cómo analizar un gigabyte de dobles por segundo

¿Cómo leer un valor doble de una cadena en código C ++?

std::stringstream in(mystring);
while(in >> x) {
   sum += x;
}
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

En Intel Skylake con el compilador GCC 8.3, este código analiza 50 MB / s. Los discos duros pueden proporcionar fácilmente lecturas secuenciales a una velocidad de varios GB / s, por lo que no hay duda de que no es la velocidad de lectura del disco lo que nos limita, sino la velocidad de análisis. ¿Cómo puedo acelerarlo?



Lo primero que se sugiere a sí mismo es abandonar la conveniencia proporcionada por los flujos en C ++ y llamar a strtod (3) directamente:

do {
    number = strtod(s, &end);
    if(end == s) break;
    sum += number;
    s = end; 
} while (s < theend);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

La velocidad aumenta hasta 90 MB / s; el perfil muestra que cuando se lee de la secuencia se ejecutan ~ 1600 instrucciones para cada número leído, cuando se usa strtod - ~ 1100 instrucciones por número. Las bibliotecas estándar C y C ++ pueden justificarse por los requisitos de universalidad y portabilidad; pero si nos limitamos a analizar solo el doble y solo en x64, entonces podemos escribir un código mucho más eficiente: 280 instrucciones por número son suficientes.

Analizar enteros

Comencemos con una subtarea: dado un número decimal de ocho dígitos, necesitamos convertirlo a int. En la escuela, a todos nos enseñaron a hacer esto en un ciclo:

int sum = 0;
for (int i = 0; i <= 7; i++)
{
	sum = (sum * 10) + (num[i] - '0');
}
return sum;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Compilado GCC 9.3.1 -O3, para mí este código maneja 3 GB / s. La forma obvia de acelerarlo es invertir el ciclo; pero el compilador lo hace por sí mismo.



Tenga en cuenta que la notación decimal de un número de ocho dígitos encaja en la variable int64_t: por ejemplo, la cadena "87654321" se almacena de la misma manera que el valor int64_t 0x3132333435363738 (el primer byte contiene los 8 bits menos significativos, el último - los más significativos). Esto significa que en lugar de ocho accesos a la memoria, podemos asignar el valor de cada dígito con un desplazamiento:

int64_t sum = *(int64_t*)num;
return (sum      & 15) * 10000000 +
    ((sum >> 8)  & 15) * 1000000 +
    ((sum >> 16) & 15) * 100000 +
    ((sum >> 24) & 15) * 10000 +
    ((sum >> 32) & 15) * 1000 +
    ((sum >> 40) & 15) * 100 +
    ((sum >> 48) & 15) * 10 +
    ((sum >> 56) & 15);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Todavía no hay aceleración; por el contrario, ¡la velocidad desciende a 1,7 GB / s! Vayamos más allá: Y con una máscara de 0x000F000F000F000F nos dará 0x0002000400060008 - dígitos decimales en posiciones pares. Combinemos cada uno de ellos con lo siguiente:

sum = ((sum & 0x000F000F000F000F) * 10) + 
      ((sum & 0x0F000F000F000F00) >> 8);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Después de eso, 0x3132333435363738 se convierte en 0x00 (12) 00 (34) 00 (56) 00 (78) - los bytes en las posiciones pares se ponen a cero, en los impares - contienen representaciones binarias de números decimales de dos dígitos.

return (sum        & 255) * 1000000 +
      ((sum >> 16) & 255) * 10000 +
      ((sum >> 32) & 255) * 100 +
      ((sum >> 48) & 255);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

- completa la conversión de un número de ocho dígitos.



El mismo truco se puede repetir con pares de números de dos dígitos:

sum = ((sum & 0x0000007F0000007F) * 100) +
      ((sum >> 16) & 0x0000007F0000007F);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

- 0x00 (12) 00 (34) 00 (56) 00 (78) se convierte en 0x0000 (1234) 0000 (5678);

y con el par resultante de cuatro dígitos:

return ((sum & 0x3FFF) * 10000) + ((sum >> 32) & 0x3FFF);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

- 0x00 (12) 00 (34) 00 (56) 00 (78) se convierte en 0x00000000 (12345678). La mitad menor del int64_t resultante es el resultado deseado. A pesar de la notable simplificación del código (tres multiplicaciones en lugar de siete), la velocidad de análisis (2.6 GB / s) sigue siendo peor que la de la implementación ingenua.



Pero la combinación de pares de números se puede simplificar incluso si observa que la siguiente acción aplicará la máscara 0x007F007F007F007F, lo que significa que cualquier basura se puede dejar en los bytes que aceptan valores NULL:

sum = ((sum & 0x0?0F0?0F0?0F0?0F) * 10) + ((sum & 0x0F??0F??0F??0F??) >> 8) =
   = (((sum & 0x0F0F0F0F0F0F0F0F) * 2560)+ (sum & 0x0F0F0F0F0F0F0F0F)) >> 8 =
    = ((sum & 0x0F0F0F0F0F0F0F0F) * 2561) >> 8;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

En términos simplificados, una máscara en lugar de dos, y no hay adición. Las dos expresiones restantes se simplifican de la misma manera:

sum = ((sum & 0x00FF00FF00FF00FF) * 6553601) >> 16;
return((sum & 0x0000FFFF0000FFFF) * 42949672960001) >> 32;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Este truco se llama SIMD dentro de un registro (SWAR): usándolo, la velocidad de análisis aumenta a 3.6 GB / s.



Se puede usar un truco SWAR similar para verificar si una cadena de ocho caracteres es un número decimal de ocho dígitos:

return ((val & 0xF0F0F0F0F0F0F0F0) |
      (((val + 0x0606060606060606) & 0xF0F0F0F0F0F0F0F0) >> 4))
            == 0x3333333333333333;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Dispositivo doble

El estándar IEEE asignó 52 bits a la mantisa y 11 al exponente dentro de estos números; esto significa que el número se almacena como $$$\pm1.m\cdot2^e$ dónde $$$0\le m<2^{52}<10^{16}$ y $$$-1022\le e\le+1023$ ... En particular, esto significa que en la notación decimal de doble, solo los 17 dígitos más significativos son significativos; los bits menos significativos no pueden afectar el valor doble de ninguna manera. Incluso 17 dígitos significativos es mucho más de lo que podría ser necesario para cualquier aplicación práctica: los números desnormalizados complican un poco el trabajo (desde







$$$2^{-1074}$ antes de $$$2^{-1022}$ con un número correspondientemente menor de bits significativos en la mantisa) y requisitos de redondeo (cualquier número decimal debe estar representado por el binario más cercano, y si el número está exactamente en el medio entre el binario más cercano, entonces se prefiere la mantisa par). Todo esto asegura que si la computadora A convierte el número X en una cadena decimal con 17 dígitos significativos, entonces cualquier computadora B, que lea esta cadena, recibirá el mismo número X, bit por bit, independientemente de si A y B son iguales. modelos de procesador, sistemas operativos y lenguajes de programación. (Imagínese lo divertido que sería depurar su código si los errores de redondeo fueran diferentes en diferentes computadoras). Debido a los requisitos de redondeo, surgen los malentendidos mencionados recientemente . en Habré: la fracción decimal 0.1 se representa como la fracción binaria más cercana $$$7205759403792794\cdot2^{-56}$ igual a exactamente 0.1000000000000000055511151231257827021181583404541015625; 0.2 - en la forma $$$7205759403792794\cdot2^{-55}$ igual a exactamente 0,200000000000000011102230246251565404236316680908203125; pero su suma no es una fracción binaria más cercana a 0.3: aproximación desde abajo $$$5404319552844595\cdot2^{-54}$ = 0.299999999999999988897769753748434595763683319091796875 resulta ser más preciso. Por lo tanto, el estándar IEEE requiere agregar 0.1 + 0.2 para producir un resultado diferente a 0.3.

Analizando doble

Una generalización sencilla del algoritmo de enteros consiste en multiplicaciones iterativas y divisiones por 10.0; a diferencia del análisis de enteros, aquí es necesario procesar dígitos de menor a mayor para que los errores de redondeo en los dígitos altos "oculten" los errores de redondeo en los bajos. Al mismo tiempo, el análisis sintáctico de la mantisa se reduce fácilmente al análisis sintáctico de números enteros: basta con cambiar la normalización para que el punto binario imaginario no esté al principio de la mantisa, sino al final; multiplicar o dividir el número entero de 53 bits resultante por la potencia deseada de diez; y finalmente, reste 52 del exponente, es decir mueva el punto 52 bits a la izquierda, donde debería estar de acuerdo con el estándar IEEE. Además, hay tres hechos importantes a tener en cuenta:

1. Basta con aprender a multiplicar o dividir por 5, y la multiplicación o división por 2 es solo un incremento o decremento de un exponente;
2. uint64_t 5 0xcccccccccccccccd 66 , , $$$\frac{\texttt{0xcccccccccccccccd}}{2^{66}} - \frac1 5 = \frac1{5\cdot2^{66}} < 2^{-68}$ 64 ( );
3. – $$$10^{-324} < 2^{-1074}$ $$$2^{1024} < 10^{309}$ ; , 309 , 324 0xcccccccccccccccd, . ( 53 ; 128- , 53- 53- .) 633 double ( , ⅕, 7205759403792794 – 0xcccccccccccccccd, 53 ), double ; 53 , . , , 64 64 , , 128- . .

La complejidad del redondeo estándar es que para saber que el resultado está exactamente en el medio entre las fracciones binarias más cercanas, no solo necesitamos 54 bits significativos del resultado (el bit cero menos significativo significa "todo está en orden", el uno significa "golpear el medio"), pero debe observar si hubo una continuación distinta de cero después del bit menos significativo. En particular, esto significa que considerar solo los 17 dígitos más significativos en la notación decimal del número no es suficiente: definen la mantisa binaria solo con una precisión de ± 1, y para seleccionar la dirección de redondeo deseada, tendrá que considerar los dígitos inferiores. Por ejemplo, 10000000000000003 y 10000000000000005 son el mismo valor doble (igual a 10000000000000004) y 10000000000000005.00 (muchos ceros) 001 es un valor diferente (igual a 10000000000000006).



El profesor Daniel Lemire de la Universidad por correspondencia de Quebec (TÉLUQ), que inventó este analizador, lo publicó en github . Esta biblioteca proporciona dos funciones from_chars



: una analiza una cadena para duplicarla y la otra para flotar. Sus experimentos demostraron que en el 99,8% de los casos basta con considerar 19 dígitos decimales significativos (64 bits): si para dos valores consecutivos de una mantisa de 64 bits se obtiene el mismo valor doble final, entonces este es el valor correcto. . Solo en el 0,2% de los casos, cuando estos dos valores no coinciden, se ejecuta un código más complejo y más universal que siempre implementa el redondeo correcto.