Wie kann man sqrt am schnellsten berechnen?

khalderon

Ich brauche sqrt nur für reelle Zahlen.

pale dog

Zdravko schrieb:

Ich brauche sqrt nur für reelle Zahlen.

reichen auch ganze zahlen?

lolz

Reicht auch die Inverse? Die geht sehr schnell http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf

float sqrt (float f)
{
    if (f == 0.0f)
        return 0.0f;

    float f1 = f;
    unsigned long *ptr = (unsigned long*)&f;
    *ptr = (0xbe6f0000-*ptr)>>1;

    return f1*f;
}

TomasRiker

Das ist aber extrem ungenau ...

khalderon

pale dog schrieb:

Zdravko schrieb:

Ich brauche sqrt nur für reelle Zahlen.

reichen auch ganze zahlen?

Nein. Also z.B. haben wir float Variablen.
Ich habe bisher nur diese ziemlich gute Variante von sqrt:

float mysqrtopt(float r)
{
register float xn=r;
xn=(xn+r/xn)/2;
xn=(xn+r/xn)/2;
xn=(xn+r/xn)/2;
xn=(xn+r/xn)/2;
xn=(xn+r/xn)/2;
xn=(xn+r/xn)/2;
xn=(xn+r/xn)/2;
return xn;
}

Ob das aber wirklich schnell ist, kann ich nicht bestätigen. Hauptsache ist, dass es gut sqrt approximiert.

rapso

also den source oben koenntest du als fangsnaehrung nehmen (statt xn=r also xn=...) und statt /2 waere *0.5f eventuell gut.

pale dog

float fake_sqrt (float f)
{
	int i = (((*(int*)&f & 0x7FFFFFFF) - 0x3f800000)>>1) + 0x3f800000;
	return *(float*)&i;
}

das hätte man doch auch selber mit google finden können.

pale dog schrieb:

float fake_sqrt (float f)
{
	int i = (((*(int*)&f & 0x7FFFFFFF) - 0x3f800000)>>1) + 0x3f800000;
	return *(float*)&i;
}

Nicht gerade genau.

@Zdravko:
Nimm doch sqrt aus math.h, das benutzt das Heronverfahren.

khalderon

Und endlich? Was soll ich benutzen? Jeder sagt etwas anderes.

proggingmania

float mysqrtopt(float r)

Das wird bei etwas grösseren Zahlen, hier ab ca. r > 3000 schnell wackelig nicht wahr.

Gucke mal dies:

#define MAGIC_SQRT_NUMBER 0x5f3759df
float magic_sqrt(float x)
{
	float h = x/2;
	int i = MAGIC_SQRT_NUMBER-(*(int*)&x>>1);
	x = *(float*)&i;
	return 1/(x=x*(1.5f-h*x*x));
}

Für x < 10⁴ ganz gut.
Für grössere x kannst du vor der Ausgabe durch Einfügen der Zeilen
x=x*(1.5f-h*x*x) die Genauigkeit steigern.

Für sehr grosse Zahlen kannst du gleich sqrt aus math.h nehmen, da gibts keinen Gewinn an Geschwindigkeit mehr.

Daniel E.

Zdravko schrieb:

Und endlich? Was soll ich benutzen? Jeder sagt etwas anderes.

Welcome to the real world.

Hinsetzen, Genauigkeitsrahmen der einzelnen Lösungen prüfen, Laufzeiten nachrechnen, vergleichen, entscheiden. Viel Spaß.

Benutze dein Hirn!

proggingmania

Du könntest auch deine Version benutzen und ein bisschen verbessern, indem du deine Näherungsschritte in eine Schleife packst.
Die Schleife wird abgebrochen, wenn eine gewünschte Genauigkeit ( als PRECISION definiert ) erreicht ist.
Bei sehr grossen Zahlen und einem zu kleinen Precisionswert gibt es eine Endlosschleife.
1E-7 scheint mir ein guter Kompromiss zu sein.
Die Variable iteration_step_counter zählt die Schleifendurchläufe, die kannst du natürlich später rausschmeissen. ( falls du diese Version benutzen willst )
Bei FLT_MAX sind es bei mir 68 Durchläufe.
Je kleiner die Zahl umso weniger Durchläufe. ( Zahl=9 Durchläufe = 5 usw. )

