Life Hack: comment analyser un gigaoctet de doubles par seconde

Comment lire une valeur double à partir d'une chaîne en code C ++?

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Sur Intel Skylake avec le compilateur GCC 8.3, ce code analyse 50 Mo / s. Les disques durs peuvent facilement fournir des lectures séquentielles à une vitesse de plusieurs Go / s, il ne fait donc aucun doute que ce n'est pas la vitesse de lecture à partir du disque qui nous limite, mais la vitesse d'analyse. Comment puis-je l'accélérer?



La première chose qui se suggère est d'abandonner la commodité fournie par les flux en C ++ et d'appeler directement strtod (3):

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

La vitesse augmente jusqu'à 90 Mo / s; le profilage montre que lors de la lecture à partir du flux, ~ 1600 instructions sont exécutées pour chaque numéro lu, lors de l'utilisation de strtod - ~ 1100 instructions par numéro. Les bibliothèques standard C et C ++ peuvent être justifiées par les exigences d'universalité et de portabilité; mais si nous nous limitons à l'analyse uniquement en double et uniquement sur x64, alors nous pouvons écrire du code beaucoup plus efficace: 280 instructions par nombre suffisent.

Analyse des entiers

Commençons par une sous-tâche: étant donné un nombre décimal à huit chiffres, nous devons le convertir en int. À l'école, on nous a tous appris à faire cela dans un cycle:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Compilé GCC 9.3.1 -O3, pour moi ce code gère 3 Go / s. La manière évidente de l'accélérer est d'inverser la boucle; mais le compilateur le fait lui-même.



Notez que la notation décimale d'un nombre à huit chiffres s'inscrit dans la variable int64_t: par exemple, la chaîne «87654321» est stockée de la même manière que la valeur int64_t 0x3132333435363738 (le premier octet contient les 8 bits les moins significatifs, le dernier - les plus significatifs). Cela signifie qu'au lieu de huit accès mémoire, nous pouvons allouer la valeur de chaque chiffre avec un décalage:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Il n'y a pas encore d'accélération; au contraire, la vitesse tombe à 1,7 Go / s! Allons plus loin: ET avec le masque 0x000F000F000F000F nous donnera 0x0002000400060008 - chiffres décimaux dans des positions paires. Combinons chacun d'eux avec ce qui suit:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Après cela, 0x3132333435363738 se transforme en 0x00 (12) 00 (34) 00 (56) 00 (78) - les octets aux positions paires sont mis à zéro, aux impairs - ils contiennent des représentations binaires de nombres décimaux à deux chiffres.

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

- termine la conversion d'un nombre à huit chiffres.



La même astuce peut être répétée avec des paires de nombres à deux chiffres:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

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

et avec la paire résultante de ceux à quatre chiffres:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

- 0x00 (12) 00 (34) 00 (56) 00 (78) devient 0x00000000 (12345678). La moitié inférieure de l'int64_t résultant est le résultat souhaité. Malgré la simplification notable du code (trois multiplications au lieu de sept), la vitesse d'analyse (2,6 Go / s) reste pire que celle de l'implémentation naïve.



Mais la combinaison de paires de nombres peut être simplifiée même si vous remarquez que l'action suivante appliquera le masque 0x007F007F007F007F, ce qui signifie que tout garbage peut être laissé dans les octets Nullable:

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 termes simplifiés, un masque au lieu de deux, et il n'y a pas d'ajout. Les deux expressions restantes sont simplifiées de la même manière:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Cette astuce s'appelle SIMD dans un registre (SWAR): en l'utilisant, la vitesse d'analyse passe à 3,6 Go / s.



Une astuce SWAR similaire peut être utilisée pour vérifier si une chaîne de huit caractères est un nombre décimal à huit chiffres:

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

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
    

Double appareil

La norme IEEE attribuait 52 bits à la mantisse et 11 à l'exposant dans ces nombres; cela signifie que le numéro est stocké sous $$$\pm1.m\cdot2^e$ , où $$$0\le m<2^{52}<10^{16}$ et $$$-1022\le e\le+1023$ . En particulier, cela signifie que dans la notation décimale de double, seuls les 17 chiffres les plus significatifs sont significatifs; les bits les moins significatifs ne peuvent en aucun cas affecter la valeur double. Même 17 chiffres significatifs sont bien plus que ce qui pourrait être nécessaire pour toute application pratique: les nombres dénormalisés compliquent un peu le travail (de







