Wurzel für integer berechnen.

*john 0

Ich habe es jetzt mit folgendem Code aus dem englischen Wikipedia Artikel zum Thema etwas in Relation zu sqrt verbessern können. Der Hauptvorteil ist, dass der Code mit beliebig breiten Integer funktioniert.

template<typename T>
requires std::unsigned_integral<T>
inline constexpr T int_sqrt (const T s) noexcept {
	if (s <= 1) return s;

	// Initial estimate (must be too high)
	T x0 = (1 << ((std::bit_width(s) / 2) + 1));
	T x1 = (x0 + s / x0) / 2; 

	// Bound check
	while (x1 < x0) {
		x0 = x1;
		x1 = (x0 + s / x0) / 2;
	}

	return x0;
}

Quiche Lorraine

@john-0 sagte in Wurzel für integer berechnen.:

// Bound check
while (x1 < x0) {
x0 = x1;
x1 = (x0 + s / x0) / 2;
}

Kommt mir nicht unbekannt vor

Ich hatte mir hierzu vor längerem ein paar Gedanken gemacht, wieviele Schritte ich für eine "gute" Näherung benötige. Und hierbei bin ich auf ln(n)/ln(2) + m Schritte gekommen.

Hierzu habe ich einen möglichst guten Startwert mittels binärer Suche gesucht. Nehmen wir mal an, ich suche die Wurzel von 10.

Die Wurzel muss im Intervall [1; 10] liegen, da 1*1 < 10 und 10 * 10 > 10 ist.
Nun betrachte ich die Mitte (1 + 10) / 2 = 5.5. Und 5.5 * 5.5 = 30.25, ist also größer als 10. Die Wurzel muss also zwischen 1 und 5.5 liegen.
Nun betrachte ich die Mitte (1 + 5.5) / 2 = 3.25. Und 3.25 * 3.25 = 10.5625, ist also größer als 10. Die Wurzel muss also zwischen 1 und 3.25 liegen.
Nun betrachte ich die Mitte (1 + 3.25) / 2 = 2.125. Und 2.125 * 2.125 = 4.515, ist also kleiner als 10. Die Wurzel muss also zwischen 2.125 und 3.25 liegen.

Ich habe bei einer Intervall-Größe kleiner 1 abgebrochen und habe die Intervall-Mitte als Startwert benutzt. Und danach zeigte sich nach ein paar Schritten die ersten stabilen Werte, wo also x(n+1) - x(n) < 1 ist und pro Schritt eine fixe Stelle hinzukam.

Ist jetzt zwar etwas am Thema vorbei, dachte aber die binäre Suche könnte hier auch helfen.

*john 0

Dann hatte ich noch diesen Code im Netz gefunden, der ohne Division auskommt. Allerdings bringt das keinen substantiellen Vorteil. An der Art der Funktionsdefinition sieht man schon wie alt der Code ist.

/*
 *	Square root by abacus algorithm, Martin Guy @ UKC, June 1985.
 *	From a book on programming abaci by Mr C. Woo.
 *	Argument is a positive integer, as is result.
 *
 *	I have formally proved that on exit:
 *		   2		   2		   2
 *		res  <= x < (res+1)	and	res  + op == x
 *
 *	This is also nine times faster than the library routine (-lm).
 */

int
sqrt(x)
int x;
{
	/*
	 *	Logically, these are unsigned. We need the sign bit to test
	 *	whether (op - res - one) underflowed.
	 */
	
	register int op, res, one;

	op = x;
	res = 0;

	/* "one" starts at the highest power of four <= than the argument. */

	one = 1 << 30;	/* second-to-top bit set */
	while (one > op) one >>= 2;

	while (one != 0) {
		if (op >= res + one) {
			op = op - (res + one);
			res = res +  2 * one;
		}
		res /= 2;
		one /= 4;
	}
	return(res);
}

Hallo @john-0 ,

das ist mit dem Intervall-Halbierungs-Verfahren (nach Isaac Newton) doch wirklich überhaupt keine Herausforderung mehr...

Beispiel:

#include <cstdint>
#include <iostream>

