Wie bei Heronverfahren abbrechen?

SeppJ

krümelkacker schrieb:

Ich würde Deine Variante vorziehen, eben weil floatingpoint-Division typischerweise (?) 4mal langsamer ist als Multiplikation (zumindest hatte ich das mal so beobachtet auf einem core2 Rechner).

Und? Wieviel ist das schon gegen den Rest der Rechenschritte? Dafür beschränkst du dich auf maximal ca. 7 Stellen Genauigkeit von den ca. 15 möglichen bei double. Wenn dir 7 Stellen ausreichen, ist es natürlich ok, aber dann kannste auch gleich mit float rechnen (mit Division) und nicht den falschen Anschein von double-Genauigkeit erwecken.

Wenn man so viel Angst hat vor Division (4 mal länger sind sind immer noch bloß ca. 8 Taktzyklen) kann man ja auch am Anfang einmal 1/a ausrechnen.

rapso

HeronVerfahren schrieb:

//xn = aktueller wert, a= Zahl, aus der die Wurzel gezogen werden soll
if( abs((xn * xn) - a) <= prec) break;
Bei Wikipedia steht diese Abschätzung:
if( (xn - xn/a) <= prec) break;
Aber ist diese Variante nicht langsamer als meine naive? Denn hier wird eine Division statt einer Multiplikation durchgeführt.
Gibt es vielleicht bessere Methoden, den Fehler abzuschätzen?

wieso nicht

if( xn - xn <= prec*a) break;

?
sollte gleich funktionieren wie zweiteres und sich von der laufzeit verhalten wie ersteres.

