Welchen Algorithmus für reele Wurzeln?

Hallo,
nach welchem Algorithmus werden eigentlich von Rechnern Wurzeln der Form

$x = a^{b}$

berechnet? a, b, x sind reelle Zahlen.
Ich kann also a, b eintippen und erhalte x.
b kann ja auch z.B. die Zahl $$ \sqrt{2} $$ sein, usw.

SeppJ

Ein möglicher Quellcode mit Erklärung:

http://www.google.com/codesearch/p?hl=en#xQlobzjVCWM/e_pow.c&q=pow&sa=N&cd=3&ct=rc

Andere Implementierungen dürfen nicht großartig anders sein. Falls du an Alternativen oder Näherungsweisen Berechnungen interessiert bist, findet man auch jede Menge wissenschaftliche Artikel dazu. Du darfst das bloß nicht "reelle Wurzel" nennen. Das was du beschreibst sind einfach allgemein Potenzen beziehungsweise Exponentialfunktionen^*.

^*: Wissenschaftliche Artikel müssen natürlich mit englischen Suchbegriffen gesucht werden.

icarus2

Man kann die Quadratwurzel einer Zahl mit dem Heron-Verfahren approximieren. Ich denke allerdings nicht, dass dieser Algorithmus auf einem Taschenrechner oder so implementiert ist.

Ich glaube man könnte auch mit der Taylor-Reihe approximieren. Das wird bei sin und cos oft gemacht.
Ich sehe gerade keinen Grund, warum Taylor nicht funktionieren würde. Aber sicher ibn ich mir nicht, du müsstest es ausprobieren. Bitte korrigiert mich, falls jemand weiss, dass es mit Taylor nicht geht oder nicht sinnvoll ist.

SeppJ

Taylorreihen sind in praktischer Numerik nicht sinnvoll und ich kenne keine praktisch benutzte Implementierung irgendeiner Funktion über ihre Taylorrreihe. Sinus & Co. werden z.B. in der Regel über Tabellen oder CORDIC implementiert. Taylorreihen konvergieren nämlich in der Regel sehr, sehr schlecht, gerade der Sinus ist ein hervorragendes Beispiel wo man ewig hohe Ordnungen für akzeptable Genauigkeit benötigt. Taylorreihen sind eher etwas für theoretische Anwendungen, Anwendungen bei denen sie tatsächlich mal gut funktionieren (z.B. Sinus für kleine Winkel) und für Anwendungen wo man wirklich nur grobe Approximationen benötigt.

icarus2

Alles klar, danke für die Info.

SeppJ

Und noch etwas, was mir gerade einfällt: Die Taylorreihe der Quadratwurzel ist auch ein typisches Beispiel für ein Reihe die nicht überall konvergiert. Wenn ich mich recht entsinne, konvergiert sie nur zwischen 0 und 1, müsste ich jetzt nachrrechnen, wenn ich weniger faul wäre

icarus2

Ich wollte mal kurz das Heron-Verfahren ausprobieren und habs kurz implementiert.

#include <iostream>

int main()
{
	using std::cout;
	using std::cin;
	using std::endl;

	double a; // The number to calculate the square root
	double x; // The square root of a
	int number_of_iterations;

	cout << "Please enter the number to compute the square root: ";
	cin >> a;
	cout << endl;

	cout << "Please enter the number of iterations: ";
	cin >> number_of_iterations;
	cout << endl;

	x = a; // Choose the start x to be a (not suitable for large a)

	for( int i = 0; i < number_of_iterations; i++ )
	{
		x = ( x + a / x ) / 2; // Heron formula
	}

	cout.precision( 12 );
	cout << "The square root of a is: " << x << endl << endl;
}

Für Werte < 1000 kommt man bereits mit 10 Iterationen etwa auf 4-5 Nachkommastellen genau heran.
Für werte wie 200'000 braucht man dann halt schon so ca. 100 Iterationen.
Was helfen könnte wäre den Startwert für x besser zu wählen.

SeppJ

icarus2 schrieb:

Was helfen könnte wäre den Startwert für x besser zu wählen.

Da kann man überhaupt sehr clevere Sachen machen. Fürs Heron Verfahren weiß ich es konkret nicht, aber für das Newtonverfahren gibt es ungeheuer clevere Startwerte welche die Feinheiten der Computerarihmetik ausnutzen, z.B. der berühmte Code für 1/sqrt(x) (fälschlicherweise den Programmierern von Quake 3 zugeschrieben):
http://en.wikipedia.org/wiki/Fast_inverse_square_root

(Anmerkung: Als ich ihn vor ein wenigen Jahren versucht habe, war er langsamerr als die normale Quadratwurzel der libm, wenn man gleiche Genauigkeit wünschte (wozu 2 Iterationen nötig waren). Gut mögiich, dass meine libm mittlerweile selbst diesen Algorithmus implementiert hatte.^*)

^*: edit2: Nach Lesen meines eigenen Links , wird die libm wohl SSE Anweisungen benutzen und daher schneller undgleichzeitig genauer sein.

icarus2

Sieht interessant aus. Ich habe leider nicht alles verstehen können. Ich hab mir den Link mal gespeichert. Vielleicht dann nochmals durchlesen wenn ich Numerik gehabt habe ^^

krümelkacker

Zunächst kann man erstmal die "Domäne" verkleinern; denn

( x * 2^2k )^0.5 = x^0.5 * 2^k

Damit kann man das Problem, eine Wurzel einer beliebigen Zahl zu berechnen, auf die Wurzel einer Zahl zwischen 0.5 und 2 reduzieren. Die übliche Fließkommazahlendarstellung erlaubt schnelles Extrahieren von Exponenten und Multiplizieren mit Zweierpotenzen. Siehe std::frexp und std::ldexp.

double wurzel(double t)
{
  int exponent;
  h = frexp(t,&exponent);
  if (exponent & 1u) { // ungerade?
    h *= 2;
    exponent--;
  }
  // Hier gilt
  // 1.  0.5 <= h <= 2.0
  // 2.  h * pow(2,exponent) = t
  // 3.  exponent ist gerade
  double w = wurzel_bei_1(h);
  return std::ldexp(w,exponent/2);
}

Jetzt muss man sich nur noch um "wurzel_bei_1" kümmern. Als Nullstellenproblem formuliert: Wir suchen eine Lösung für ww-h=0. Da 0.5<=h<=2.0 gilt, ist w_est=1.0 gar keine so schlechte Näherung. Wir können aber die Wurzelfunktion bei 1 auch linearisieren und erhalten eine etwas bessere Schätzung: w_est=0.5(1+h) Diese können wir über ein paar Newton-Iterationen verfeinern.

double wurzel_bei_1(double h)
{
  double w = 0.5*(1+h); // Linearisierung der Wurzelfunktion bei 1
  for (int iter=0; iter<6; ++iter) {
    double y = w*w-h; // Polynom bei w auswerten
    double m = 2*w;   // Ableitung bei w auswerten
    w -= y/m;
  }
  return w;
}

Das Newton-Vervahren konvergiert quadratisch. Das heißt, mit jeder Iteration wird die Anzahl der korrekten Dezimalstellen verdoppelt. (Das ist gut so)

Leider ist die Division üblicherweise sehr langsam. Wenn man stattdessen 1/sqrt(h) iterativ berechnet, kommt man um die Division herum und man kann hinterher einmal den Kehrwert berechnen. Vielleicht gibt es auch noch einen anderen effizienteren Weg. Ohne Recherche komme ich nur auf diese Ansätze.