sqrt funktion nachbauen

Hi,

ich moechte die wurzel einer zahl ohne einer inbuilt funktion berechnen... ich verwende binary search...bitte um feedback.

die zeitkomplexitaet ist O(n log n) ?

double square_root(double num) {
	double l = 0;
	double r = num;
	double percision = 0.001;

	if(num < 0) {
		return -1; // error
	}

	if(num == 0 || num == 1) {
		return num;
	}

	if(num < 1.0) {
		r = 1.0;
	}

	while((r - l) > percision) {
		double mid = l + (r - l) / 2.0;
		double midsq = mid * mid;

		if(midsq == num) {
			return mid;
		}
		else if(midsq > num) {
			r = mid;
		}
		else {
			l = mid;
		}
	}

	return l;
}

meinte: O(log n) ... fuer worst case

SeppJ

Dass das an sich keine gute Methode ist, ist dir klar, oder?

precision schreibt man anders, als du denkst.

Ansonsten ganz ok, denke ich.

Zeitkomplexität: Sieht kompliziert aus, ich gebe die Frage mal ab, da ich gerade faul bin.

Th69

if(midsq == num)

Ganz schlecht, (berechnete) Fließkommazahlen direkt miteinander zu vergleichen.
Die Abfrage wird also nur ganz selten 'true' zurückliefern.

Aber sie ist auch sowieso überflüssig, da du ja anhand der Precsion die Schleife durchläufst und am Ende die Näherung zurückgibst.

Und ein paar andere Verbesserungen siehst du noch unter http://code.antonio081014.com/2013/06/square-root-from-binary-search-to.html (wenn auch Java-Code, aber das ändert ja nichts am Algorithmus)

hi, welche methode ist schneller um die wurzel zu berechnen?

SeppJ

algoman1 schrieb:

hi, welche methode ist schneller um die wurzel zu berechnen?

Erst einmal Newtons Methode probieren, damit liegt man sicherlich nicht schlecht. Das eingebaute sqrt dürfte aber schneller sein. Die benutzen zwar wahrscheinlich auch eine Variation von Newtons Methode, aber noch zusätzlich irgendwelche Mikrooptimierungen. Oder direkt Hardwareroutinen der CPU, die zwar auch ähnlich funktionieren sollten, aber natürlich nochmal viel schneller sind.

Es wird häufig auch CORDIC eingesetzt.

hi, eine moeglichkeit ist die newton raphson's method...
warum wird am anfang x0 = a/2 zugewiesen?

double sqrt(double a) {
   double x0 = a/2;
   double err = 1, eps = Math.pow(10,-11);
   while (err >= eps) {
       double x = x0 - ((x0 * x0 - a)/(2 * x0));
       err = x-x0;
       x0 = x;
   }
   return x0;
}

knoppi

Ehrlich gesagt ist das Scheißegal, weil das Teil sowieso ab ca.der dritten Nachkommastelle falsche Ergebnisse liefert:

2: 1.41666666666666674068

mit der obigen Funktion

2: 1.41421356237309514547

mit der C-Funktion sqrt

5: 2.25000000000000000000

mit der obigen Funktion

5: 2.23606797749978980505

mit der C-Funktion sqrt

Das scheinen alles, inklusive Newtons Formel, mehr oder weniger Variationen der Iterationsformel von Archimedes zu sein. Probiers doch einfach mal mit dem Original

Viele Grüße

knoppi

Bashar

knoppi schrieb:

Ehrlich gesagt ist das Scheißegal, weil das Teil sowieso ab ca.der dritten Nachkommastelle falsche Ergebnisse liefert:

Ob das an der 'percision' liegen könnte? :p

Das scheinen alles, inklusive Newtons Formel, mehr oder weniger Variationen der Iterationsformel von Archimedes zu sein.

Die Archimedes-Formel kenn ich zwar nicht, aber vielleicht meinst du das Heron-Verfahren. Das ist in der Tat ein Spezialfall des Newton-Verfahrens.

knoppi

Iterationsformel von Archimedes von Syrakus

xneu = 0.5*(xalt+a/xalt)

#include <math.h>
#include <stdio.h>

/* Errechnen der Quadradwurzel nach der Iterationsformel von Archimedes */

double archimed(double);

int main(int argc, char *argv[]) {
	double a = 0;

	printf("Von welcher Zahl soll die Quadratwurzel gezogen werden ?\n");
	scanf("%lf", &a);

	printf("\nDie Quadratwurzel von %.2lf ist %.11lf\n", a, archimed(a));
	printf("\nDie Quadratwurzel von %.2lf ist %.11lf\n", a, sqrt(a));

	return 0;
}

double archimed(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;
}

Ich setze hier mal eine Eingabe für a >= 1 voraus, um grad das Wesentliche zu zeigen.

Ich hab den Vergleich mit sqrt() mit einbezogen. Außerdem habe ich noch Vergleiche mit Taschenrechner und dem hier auf meinem PC gemacht. Ist zwar nicht unbedingt ein Beweis, aber bisher stimmte alles.

Viele Grüße

knoppi

SeppJ

Was ist, wenn ich die Wurzel von 10^-50 bestimmen will?

Bashar

knoppi schrieb:

Iterationsformel von Archimedes von Syrakus
xneu = 0.5*(xalt+a/xalt)

Auch bekannt als Heron-Verfahren. Hatte ich glaub ich schon erwähnt... manche nennen es auch babylonisches Verfahren.

SeppJ

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