Hack de vida: como analisar um gigabyte de duplos por segundo

Como ler um valor duplo de uma string no código C ++?

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

No Intel Skylake com compilador GCC 8.3, este código analisa 50 MB / s. Os discos rígidos podem fornecer facilmente leituras sequenciais a uma velocidade de vários GB / s, portanto, não há dúvida de que não é a velocidade de leitura do disco que nos limita, mas a velocidade de análise. Como posso acelerar isso?



A primeira coisa que se sugere é abandonar a conveniência fornecida por streams em C ++ e chamar strtod (3) diretamente:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

A velocidade cresce até 90 MB / s; o perfil mostra que, ao ler o fluxo, ~ 1600 instruções são executadas para cada número lido, ao usar strtod - ~ 1100 instruções por número. As bibliotecas padrão C e C ++ podem ser justificadas pelos requisitos de universalidade e portabilidade; mas se nos limitarmos a analisar apenas duas vezes e apenas em x64, poderemos escrever um código muito mais eficiente: 280 instruções por número são suficientes.

Analisando inteiros

Vamos começar com uma subtarefa: dado um número decimal de oito dígitos, precisamos convertê-lo para int. Na escola, todos fomos ensinados a fazer isso em um 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 mim este código lida com 3 GB / s. A maneira óbvia de acelerá-lo é inverter o ciclo; mas o compilador faz isso sozinho.



Observe que a notação decimal de um número de oito dígitos se encaixa na variável int64_t: por exemplo, a string "87654321" é armazenada da mesma maneira que o valor int64_t 0x3132333435363738 (o primeiro byte contém os 8 bits menos significativos, o último - os mais significativos). Isso significa que em vez de oito acessos à memória, podemos alocar o valor de cada dígito com um deslocamento:

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);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Ainda não há aceleração; pelo contrário, a velocidade cai para 1,7 GB / s! Vamos além: E com a máscara 0x000F000F000F000F nos dará 0x0002000400060008 - dígitos decimais em posições pares. Vamos combinar cada um deles com o seguinte:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Depois disso, 0x3132333435363738 se transforma em 0x00 (12) 00 (34) 00 (56) 00 (78) - bytes nas posições pares são zerados, nas ímpares - eles contêm representações binárias de números decimais de dois dígitos.

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

- completa a conversão de um número de oito dígitos.



O mesmo truque pode ser repetido com pares de números de dois dígitos:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

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

e com o par resultante de quatro dígitos:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

- 0x00 (12) 00 (34) 00 (56) 00 (78) torna-se 0x00000000 (12345678). A menor metade do int64_t resultante é o resultado desejado. Apesar da notável simplificação do código (três multiplicações em vez de sete), a velocidade de análise (2,6 GB / s) permanece pior do que a implementação ingênua.



Mas a combinação de pares de números pode ser simplificada, mesmo se você perceber que a seguinte ação aplicará a máscara 0x007F007F007F007F, o que significa que qualquer lixo pode ser deixado nos bytes anuláveis:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Em termos simplificados, uma máscara em vez de duas e não há adição. As duas expressões restantes são simplificadas da mesma maneira:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Esse truque é chamado de SIMD dentro de um registrador (SWAR): usando-o, a velocidade de análise cresce para 3,6 GB / s.



Um truque SWAR semelhante pode ser usado para verificar se uma string de oito caracteres é um número decimal de oito dígitos:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Dispositivo duplo