Ob das schneller als die sqrt-Variante ist, keine Ahnung, das auszuprobieren überlasse ich dir.

#define PRECISION 1E-7
#define FLT_MAX         3.402823466e+38F   

int iteration_step_counter = 0;

float mysqrtopt_opt( float r )
{
	float x = r;
	do
	{
		x  = ( x + r / x ) / 2;
		iteration_step_counter++;

	}while( ( x*x / r ) > (1.0 + PRECISION) );

	return x;
}

int main()
{
	float z = FLT_MAX;
	float x = mysqrtopt_opt(z);
	printf("%G Anzahl der Schleifendurchlaeufe: %d\n", x, iteration_step_counter );
	return 0;
}

khalderon

proggingmania schrieb:

float mysqrtopt(float r)
Das wird bei etwas grösseren Zahlen, hier ab ca. r > 3000 schnell wackelig nicht wahr.

Gucke mal dies:
#define MAGIC_SQRT_NUMBER 0x5f3759df
float magic_sqrt(float x)
{
	float h = x/2;
	int i = MAGIC_SQRT_NUMBER-(*(int*)&x>>1);
	x = *(float*)&i;
	return 1/(x=x*(1.5f-h*x*x));
}
Für x < 10⁴ ganz gut.
Für grössere x kannst du vor der Ausgabe durch Einfügen der Zeilen
x=x*(1.5f-h*x*x) die Genauigkeit steigern.

Für sehr grosse Zahlen kannst du gleich sqrt aus math.h nehmen, da gibts keinen Gewinn an Geschwindigkeit mehr.

Ich werde das benutzen. Vielen Dank!

Mr. N

float sqrtf(float);
double sqrt(double);

Die beiden Funktionen sollten verdammt schnell sein, zumindest schneller als irgendwelche Bithacks. (Wozu gibts die FPU?)

khalderon

Mr. N schrieb:

float sqrtf(float);
double sqrt(double);
Die beiden Funktionen sollten verdammt schnell sein, zumindest schneller als irgendwelche Bithacks. (Wozu gibts die FPU?)

Weißt du wie schnell die FPU von AT91 oder TI MSP ist?

Mr. N

Zdravko schrieb:

Mr. N schrieb:
float sqrtf(float);
double sqrt(double);
Die beiden Funktionen sollten verdammt schnell sein, zumindest schneller als irgendwelche Bithacks. (Wozu gibts die FPU?)
Weißt du wie schnell die FPU von AT91 oder TI MSP ist?

Nö. Woher soll ich das wissen? Haben die überhaupt eine? Wenn nicht, wird das Standard-sqrt ja wohl per Bibliothek sein, und vermutlich wird die Bibliothek sogar recht gut sein.

Aber du kannst ja Benchmarks machen und vergleichen.

pale dog

Mr. N schrieb:

Zdravko schrieb:
Mr. N schrieb:
float sqrtf(float);
double sqrt(double);
Die beiden Funktionen sollten verdammt schnell sein, zumindest schneller als irgendwelche Bithacks. (Wozu gibts die FPU?)
Weißt du wie schnell die FPU von AT91 oder TI MSP ist?
Nö. Woher soll ich das wissen? Haben die überhaupt eine?

haben sie nicht.
ich würde dem op ja zu festkommaarithmetik raten, mit floats tut er sich keinen gefallen.
bei AT91SAM7's isses noch machbar, die haben genug power (ARM7), aber 'nem MSP430 würde ich das nicht zumuten wollen. der ist selbst mit C++ überfordert. assembler oder C schmecken ihm am besten