Angenommen, es sind korrekt gerundete Standardbibliotheksfunktionen verfügbar, wie sie in CRlibm zu finden sind. Wie würde man dann die korrekt gerundete Kubikwurzel einer Eingabe mit doppelter Genauigkeit berechnen?
Diese Frage ist kein "tatsächliches Problem, dem [ich] gegenüberstehe", um die FAQ zu zitieren. Es ist ein bisschen wie Hausaufgaben auf diese Weise. Aber die kubische Wurzel ist eine häufig vorkommende Operation und man könnte sich vorstellen, dass diese Frage ein tatsächliches Problem ist, dem jemand gegenübersteht.
Da "die besten Stack-Überlauf-Fragen ein bisschen Quellcode enthalten", hier ein bisschen Quellcode:
%Vor% Das obige berechnet keine korrekt gerundete kubische Wurzel, weil 1/3 nicht genau als double darstellbar ist.
ZUSÄTZLICHE HINWEISE:
Ein Artikel beschreibt, wie man eine kubische Gleitkomma-Wurzel berechnet, aber die letzte Iteration ( s) des empfohlenen Newton-Raphson-Algorithmus müsste mit höherer Genauigkeit durchgeführt werden, damit der Algorithmus eine korrekt gerundete kubische Wurzel mit doppelter Genauigkeit berechnet. Das mag der beste Weg sein, um es zu berechnen, aber ich bin immer noch auf der Suche nach einer Abkürzung, die vorhandene korrekt gerundete standardisierte Funktionen nutzen würde.
C99 enthält eine Funktion cbrt() , aber es kann nicht erwartet werden, dass sie korrekt gerundet wird oder sogar treu für alle Compiler. Die Designer von CRlibm hätten auswählen können, dass cbrt() in die Liste der bereitgestellten Funktionen aufgenommen wird, was sie jedoch nicht taten. Verweise auf Implementierungen, die in anderen Bibliotheken mit korrekt gerundeten mathematischen Funktionen verfügbar sind, sind willkommen.
Ich fürchte, ich weiß nicht, wie man eine korrekt gerundete doppelpräzise Kubikwurzel garantiert, sondern kann eine solche anbieten, die sehr genau richtig gerundet ist, wie auch in der Frage erwähnt. Mit anderen Worten, der maximale Fehler scheint sehr nahe bei 0,5 ul zu liegen.
Peter Markstein, "IA-64 und Elementarfunktionen: Geschwindigkeit und Präzision" (Prentice-Hall 2000)
stellt effiziente FMA-basierte Techniken zum korrekten Runden der reziproken, der Quadratwurzel und der reziproken Quadratwurzel dar, deckt jedoch in dieser Hinsicht nicht die Kubikwurzel ab. Im Allgemeinen erfordert Marksteins Ansatz ein vorläufiges Ergebnis, das innerhalb von 1 ulp vor der abschließenden Rundungssequenz korrekt ist. Ich habe nicht die mathematischen Mittel, um seine Technik auf die Rundung von Würfelwurzeln auszudehnen, aber es scheint mir, dass dies prinzipiell möglich sein sollte, eine Herausforderung, die der reziproken Quadratwurzel etwas ähnlich ist.
Bit-weise Algorithmen eignen sich leicht für die Berechnung von Wurzeln mit korrekter Rundung. Da Bindefälle für den IEEE-754-Rundungsmodus "am nächsten" oder "gerade" nicht auftreten können, muss lediglich die Berechnung durchgeführt werden, bis sie alle Mantissenbits plus ein rundes Bit erzeugt hat. Der bitweise Algorithmus für die Quadratwurzel, basierend auf dem Binomialsatz, ist sowohl in nicht-wiederherstellenden als auch in wiederherstellenden Varianten gut bekannt und war die Grundlage von Hardware-Implementierungen. Derselbe Ansatz über Binomialsatz funktioniert für die Kubikwurzel, und es gibt ein wenig bekanntes Papier, das die Details einer nicht wiederherstellenden Implementierung beschreibt:
H. Peng, "Algorithmen für die Extraktion von Quadratwurzeln und Würfelwurzeln", Proceedings 5. IEEE International Symposium on Computer Arithmetic, pp. 121-126, 1981.
Am besten kann ich sagen, dass es beim Experimentieren gut genug ist, um Kubikwurzeln aus ganzen Zahlen zu extrahieren. Da jede Iteration nur ein einziges Ergebnisbit erzeugt, ist sie nicht genau schnell. Für Anwendungen in der Gleitkomma-Arithmetik hat sie den Nachteil, einige Buchführungsvariablen zu verwenden, die ungefähr die doppelte Anzahl von Bits des Endergebnisses erfordern. Dies bedeutet, dass eine 128-Bit-Integer-Arithmetik verwendet werden muss, um einen Würfelwürfel mit doppelter Genauigkeit zu implementieren.
Mein C99-Code basiert auf Halley's rationaler Methode für die Kubikwurzel Die kubische Konvergenz bedeutet, dass die anfängliche Approximation nicht sehr genau sein muss, da die Anzahl der gültigen Zahlen in jeder Iteration verdreifacht. Die Berechnung kann auf verschiedene Arten erfolgen. Im Allgemeinen ist es von numerischem Vorteil, iterative Schemata als
anzuordnennew_guess: = old_guess + Korrektur
da correction für eine genügend enge Anfangsschätzung signifikant kleiner ist als old_guess . Dies führt zu dem folgenden Iterationsschema für die Kubikwurzel:
x: = x - x * (x 3 - a) / (2 * x 3 + y)
Diese besondere Anordnung ist auch in Kahans Notizen zur Kubikwurzel aufgeführt. Es hat den weiteren Vorteil, sich natürlich der Verwendung von FMA (fusion-multiply_operation) zu bedienen . Ein Nachteil besteht darin, dass die Berechnung von 2 · x 3 zu einem Überlauf führen könnte, weshalb ein Argumentreduktionsschema für mindestens einen Teil der Eingangsdomäne mit doppelter Genauigkeit erforderlich ist. In meinem Code verwende ich einfach die Argumentreduktion auf alle nicht-außergewöhnlichen Eingaben, basierend auf einer einfachen Manipulation der Exponenten von IEEE-754-Operanden mit doppelter Genauigkeit.
Das Intervall [0,125,1) dient als primäres Näherungsintervall. Eine Polynom-Minimax-Approximation wird verwendet, die eine anfängliche Schätzung in [0,5,1] zurückgibt. Der enge Bereich erleichtert die Verwendung von Berechnungen mit einfacher Genauigkeit für die Abschnitte mit geringer Genauigkeit der Berechnung.
Ich kann nichts über die Fehlergrenzen meiner Implementierung beweisen, aber das Testen mit mehr als 200 Millionen zufälligen Testvektoren gegen eine Referenzimplementierung (genau auf etwa 200 Bit) hat keine Fehler über eine halbe ul erkannt, was darauf hindeutet, dass Der maximale Fehler muss sehr nahe bei 0,5 ul liegen.
%Vor%Da es viele leicht berechenbare rationale Punkte auf der Kurve x = y ^ 3 gibt, bin ich versucht, ungefähr s ^ 3 ~ x zu reduzieren, mit s rational und nur ein paar Bits breit. Dann hast du:
%Vor%Offensichtlich ist es dann, den Korrekturterm mit Hilfe Ihrer bevorzugten Serienapproximation zu berechnen und eine Residualdarstellung über Head-Tail-FMA zu berechnen, um das Ergebnis bei Bedarf nach oben oder unten zu bringen (Sie werden nicht die volle Berechnung benötigen der Zeit, offensichtlich).
Das ist nicht ganz im Geiste der Frage, aber es kann definitiv zum Funktionieren gebracht werden, und es ist sehr einfach, die notwendigen Grenzen auf diese Weise zu beweisen. Hoffentlich kann jemand anderes etwas schlauer vorschlagen (ich habe meine Schlauheit bereits für den Monat verbraucht).
Tags und Links algorithm floating-point ieee-754