eine division auf x87 sollte nicht viel kuerzer sein als eine wurzelberechnung. der gewinn ist vermutlich nichtig wenn du das so machst (ich wuerde wegen den RAW hazards eher vermuten dass es langsammer ist als fsqrt.

SeppJ

rapso schrieb:

wieso nicht

if( xn - xn <= prec*a) break;

?

Ich glaube du hast da einen kleinen Fehler drin...

krümelkacker

SeppJ schrieb:

Und? Wieviel ist das schon gegen den Rest der Rechenschritte?

30% oder mehr bei naiver Implementierung.

Da man a/x aber wieder in der nächsten Iteration braucht, muss man das ja nicht doppelt berechnen. Hier eine weniger naive Implementierung:

a = ...;
x = ...;
schwellwert = ...;
double t = a/x;
for (;;) {
  x = 0.5*(x+t);
  t = a/x;
  if (abs(t-x)<=schwellwert) break;
}
// x = 0.5*(x+t); // eventuell noch eine runde, um die korrekte zahl von stellen zu verdoppeln soweit möglich
return x;

Aber ich wüsste nicht, wie man schwellwert geschickt wählen soll, da man sqrt(a) ja noch nicht kennt. Idealerweise sqrt(a)*sqrt(epsilon) mit der zusätzlichen Iteration am Ende.

Wenn die Abbruchbedingung anders aussieht

abs(x*x-a)<=anderer_schwellwert

fällt es mir zumindest einfacher, hier einen Schwellwert anzugeben, ohne sqrt(a) kennen zu müssen; denn abs((x*x-a)/a) ist ungefähr 2*abs((x-sqrt(a))/sqrt(a)) und a ist im Gegensatz zu sqrt(a) schon bekannt. Bei der Abbildung x -> x² verdoppelt sich nämlich der relative Fehler.

SeppJ schrieb:

[...]Dafür beschränkst du dich auf maximal ca. 7 Stellen Genauigkeit von den ca. 15 möglichen bei double[...]

Nein, eigentlich nicht. Ich glaube, da gibt es ein paar Missverständnisse.

Michael E.

krümelkacker schrieb:

SeppJ schrieb:

[...]Dafür beschränkst du dich auf maximal ca. 7 Stellen Genauigkeit von den ca. 15 möglichen bei double[...]

Nein, eigentlich nicht. Ich glaube, da gibt es ein paar Missverständnisse.

Kannst du das ein wenig ausführen? Wo ist mein Fehler?

Ich denke, dass bei der Methode |x² - a| < ε mit ε ≈ 10^-14 der Fehler von x zu sqrt(a) ungefähr 10^-7 ist, während bei der anderen Methode der Fehler von x zu sqrt(a) unter 10^-14 gedrückt werden kann, da sqrt(a) zwischen x und a/x liegt.

Michael E. schrieb:

Wenn ich das richtig sehe, ist bei der Methode |x² - a| < ε mit ε ≈ 10^-14 der Fehler von x zu sqrt(a) ungefähr 10^-7, während bei der anderen Methode der Fehler von x zu sqrt(a) unter 10^-14 gedrückt werden kann.

Ah, das hatte ich garicht berücksichtigt. Sehr guter Einwand!

krümelkacker schrieb:

Ich würde Deine Variante vorziehen, eben weil floatingpoint-Division typischerweise (?) 4mal langsamer ist als Multiplikation (zumindest hatte ich das mal so beobachtet auf einem core2 Rechner).

Leider nicht. Ich bin gezwungen eine Fließzahlenklasse aus einer Numerischen Library zu benutzen, da ist Division deutlich langsamer...

krümelkacker schrieb:

Da man a/x aber wieder in der nächsten Iteration braucht, muss man das ja nicht doppelt berechnen. Hier eine weniger naive Implementierung:

Genau nach so etwas habe ich gesucht. Vielen Danke Leute, ihr seid die Besten!

krümelkacker

Michael E. schrieb:

krümelkacker schrieb:

SeppJ schrieb:

[...]Dafür beschränkst du dich auf maximal ca. 7 Stellen Genauigkeit von den ca. 15 möglichen bei double[...]

Nein, eigentlich nicht. Ich glaube, da gibt es ein paar Missverständnisse.

Kannst du das ein wenig ausführen? Wo ist mein Fehler?

Ich denke, dass bei der Methode |x² - a| < ε mit ε ≈ 10^-14 der Fehler von x zu sqrt(a) ungefähr 10^-7 ist,

Du sprichst von absoluten Fehlern, ich von relativen. Außerdem ist das falsch, was Du sagst. Wenn Du für a eine Zahl in der Nähe von 1 wählst, dann gilt sqrt(a) ≈ (1+a)/2. Der absolute Fehler halbiert sich also in diesem Fall statt sich zu ver-10⁷-fachen.

Mit ε bezeichnen wir in der Numerik typischerweise die Maschinengenauigkeit. Das ist eine obere Schranke für den relativen Fehler, wenn wir versuchen, eine Zahl als Fließkommazahl mit endlicher Mantisse darzustellen. Wenn Du |x² - a| ≤ g₁ oder |x - a/x| ≤ g₂ als Abbruchbedingung benutzen willst und für x einen relativen Fehler von ε oder weniger haben möchtest, musst Du dementsprechend g₁=2 ε a und g₂=ε a^1/2 wählen -- wenn man weitere Rechenungenauigkeiten ignoriert. Und da wir a^1/2 nicht kennen, macht IMHO die Abbruchbedingung |x² - a| ≤ g₁ mehr Sinn.

In der Praxis gibt es dann noch das Problem, dass weitere Rundungsfehler auftreten und dass die Schleife möglicherweise deswegen nicht terminiert. Deswegen hatte ich vorgeschlagen, statt ε einfach ε^1/2 zu benutzen und eine finale Iteration dranzuhängen, die nach dem "Abbruch" trotzdem ausgeführt wird. ε^1/2 deswegen, weil ich erwarte, dass die finale Iteration die Genauigkeit verdoppelt und den relativen Fehler wieder auf ε drückt.

Es kommt wirklich drauf an, was man ausrechnen will. Möchtest Du die Abbruchbedingung bzgl der absoluten Genauigkeit g2 von x ausdrücken, dann macht |x - a/x| ≤ g₂ mehr Sinn. Man kann hier aber g2 nur sinnvoll wählen, wenn man die Größenordnung von x einigermaßen kennt.

Ich bin einfach mal davon ausgegangen, dass man unabhängig von der Größenordnung von a ein x bestimmen will, welches immer mindestens soundso viele signifikante korrekte Stellen hat.

nachtfeuer

wäre es nicht irgendwie auch sinnvoller, das ganze integermäßig zu berechnen?
Man kann ja die Nachkommastellen einfach abschneiden oder bis zu einem bestimmten Punkt kleinhalten.
Ich glaube aber, bis man bis zur Genauigkeit X berechnet hat, könnte das ganze per handschriftlicher Umsetzung (zumindest quatratisch) - eventuell auf SSE - durchaus schon fertig sein.

krümelkacker

Hatte mal Lust drauf...

#include <iostream>
#include <ostream>
#include <cmath>
#include <cfloat>

const double sqrt_eps = std::sqrt(DBL_EPSILON);

double wurzel_schaetzung(double a)
{
  int exp;
  double m = std::frexp(a,&exp);
  // 0.5<=m<=1, m*pow(2,exp)=a
  if (exp&1u) { --exp; m*=2; }
  // 0.5<=m<=2, m*pow(2,exp)=a, exp%2=0
  return std::ldexp(0.5+0.5*m,exp/2);
}

double wurzel(double a)
{
  using std::abs;
  if (a==0) return 0;
  double x = wurzel_schaetzung(a);
  double s = x*sqrt_eps;
  double t = a/x;
  for (;;) {
    x = 0.5*(x+t);
    t = a/x;
    if (abs(x-t)<=s) break;
  }
  x = 0.5*(x+t);
  return x;
}

int main()
{
  using namespace std;
  cout.precision(15);
  cout << wurzel(2) << endl;
}

frexp und ldexp sind erwartungsgemäß schnell (Manipulation der rohen Fließkomma-Bits) und bieten damit einen portablen Weg, die Wurzel schnell zu schätzen. Eigentlich kann man hier auch für ein bekanntes DBL_EPSILON eine bestimmte feste Iterationsanzahl benutzen, da der relative Fehler der hier verwendeten Schätzung beschränkt ist (6%).

hustbaer

SeppJ schrieb:

Wenn man so viel Angst hat vor Division (4 mal länger sind sind immer noch bloß ca. 8 Taktzyklen) kann man ja auch am Anfang einmal 1/a ausrechnen.

Naja, ne, eher mul ~ 1-2 und div ~ 20-50 je nach CPU (Throughput).

Latenz kommt schon eher hin mit Faktor 4 (mul ~5, div ~ 20-50)