Arrondi y = x * x au plus proche

Question

Arrondir le résultat de la division à l'entier le plus proche est pretty simple. Mais j'essaie d'arrondir un résultat de division pour que l'opération suivante donne la meilleure approximation. Ce qui est le mieux expliqué par une fonction simple:

const int halfbits = std::numeric_limits<unsigned int>::digits/2; 
unsigned x = foo(); // Likely big. 
unsigned x_div = x >> halfbits; // Rounded down 
unsigned y = x_div * x_div; // Will fit due to division. 

Je peux arrondir x_div à la plus proche en ajoutant 1<<(halfbits-1). Mais puisque x² n'est pas une fonction linéaire, y en général n'est pas arrondi correctement. Existe-t-il un moyen simple et plus précis de calculer (x*x) >> (halfbits*2) sans utiliser de plus grands types?

Je pense qu'ajouter 3<<(halfbits-3) à x_div améliore l'arrondi, mais ne peut pas prouver que c'est la meilleure solution. Aussi, cela peut-il être généralisé pour xⁿ?

Édition: à la demande générale, je prends la liberté de "traduire" la question en termes purement arithmétiques (aucune de ces choses de décalage en bits C ...).
Note: Toutes les divisions suivantes sont des divisons entières, par exemple 13/3 serait 4.
Problème: Nous ne pouvons pas calculer x^2 parce que x est grand, donc nous aimerions calculer (x^2)/(2^N).
Pour ce faire, nous calculons
x_div = x/sqrt (2^N)
que nous avons ensuite place:
y = x_div * x_div

Cependant, ce résultat est généralement peu de la valeur exacte de (x^2)/(2^N) et l'OP suggère d'ajouter 0,5 * sqrt (2^N) ou peut-être 0,375 * sqrt (2^N) pour mieux approximer le résultat ...

Comme Oli Charlesworth La réponse suggère qu'il y a une bien meilleure façon d'arriver à valeur réelle, en pensant à x^2 as (x_hi + x_lo)^2.

Answer 1

Considérant qu'avec x_div troncature provoquera une erreur de magnitude d'au plus 1, avec x_div*x_div l'erreur pourrait aller jusqu'à 1<<(half_digits+2).

Pour comprendre pourquoi, considérons que nous pouvons exprimer cette élévation au carré comme suit (en utilisant long multplication):

x * x = (x_lo + x_hi) * (x_lo + x_hi) 
     = x_lo^2 + x_hi^2 + 2*x_lo*x_hi 

où x_lo et x_hi sont les moitiés inférieure et supérieure de x, respectivement. Avec un peu d'art ASCII agréable, on peut considérer la façon dont ces toutes les lignes jusqu'à:

MSB  :  :  :  LSB 
    +------+------+  :  : 
    | x_hi^2 |  :  : 
    +-----++-----++-----+:  : 
    :  | 2*x_lo*x_hi |:  : 
    :  +------++-----++------+ 
    :  :  | x_lo^2 | 
    :  :  +------+------+ 
    :<----------->:  :  : 
    Our result 

Comme on peut le voir, les termes d'ordre inférieur affectent plusieurs bits dans le résultat final. Toutefois, chacun de ces termes doit correspondre au type d'origine sans débordement. Donc, avec un décalage et un masquage appropriés, nous pouvons obtenir le résultat souhaité.

Bien sûr, le compilateur/matériel fait tout cela pour vous si vous utilisez un type plus grand, donc vous devriez le faire si vous avez l'option.

Answer 2

Utilisez TYPE int pour «y». Je suppose que cela permettrait de résoudre le problème.