Rund y = x * x zum nächsten

Question

Runden des Ergebnisses der Division auf die nächste ganze Zahl ist pretty simple. Aber ich versuche ein Ergebnis der Division zu runden, damit die nachfolgende Operation die beste Annäherung ergibt. Das ist am besten erklärt durch eine einfache Funktion:

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. 

I x_div zum nächsten Runde kann 1<<(halfbits-1) durch Zugabe. Aber da x² keine lineare Funktion ist, wird y im Allgemeinen nicht richtig gerundet. Gibt es eine einfache und genauere Möglichkeit, (x*x) >> (halfbits*2) zu berechnen, ohne größere Typen zu verwenden?

Ich denke, dass das Hinzufügen von 3<<(halfbits-3) zu x_div die Rundung verbessert, kann aber nicht beweisen, dass das die beste Lösung ist. Kann dies auch für xⁿ verallgemeinert werden?

Bearbeiten: auf vielfachen Wunsch erlaube ich mir, die Frage in rein rechnerischen Begriffen "übersetzen" (keine dieser C Bit-Verschiebung Sache ...).
Hinweis: Alle Divisionen danach sind ganzzahlige Divisionen, zum Beispiel 13/3 wäre 4.
Problem: Wir können nicht x^2 berechnen, weil x groß ist, also möchten wir stattdessen berechnen (x^2)/(2^N).
Dazu berechnen wir
x_div = x/sqrt (2^N)
die wir dann Quadrat:
y = x_div * x_div

Dies ist jedoch aufgrund der Regel kurz vor dem exakten Wert von (x^2)/(2^N) und das OP schlägt vor, 0.5 * sqrt (2^N) oder vielleicht 0.375 * sqrt (2^N) hinzuzufügen, um das Ergebnis besser zu nähern ...

Oli Charlesworth 's Antwort schlägt vor, es gibt eine viel bessere Weise, zum zu gelangen tatsächlicher Wert, indem an x^2 als (x_hi + x_lo)^2 gedacht wird.

Answer 1

Während mit x_div Kürzung eine Fehlergröße von höchstens 1 verursacht, mit x_div*x_div könnte der Fehler bis zu 1<<(half_digits+2) sein.

Um zu sehen, warum der Ansicht, dass wir diese Quadratur ausdrücken kann wie folgt (mit long multplication):

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

wo x_lo und x_hi sind die unteren und oberen Hälften x sind. Mit einigen netten ASCII-Art, können wir überlegen, wie diese alle Line-Up:

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

Wie wir sehen können, die niedriger Ordnung Bedingungen mehrere Bits in das Endergebnis beeinflussen.

Allerdings sollte jeder dieser Begriffe in den ursprünglichen Typ ohne Überlauf passen. Mit entsprechenden Verschiebungen und Maskierungen können wir das gewünschte Ergebnis erzielen.

Natürlich tut der Compiler/Hardware das alles für Sie, wenn Sie einen größeren Typ verwenden, also sollten Sie das tun, wenn Sie die Option haben.

Answer 2

Verwenden Sie TYPE int für "y". Ich denke, das würde den Zweck lösen.