SeppJ schrieb:
Was ist, wenn ich die Wurzel von 10^-50 bestimmen will?
Um mal das Problem zu verdeutlichen habe ich ein kleines Testprogramm geschrieben. Wenn wir die Standardwurzel mal als exakt ansehen, dann vergleichen wir mal Ergebnisse. Eine absolute Präzisionsvorgabe ist schlecht.
Als Vergleich habe ich mal die berühmt berüchtigte Q3-Wurzel genommen (auch ein Newtonverfahren) mit einmal 1 Schritt und einmal 2 Schritten, Schrittzahl fest, ohne feste Präzisionsvorgabe. Man beachte die absoluten und relativen Abweichungen von der Standardwurzel der verschiedenen Verfahren für verschiedene Größenordnungen der Eingabe:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
double sqrt_knoppi(double a) {
double xneu = 0, xalt = 0, epsilon = 0.00000000001;
xneu = a;
do {
xalt = xneu;
xneu = 0.5*(xalt+a/xalt);
} while((xalt-xneu) > epsilon);
return xneu;
}
const int64_t magic_number = 0x5fe6eb50c7b537a9;
double sqrt_q3(const double x)
{
const double xhalf = 0.5f*x;
union
{
double x;
int64_t i;
} u;
u.x = x;
u.i = magic_number - (u.i >> 1);
return x*u.x*(1.5f - xhalf*u.x*u.x);
}
double sqrt_q3_2(const double x)
{
const double xhalf = 0.5f*x;
union
{
double x;
int64_t i;
} u;
u.x = x;
u.i = magic_number - (u.i >> 1);
u.x = u.x*(1.5f - xhalf*u.x*u.x); // Mehr von dieser Zeile für höhere Präzision
return x*u.x*(1.5f - xhalf*u.x*u.x);
}
void print_diff(double result, double sqrt)
{
printf("%.15f \tdiff = %.15f\tRel = %.15f\n", result, fabs(sqrt -result), fabs(sqrt -result)/sqrt);
}
void sqrt_test(double d)
{
double result = sqrt(d);
printf("sqrt(%.15f) = %.15f\n", d, result);
printf("knoppi: sqrt(%.15f) = ", d); print_diff(sqrt_knoppi(d), result);
printf("knoppi: sqrt_q3(%.15f) = ", d); print_diff(sqrt_q3(d), result);
printf("knoppi: sqrt_q3_2(%.15f) = ", d); print_diff(sqrt_q3_2(d), result);
putchar('\n');
}
int main()
{
double d = 15;
for(int i = 0; i < 30; ++i)
{
sqrt_test(d);
d /= 3;
}
}
https://ideone.com/X65xgk
