02 sept. 2013

Racine d'un entier 64 bits non signé en C

Ce billet est un peu spécial, dans le sens ou je me contente de poster un morceau de code sans beaucoup de commentaires. J'ai écrit ce code récemment pour résoudre un problème un peu agaçant : calculer une racine carrée sur un nombre entier 64 bits.

Pour ceux qui se disent que le problème est simple (appeler la fonction sqrt() de la librairie C ou std::sqlt() de la librairie standard C++), je répond "oui, presque". Parce qu'un nombre non signé de 64 bits ne peut étre représenté correctement sur un double de 64 bits (sqrt() utilise le co-processeur arithmétique de votre x86). On peut utiliser un long double - mais dans ce cas, les performances sont désastreuses (parce que votre processeur n'a pas de registre générique travaillant en 80 bits).

Tout ça pour dire que la fonction ci-dessous peut aussi vous être utile (même si elle n'est pas parfaite : elle dépends principalement de la réponse à une devinette : quelle est la valeur de la racine carrée que je cherche ? Plus je sais répondre précisément à cette devinette, plus l'exécution est rapide).

Le code

#include <stdint.h>
static inline uint64_t newton_raphson(uint64_t prev, uint64_t number) { return (prev >> 1) + (number / (prev << 1)); }
uint64_t sqrt64(uint64_t number) { uint64_t current; uint64_t prev; uint64_t distance;
if (!number) return 0;
/* * gross approximation ; depends on sizeof(number) because we * are using some gcc builtins (which maps to ASM instructions * in this case) */ int hob = 8 * sizeof(number) - __builtin_clzll(number); current = 1ULL << (hob >> 1);
do { prev = current; current = newton_raphson(current, number); distance = prev < current ? current - prev : prev - current; } while (distance > 1);
return current; }

L'explication sommaire

La méthode générique est relativement simple - c'est l'application pure et simple de la convergence de suite dite de Newton-Raphson (d'où le nom de la fonction) :

x[n] = x[n-1] - f(x[n-1]) / f'(x[n-1])

Pour la racine carrée, la fonction f(x) = x^2 - k, donc f'(x) = 2x, k étant la valeur dont on souhaite obtenir la racine carrée.

     xn = x[n-1] - (x[n-1] * x[n-1] - k) / 2 * x[n-1]
==>  x[n] = x[n-1] - (x[n-1] / 2) + (k / 2 * x[n-1])
==>  x[n] = x[n-1] / 2 + (k / 2 * x[n-1])

Le problème avec cette méthode est que si la première valeur de la suite (x[0]) est mal devinée, on peut faire pas mal de cycles. Dans le code ci-dessous, on fait très simple (et - de fait - très mal).

L'algorithme peut être adapté à tout type d'entiers (128, 256 bits) Par contre, il peut être utile d'implémenter une fonction de grossière approximation d'une racine carrée n'un nombre entier (dans l'idéal, il faut que ça se fasse en une, deux opérations maximum, sans quoi ça va rapidement perdre de l'intérêt).

Dans le code d'exemple, on utilise la builting gcc __builtin_clzll() (count leading zeros). On se base sur le fait que quelque soit n, il existe un entier de la forme 2^k tel que n >= 2^k. Hors la racine carrée de ce nombre est égale à 2^(k/2). Du coup, on a une première valeur approximée de sqrt(n) en prenant 2^(k/2) (avec k^2 arrondi à l'entier inférieur). Avec cette approximation grossière, on obtient la racine carrée en maximum 6 itérations NR (sur un Core i7, ça représente environ 200 ns).

Je ne prétends pas que ce code soit véritablement production ready. Il est très certainement possible de l'optimiser pour obtenir de meilleures performances.

Pour information, avec une version récente de gcc, on obtient le code assembleur suivant (compilé en -O3 en x86_64):

00000000004006d0 <sqrt64>:
  4006d0:    31 c0                    xor    %eax,%eax
  4006d2:    48 85 ff                 test   %rdi,%rdi
  4006d5:    74 56                    je     40072d <sqrt64+0x5d>
  4006d7:    48 0f bd c7              bsr    %rdi,%rax
  4006db:    b9 40 00 00 00           mov    $0x40,%ecx
  4006e0:    be 01 00 00 00           mov    $0x1,%esi
  4006e5:    48 83 f0 3f              xor    $0x3f,%rax
  4006e9:    29 c1                    sub    %eax,%ecx
  4006eb:    d1 f9                    sar    %ecx
  4006ed:    48 d3 e6                 shl    %cl,%rsi
  4006f0:    eb 09                    jmp    4006fb <sqrt64+0x2b>
  4006f2:    66 0f 1f 44 00 00        nopw   0x0(%rax,%rax,1)
  4006f8:    48 89 d6                 mov    %rdx,%rsi
  4006fb:    48 8d 0c 36              lea    (%rsi,%rsi,1),%rcx
  4006ff:    31 d2                    xor    %edx,%edx
  400701:    48 89 f8                 mov    %rdi,%rax
  400704:    49 89 f0                 mov    %rsi,%r8
  400707:    48 f7 f1                 div    %rcx
  40070a:    49 d1 e8                 shr    %r8
  40070d:    48 89 f1                 mov    %rsi,%rcx
  400710:    4a 8d 14 00              lea    (%rax,%r8,1),%rdx
  400714:    48 89 d0                 mov    %rdx,%rax
  400717:    48 29 d1                 sub    %rdx,%rcx
  40071a:    48 29 f0                 sub    %rsi,%rax
  40071d:    48 39 d6                 cmp    %rdx,%rsi
  400720:    48 0f 42 c8              cmovb  %rax,%rcx
  400724:    48 83 f9 01              cmp    $0x1,%rcx
  400728:    77 ce                    ja     4006f8 <sqrt64+0x28>
  40072a:    48 89 d0                 mov    %rdx,%rax
  40072d:    f3 c3                    repz retq
  40072f:    90                       nop

Ajouter un commentaire

Les commentaires peuvent être formatés en utilisant une syntaxe wiki simplifiée.

Fil des commentaires de ce billet