double wurzel(int32_t i32) {
  i32 = abs(i32);
  double a = 0;
  double b = i32;
  double c = 0;
  double d = 0;
  double mid = 0;
  do {
    c = d;
    mid = (a + b) / 2.0;
    d = mid * mid;
    if (d == i32) {
      return mid;
    }
    if (d < i32) {
      a = mid;
    } else {
      b = mid;
    }
  } while (c != d);
  return mid;
}

int main(int argc, char *argv[]) {
  for (int32_t i = -20; i < 20; i++) {
    std::cout << i << "  " << wurzel(i) << std::endl;
  }
  return 0;
}

Klar, es mag bessere und effizientere und genauere Verfahren geben, aber du wolltest glaube ich einen einfachen Algo...

Sorry, ich hatte jetzt nicht daran gedacht, dass es für negative Zahlen ja gar keine Wurzel gibt ... mathematisch also eher stümperhaft.

*john 0

@Fragender sagte in Wurzel für integer berechnen.:

Sorry, ich hatte jetzt nicht daran gedacht, dass es für negative Zahlen ja gar keine Wurzel gibt ... mathematisch also eher stümperhaft.

Es gibt für negative Zahlen natürlich eine Wurzel, aber die ist dann imaginär, und dann passt der Typ für den Rückgabewert nicht.

Zum eigentlichen Aspekt: Mir geht es darum, die Wurzel für Integerzahlen zu berechnen. Da ist die Intervallhalbierung zu Anfangs relativ ineffizient und das kommt immer mehr zum Tragen je größer die Zahlen werden. Dann ist der Weg über den Datentyp double, um die Wurzel für ein int32_t zu berechnen, nicht das goldene vom Ei. Es funktionier zwar, weil die Mantisse von 53Bit von der üblichen double Implementation größer ist als die 31Bit von int32_t. Nur bei größeren Ints führt das zu einem Problem. Dann ist die Division relation langsam, so dass ein Algorithmus schneller läuft, wenn man die Division eliminieren kann z.B. durch Bitshifts (deshalb die beiden Postings, der Compiler kann das relativ gut substituieren) oder durch Restklassenringmultiplikationen. Denn die sind deutlich schneller.

@Fragender sagte in Wurzel für integer berechnen.:

es mag bessere und effizientere und genauere Verfahren geben, aber du wolltest glaube ich einen einfachen Algo...

Du hast das schon gelesen, oder? @john-0

Ich nehme mal an, für deinen Anwendungsfall ist ein geringfügiger Laufzeitunterschied völlig egal und du kennst deine Anforderungen einfach nicht, oder?

Was ist denn die höchste Zahl, dessen Wurzel du berechnen willst? Und danke für die "Belehrung", die ich nicht brauche.

Leon0402

@Fragender sagte in Wurzel für integer berechnen.:

@Fragender sagte in Wurzel für integer berechnen.:

es mag bessere und effizientere und genauere Verfahren geben, aber du wolltest glaube ich einen einfachen Algo...

Du hast das schon gelesen, oder? @john-0

Ich nehme mal an, für deinen Anwendungsfall ist ein geringfügiger Laufzeitunterschied völlig egal und du kennst deine Anforderungen einfach nicht, oder?

Was ist denn die höchste Zahl, dessen Wurzel du berechnen willst? Und danke für die "Belehrung", die ich nicht brauche.

John hat doch in seinem Eingangspost geschrieben, dass er ein sehr schnelles Verfahren braucht. Ich weiß nicht, was sein Anwendungsfall konkret ist, aber basierend auf meiner Erfahrung hier im Forum, wird er das nicht zum Spaß geschrieben haben
Jemanden mit viel Erfahrung jedenfalls vorzuwerfen, dass er seine Anforderungen nicht kennt, halte ich generell für gewagt. Die Posts in diesem Thread haben auch nicht den Eindruck erweckt bisher.

Quiche Lorraine

@Leon0402 sagte in Wurzel für integer berechnen.:

Ich weiß nicht, was sein Anwendungsfall konkret ist

Ich bin hier einfach mal von Primzahlen für RSA, also Primzahlen > 2048 Bit ausgegangen.

Ist log_2(n) bereits relativ schnell,
sind die bisher gemachten Vorschläge für Zahlen > 2048 Bit ohnehin Murks.