// Programm: math.C
// Kleine Bibliothek mathematischer Funktionen.
#include
#include
#include
namespace ifm {
double sqrt(double n)
{
// PRE: n >= 0.
// POST: Gibt Approximation der Quadratwurzel
// von n zurueck.
assert(n >= 0);
if (n == 0) return 0;
// Newton-Iteration: x_{i+1} = 1/2 (x_i + n/x_i) Wir beginnen mit
// dem Startwert x_0:=(n+1)/2. Da (n+1)/2-sqrt(n) =
// (n-2sqrt(n)+1)/2 = (sqrt(n)-1)^2/2 > 0, gilt x_0 > sqrt(n).
//
// Wie in der Theorie verlangen wir, dass der Approximationswert
// in jedem Schritt kleiner wird. Nur muss sich hier die Aenderung
// auch innerhalb der double Mantisse auswirken.
//
// Damit terminiert das Programm, und wir koennen sogar beweisen,
// dass "ungefaehr" das richtige herauskommt
double curr = (n + 1) / 2; // Startwert
double prev; // voriger Wert
const double roundoff =
curr * std::numeric_limits::epsilon();
do {
prev = curr;
assert(curr > 0);
curr = (curr + n / curr) / 2;
} while (prev - curr > roundoff);
return curr;
} // double sqrt(double)
} // namespace ifm
