Wurzel für integer berechnen.

*john 0

Kennt jemand ein brauchbares Verfahren für integer Zahlen die Wurzel brauchbar exakt zu berechnen? Auf das Problem bin ich gestoßen, weil ich Programm geschrieben hat, welches Primzahlen berechnet. Ok, das ist trivial. Dabei teste ich nun jede neue Zahl mit den bereits berechneten Primzahlen. Der Punkt ist, dass man das nur bis zur Wurzel der zu testenden Zahl tun muss, dann weiß man dass das eine neue Primzahl ist. Die Version die Primzahl zu quadrieren und zu testen, ob diese größer der zutestenden Zahl ist, funktioniert. Allerdings berechnet man so für jede Primzahl ein Produkt. Für „kleine“ Zahlen kann man sqrt nehmen, das funktioniert auch ganz gut, und ist schneller. Das Problem ist dass die Mantisse begrenzt ist, und somit das nur für „kleine“ Zahlen funktioniert.

Die Abschätzung die Bitgröße der Zahl zu halbieren, ggf. ein Bit aufzurunden und alle Bit zu setzen funktioniert als Schranke auch, ist aber zu langsam, da damit zu viele Zahlen getestet werden. D.h. es bräuchte eine bessere Approximation der Wurzel. Kennt hier jemand eine bessere Lösung?

Quiche Lorraine

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

Der Punkt ist, dass man das nur bis zur Wurzel der zu testenden Zahl tun muss, dann weiß man dass das eine neue Primzahl ist.

Hast du dir schon einmal Probabilistische Primzahltests angeschaut?

https://de.wikipedia.org/wiki/Primzahltest#Probabilistische_Primzahltests
https://de.wikipedia.org/wiki/Miller-Rabin-Test

*john 0

Ich weiß, dass es da einige weitere Algorithmen gibt, das ist noch Todo. Mir ging es darum den vorhandenen Algorithmus erstmal so schnell wie möglich zu machen. Da erwies sich bisher der Ersatz der Divisionen durch Restklassenring Multiplikationen als sehr gewinnbringend.

Danke für den Lesehinweis, kommt auf die Todo-Liste, auch wenn es mir um 100% Berechnung ging und nicht um eine Abschätzung, ob eine Zahl prim ist oder nicht.

Quiche Lorraine

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

eine Abschätzung, ob eine Zahl prim ist oder nicht.

Wenn man aber den Miller-Rabin-Test wiederholt, so sinkt die Fehlerwahrscheinlichkeit exponentiell mit 4**(-n).

Alternativ kenne ich noch eine rekursive Formel für die Wurzelberechnung. Aber diese besteht pro Iteration aus einer Multipikation, einer Addition und einer Division. Dürfte also hier nicht passen.

*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.