Round y = x * x al più vicino

Question

Arrotondare il risultato della divisione al numero intero più vicino è pretty simple. Ma sto cercando di arrotondare un risultato della divisione in modo che l'operazione successiva dia la migliore approssimazione. Che è meglio spiegata con una semplice funzione:

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. 

posso arrotondare x_div al più vicino con l'aggiunta di 1<<(halfbits-1). Ma poiché x² non è una funzione lineare, y in generale non è arrotondato correttamente. Esiste un modo semplice e più accurato per calcolare (x*x) >> (halfbits*2) senza utilizzare tipi più grandi?

Penso che l'aggiunta di 3<<(halfbits-3) a x_div migliora l'arrotondamento, ma non può dimostrare che sia la soluzione migliore. Inoltre, può essere generalizzato per xⁿ?

Modifica: a grande richiesta mi prendo la libertà di "tradurre" la domanda in termini aritmetici puri (nessuna di queste cose che spostano il bit C ...).
Nota: tutte le divisioni successive sono le divisioni di interi, ad esempio 13/3 sarebbe 4.
Problema: non possiamo calcolare x^2 perché x è grande, quindi vorremmo calcolare (x^2)/(2^N).
A questo scopo si calcola
x_div = x/sqrt (2^N)
che poi quadrato:
y = x_div * x_div

Tuttavia, questo risultato è tipicamente breve del valore esatto di (x^2)/(2^N) e l'OP suggerisce di aggiungere 0,5 * sqrt (2^N) o forse 0,375 * sqrt (2^N) per approssimare meglio il risultato ...

come risposta Oli Charlesworth come suggerisce c'è un modo molto migliore per arrivare al valore effettivo, pensando a x^2 come (x_hi + x_lo)^2.

Answer 1

Considerando che con il troncamento x_div causerà un errore di grandezza al massimo 1, con x_div*x_div l'errore potrebbe essere fino a 1<<(half_digits+2).

capire perché, consideriamo che possiamo esprimere questa quadratura (usando long multplication) come segue:

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

dove x_lo e x_hi sono le metà inferiore e superiore x, rispettivamente. Con qualche bella arte ASCII, possiamo considerare come questi tutti in fila:

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

Come possiamo vedere, i termini di ordine inferiore colpiscono più bit nel risultato finale.

Tuttavia, ognuno di questi termini deve rientrare nel tipo originale senza overflow. Quindi, con il cambio e il mascheramento appropriati, possiamo ottenere il risultato desiderato.

Naturalmente, il compilatore/hardware fa tutto questo per te se usi un tipo più grande, quindi dovresti farlo solo se ne hai l'opzione.

Answer 2

Utilizzare TYPE int per "y". Immagino che risolverebbe lo scopo.