O padrão IEEE atribuiu 52 bits à mantissa e 11 ao expoente dentro desses números; isso significa que o número é armazenado como $$$\pm1.m\cdot2^e$ Onde $$$0\le m<2^{52}<10^{16}$ e $$$-1022\le e\le+1023$ ... Em particular, isso significa que na notação decimal dupla, apenas os 17 dígitos mais significativos são significativos; os bits menos significativos não podem afetar o valor duplo de forma alguma. Mesmo 17 dígitos significativos é muito mais do que pode ser necessário para qualquer aplicação prática: Números desnormalizados complicam um pouco o trabalho (de







$$$2^{-1074}$ antes $$$2^{-1022}$ com um número correspondentemente menor de bits significativos na mantissa) e requisitos de arredondamento (qualquer número decimal deve ser representado pelo binário mais próximo, e se o número estiver exatamente no meio entre o binário mais próximo, então a mantissa par é preferida). Tudo isso garante que se o computador A converter o número X em uma string decimal com 17 dígitos significativos, qualquer computador B, lendo essa string, receberá o mesmo número X, bit por bit - independentemente de A e B serem iguais. modelos de processador, sistemas operacionais e linguagens de programação. (Imagine como seria divertido depurar seu código se os erros de arredondamento fossem diferentes em computadores diferentes.) Devido aos requisitos de arredondamento, surgem os mal-entendidos mencionados recentemente . em Habré: a fração decimal 0,1 é representada como a fração binária mais próxima a ela $$$7205759403792794\cdot2^{-56}$ igual a exatamente 0,1000000000000000055511151231257827021181583404541015625; 0,2 - na forma $$$7205759403792794\cdot2^{-55}$ igual a exatamente 0,200000000000000011102230246251565404236316680908203125; mas sua soma não é uma fração binária mais próxima de 0,3: aproximação de baixo $$$5404319552844595\cdot2^{-54}$ = 0,299999999999999988897769753748434595763683319091796875 acaba sendo mais preciso. Portanto, o padrão IEEE requer a adição de 0,1 + 0,2 para produzir um resultado diferente de 0,3.

Análise dupla

Uma generalização direta do algoritmo de inteiros consiste em multiplicações e divisões iterativas por 10,0 - ao contrário da análise de inteiros, aqui é necessário processar dígitos de baixo para cima para que erros de arredondamento em dígitos altos "escondam" erros de arredondamento em dígitos baixos. Ao mesmo tempo, a análise da mantissa é facilmente reduzida à análise de inteiros: basta alterar a normalização para que o ponto binário imaginário não esteja no início da mantissa, mas no final; o inteiro de 53 bits resultante se multiplica ou divide pela potência desejada de dez; e, finalmente, subtraia 52 do expoente, ou seja, mova o ponto 52 bits para a esquerda - onde deveria estar de acordo com o padrão IEEE. Além disso, existem três fatos importantes a serem observados:

1. Basta aprender a multiplicar ou dividir por 5, e a multiplicação ou divisão por 2 é apenas um incremento ou decremento de um expoente;
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- . .

A complexidade do arredondamento padrão é que, para descobrir que o resultado está exatamente no meio entre as frações binárias mais próximas, não precisamos apenas de 54 bits significativos do resultado (o bit zero menos significativo significa "está tudo em ordem", o um significa "atingiu o meio"), mas e você precisa observar se houve uma continuação diferente de zero após o bit menos significativo. Em particular, isso significa que considerar apenas os 17 dígitos mais significativos na notação decimal do número não é suficiente: eles definem a mantissa binária apenas com uma precisão de ± 1, e para selecionar a direção de arredondamento desejada, você terá que considerar os dígitos mais baixos. Por exemplo, 10000000000000003 e 10000000000000005 têm o mesmo valor duplo (igual a 10000000000000004) e 10000000000000005,00 (muitos zeros) 001 é um valor diferente (igual a 10000000000000006).



O professor Daniel Lemire, da Correspondence University of Quebec (TÉLUQ), que inventou este parser, publicou-o no github . Esta biblioteca oferece duas funções from_chars



: uma analisa uma string para duplicar e a outra para flutuar. Seus experimentos mostraram que em 99,8% dos casos é suficiente considerar 19 dígitos decimais significativos (64 bits): se para dois valores consecutivos de uma mantissa de 64 bits o mesmo valor duplo final é obtido, então este é o valor correto . Somente em 0,2% dos casos, quando esses dois valores não coincidem, é executado um código mais complexo e universal que implementa o arredondamento sempre correto.