sqrt, perfekte Quadrate und Gleitkommafehler

7

In der Funktion sqrt der meisten Sprachen (obwohl ich mich hier hauptsächlich für C und Haskell interessiere), gibt es irgendwelche Garantien, dass die Quadratwurzel eines perfekten Quadrats genau zurückgegeben wird? Zum Beispiel, wenn ich sqrt(81.0) == 9.0 mache, ist das sicher oder gibt es eine Chance, dass sqrt 8.999999998 oder 9.00000003 zurückgibt?

Wenn die numerische Genauigkeit nicht garantiert ist, was wäre der beste Weg, um zu überprüfen, ob eine Zahl ein perfektes Quadrat ist? Nehmen Sie die Quadratwurzel, erhalten Sie den Boden und die Decke und stellen Sie sicher, dass sie sich auf die ursprüngliche Zahl zurücksetzen.

Danke!

    
gnuvince 04.01.2013, 05:45
quelle

4 Antworten

12

Wenn in IEEE 754 der Gleitkommawert mit doppelter Genauigkeit x das Quadrat einer nichtnegativ darstellbaren Zahl y ist (dh y * y == x und die Berechnung von y * y keine Rundung, Überlauf oder Unterlauf), dann wird sqrt (x) y zurückgeben.

Das ist alles, weil sqrt vom Standard IEEE 754 korrekt gerundet werden muss. Das heißt, sqrt (x) für any x ist das nächstliegende Doppel der tatsächlichen Quadratwurzel von x. Dieser sqrt funktioniert für perfekte Quadrate ist eine einfache Folge dieser Tatsache.

Wenn Sie überprüfen möchten, ob eine Fließkommazahl ein perfektes Quadrat ist, ist hier der einfachste Code, den ich mir vorstellen kann:

%Vor%

Ich brauche den leeren asm volatile -Block, der von dd abhängt, weil dein Compiler ansonsten clever sein könnte und die Berechnung von dd "wegoptimiere".

Ich habe einige seltsame Funktionen aus fenv.h verwendet, nämlich feclearexcept und fetestexcept . Es ist wahrscheinlich eine gute Idee, sich ihre man Seiten anzusehen.

Eine andere Strategie, die Sie möglicherweise in der Lage sind, Arbeit zu machen, besteht darin, die Quadratwurzel zu berechnen, zu prüfen, ob sie Bits in den unteren 26 Bits der Mantisse gesetzt hat, und sich zu beschweren, wenn dies der Fall ist. Ich versuche diesen Ansatz unten.

Und ich musste überprüfen, ob d 0 ist, weil sonst true für -0.0 zurückgegeben werden kann.

BEARBEITEN : Eric Postpischil hat vorgeschlagen, dass das Hacken mit der Mantisse besser sein könnte. Da das obige issquare in einem anderen populären Compiler, clang , nicht funktioniert, stimme ich dem zu. Ich denke der folgende Code funktioniert:

%Vor%

Das Hinzufügen und Subtrahieren von 67108864.0 von a hat den Effekt, die niedrigen 26 Bits der Mantisse zu löschen. Wir werden a genau zurückbekommen, wenn diese Bits überhaupt klar waren.

    
tmyklebu 04.01.2013, 07:28
quelle
6

Nach diesem Papier , in dem der Nachweis der Korrektheit von IEEE-Gleitkommazahlen diskutiert wird Quadratwurzel:

  

Der IEEE-754-Standard für binäre Fließkommazahl   Arithmetik [1] erfordert das Ergebnis einer Division oder eines Quadrats   Wurzeloperation wird berechnet als ob in unendlicher Genauigkeit und   dann zu einem der beiden nächsten Fließkommawerte gerundet   Nummern der angegebenen Genauigkeit, die die   unendlich präzises Ergebnis

Da ein perfektes Quadrat, das genau im Fließkomma dargestellt werden kann, eine ganze Zahl ist und seine Quadratwurzel eine Ganzzahl ist, die genau dargestellt werden kann, sollte die Quadratwurzel eines perfekten Quadrats immer genau richtig sein.

Natürlich gibt es keine Garantie, dass Ihr Code mit einer konformen IEEE Fließkomma-Bibliothek ausgeführt wird.

    
Gabe 04.01.2013 07:32
quelle
1

@tmyklebu hat die Frage perfekt beantwortet. Lassen Sie uns als Ergänzung eine möglicherweise weniger effiziente Alternative zum Testen des perfekten Quadrats von Fraktionen ohne asm-Richtlinie betrachten.

Nehmen wir an, wir haben ein IEEE 754-konformes sqrt, das das Ergebnis korrekt rundet.
Nehmen wir an, dass außergewöhnliche Werte (Inf / Nan) und Nullen (+/-) bereits behandelt werden.
Lassen Sie uns sqrt(x) in I*2^m zerlegen, wobei I eine ungerade ganze Zahl ist.
Und wo I über n Bits geht: 1+2^(n-1) <= I < 2^n .

Wenn n > 1+floor(p/2) wobei p Gleitkomma-Genauigkeit ist (z. B. p = 53 und n & gt; 27 in doppelter Genauigkeit)
Dann 2^(2n-2) < I^2 < 2^2n .
Da I ungerade ist, ist I^2 auch ungerade und überspannt daher & gt; p Bits.
Daher ist I nicht die exakte Quadratwurzel eines darstellbaren Fließkommas mit dieser Genauigkeit.

Aber gegeben I^2<2^p , könnten wir sagen, dass x ein perfektes Quadrat war?
Die Antwort ist offensichtlich nein. Eine Taylor-Erweiterung würde

geben %Vor%

Daher wird für e=ulp(I^2) bis sqrt(ulp(I^2)) die Quadratwurzel korrekt auf rsqrt(I^2+e)=I ... gerundet (rund zum nächsten geraden oder abgeschnittenen oder Bodenmodus).

Also müssten wir das sqrt(x)*sqrt(x) == x bestätigen.
Aber oben Test ist nicht ausreichend, zum Beispiel unter der Annahme, IEEE 754 mit doppelter Genauigkeit, sqrt(1.0e200)*sqrt(1.0e200)=1.0e200 , wobei 1.0e200 genau ist 99999999999999996973312221251036165947450327545502362648241750950346848435554075534196338404706251868027512415973882408182135734368278484639385041047239877871023591066789981811181813306167128854888448 deren erster Primfaktor 2^613 , kaum ein perfektes Quadrat jeder Fraktion ...

So können wir beide Tests kombinieren:

%Vor%

BEARBEITEN: Wenn wir uns auf den Fall von Ganzzahlen beschränken wollen, können wir auch überprüfen, ob floor(sqrt(x))==sqrt(x) oder dreckige Bit-Hacks in squared_significand_fits_in_precision verwenden ...

    
aka.nice 04.01.2013 23:14
quelle
0

Anstatt sqrt(81.0) == 9.0 zu versuchen, probiere 9.0*9.0 == 81.0 . Dies funktioniert immer so lange, wie das Quadrat innerhalb der Grenzen der Fließkommagröße liegt.

Bearbeiten: Ich war wahrscheinlich nicht sicher, was ich unter "Fließkomma-Größe" verstehe. Was ich meine ist, die Zahl innerhalb des Bereichs von ganzzahligen Werten zu halten, die ohne Präzisionsverlust gehalten werden können, weniger als 2 ** 53 für ein IEEE-Double. Ich erwartete auch, dass es eine separate Operation geben würde, um sicherzustellen, dass die Quadratwurzel eine Ganzzahl ist.

%Vor%     
Mark Ransom 04.01.2013 05:50
quelle

Tags und Links