$$$2^{-1074}$ à $$$2^{-1022}$ avec un nombre proportionnellement plus petit de bits significatifs dans la mantisse) et des exigences d'arrondi (tout nombre décimal doit être représenté par le binaire le plus proche, et si le nombre est exactement au milieu entre le binaire le plus proche, alors la mantisse paire est préférée ). Tout cela garantit que si l'ordinateur A convertit le nombre X en une chaîne décimale avec 17 chiffres significatifs, alors tout ordinateur B, lisant cette chaîne, recevra le même nombre X, bit pour bit - indépendamment du fait que A et B soient identiques. modèles de processeur, systèmes d'exploitation et langages de programmation. (Imaginez à quel point il serait amusant de déboguer votre code si les erreurs d'arrondi étaient différentes sur différents ordinateurs.) En raison des exigences d'arrondi, les malentendusrécemment mentionnéssurviennent. sur Habré: la fraction décimale 0,1 est représentée comme la fraction binaire la plus proche $$$7205759403792794\cdot2^{-56}$ , qui est exactement 0,1000000000000000055511151231257827021181583404541015625; 0,2 - sous la forme $$$7205759403792794\cdot2^{-55}$ , qui est exactement 0,200000000000000011102230246251565404236316680908203125; mais leur somme n'est pas une fraction binaire la plus proche de 0,3: approximation par le bas $$$5404319552844595\cdot2^{-54}$ = 0. 2999999999999988897769753748434595763683319091796875 se révèle être plus précis. Par conséquent, la norme IEEE nécessite d'ajouter 0,1 + 0,2 pour produire un résultat autre que 0,3.

Analyse double

Une généralisation simple de l'algorithme des nombres entiers consiste en des multiplications et des divisions itératives par 10,0 - contrairement à l'analyse des entiers, ici il est nécessaire de traiter les chiffres de bas en haut pour que les erreurs d'arrondi dans les chiffres élevés «masquent» les erreurs d'arrondi dans les chiffres faibles. En même temps, l'analyse de la mantisse se réduit facilement à l'analyse d'entiers: il suffit de changer la normalisation pour que le point binaire imaginaire ne soit pas au début de la mantisse, mais à la fin; l'entier de 53 bits résultant se multiplie ou se divise par la puissance désirée de dix; et enfin, soustrayez 52 de l'exposant, i.e. déplacez le point de 52 bits vers la gauche - où il devrait être selon la norme IEEE. De plus, il y a trois faits importants à noter:

1. Il suffit d'apprendre à multiplier ou à diviser par 5, et la multiplication ou la division par 2 n'est qu'un incrément ou un décrément d'un exposant;
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 complexité de l'arrondi standard est que pour découvrir que le résultat est exactement au milieu entre les fractions binaires les plus proches, nous n'avons pas seulement besoin de 54 bits significatifs du résultat (le bit zéro le moins significatif signifie "tout est en ordre", l'un signifie "frapper le milieu"), mais et vous devez regarder s'il y a eu une continuation non nulle après le bit le moins significatif. En particulier, cela signifie que ne considérer que les 17 chiffres les plus significatifs dans la notation décimale du nombre ne suffit pas: ils ne définissent la mantisse binaire qu'avec une précision de ± 1, et pour sélectionner le sens d'arrondi souhaité, vous devrez considérer les chiffres inférieurs. Par exemple, 10000000000000003 et 10000000000000005 sont la même valeur double (égale à 10000000000000004) et 10000000000000005.00 (plusieurs zéros) 001 est une valeur différente (égale à 10000000000000006).



Le professeur Daniel Lemire de l'Université par correspondance du Québec (TÉLUQ), qui a inventé cet analyseur, l'a publié sur github . Cette bibliothèque fournit deux fonctions from_chars



: l'une analyse une chaîne pour doubler, l'autre pour flotter. Ses expériences ont montré que dans 99,8% des cas, il suffit de considérer 19 chiffres décimaux significatifs (64 bits): si pour deux valeurs consécutives d'une mantisse de 64 bits, la même valeur double finale est obtenue, alors c'est la valeur correcte . Seulement dans 0,2% des cas, lorsque ces deux valeurs ne coïncident pas, un code plus complexe et plus universel est exécuté qui implémente toujours un arrondi correct.