Hypot entier sans débordement

Question

Comment calculer hypot de deux entiers, chacun inférieur à 2^63, de sorte qu'à aucun moment aucun calcul intermédiaire ne déborde 64 bits? (comme x^2+y^2 dans l'approche traditionnelle).

L'article lié mentionne un algorithme à virgule flottante, qu'il n'est pas possible d'utiliser à cause de cela t = t/x; car il est 0 pour les entiers.

L'algorithme la plus approchante est de here mais malheureusement il est pas assez précis:

int ihypot(xd1, yd1) 
    double xd1, yd1; 
{ 
    register  x1 = (int)xd1, 
       y1 = (int)yd1, 
       x2 = 0, 
       y2 = 0; 

    if ((x2 -= x1) < 0) x2 = -x2; 
    if ((y2 -= y1) < 0) y2 = -y2; 
    return (x2 + y2 - (((x2>y2) ? y2 : x2) >> 1)); 
} 

Answer 1

Vous pouvez toujours commencer à partir de la formule mathématique équivalente mise en œuvre avec des nombres entiers:

r = 2 ^-30 * (x * sqrt (2 + 2 * y/x))

A Le processeur 32 bits typique devrait vous donner accès à une division 64/32 -> 32 avec l'entrée fournie dans deux registres. Cette division peut être utilisée pour calculer 2 * y/x. Votre langage de programmation peut ou non vous donner accès à celui-ci. Ne sous-estimez pas la compétence d'optimisation des compilateurs lors de la génération de code 32 bits pour les calculs impliquant des résultats intermédiaires de 64 bits.

De même, le processeur 32 bits typique devrait fournir une multiplication 32 * 32 -> 64 pour "x * ...", avec pour résultat deux registres. La multiplication finale par 2 ^-30 équivaut à déplacer et/ou orienter les deux registres dans lesquels la multiplication 32 * 32 -> 64 a quitté son résultat.

GCC almost parvient à générer un code simple en utilisant uniquement des instructions 32 bits, mais il laisse tomber la balle à un point et appelle une fonction externe de la division multi-précision:

#include <stdint.h> 

uint32_t integer_sqrt(uint32_t); 

/*@ requires x >= y; */ 
uint32_t hypot(uint32_t x, uint32_t y){ 
    return integer_sqrt(0x40000000 + (uint32_t) ((uint64_t)y * 0x40000000/x)) * (uint64_t) x/0x40000000 ; 
} 

résultat de l'assemblage 32 bits:

hypot: 
     pushl %edi 
     pushl %esi 
     xorl %edi, %edi 
     pushl %ebx 
     movl 16(%esp), %ebx 
     movl %edi, %edx 
     xorl %edi, %edi 
     subl $16, %esp 
     movl 36(%esp), %esi 
     pushl %edi 
     pushl %ebx 
     shldl $30, %esi, %edx 
     movl %esi, %eax 
     sall $30, %eax 
     pushl %edx 
     pushl %eax 
     call __udivdi3 
     addl $20, %esp 
     addl $1073741824, %eax 
     pushl %eax 
     call integer_sqrt 
     mull %ebx 
     addl $16, %esp 
     popl %ebx 
     popl %esi 
     shrdl $30, %edx, %eax 
     popl %edi 
     ret 

Edit:

Si vous souhaitez utiliser uniquement 32 * 32-> 32 multiplication, vous devrez calculer N = log2 (x) et, si N> 15, normaliz ex et y en les décalant à droite par N-15 (décalage du résultat d'extrémité gauche de la même valeur), en effet la mise en oeuvre de formule:

r = 2 ^N-15 * sqrt ((x/2 ^{N -15}) + (y/2 ^N-15))

Si N≤1, il suffit d'utiliser la formule habituelle r = sqrt (x ^{2 + y2)}

Answer 2

Je pense que vous pouvez réécrire l'expression donnée de sorte que chaque après chaque opération soit inférieure à 2^64.Par exemple:

En prenant la racine carrée avant la multiplication, vous pouvez vous assurer que les numéros individuels ne dépassent jamais 2^64.

Answer 3

Si vous êtes d'accord avec l'utilisation de virgule flottante, vous pouvez utiliser l'expression

a/cos(atan2(b,a)) 

où a et b sont les entrées