Container verschnellern....

Columbo

Mit Approximation meinte ich ja auch den Seed für den Algo, nicht die eigentliche Wurzel.

Hoffe, uint64 ist gemeint.

int64 ist ein Typedef auf std::uint64_t ; ist ungünstig, aber so muss ich ein u weniger schreiben (und das es vorzeichenlos ist, ist irgendwie evident).

Und innendrin nur double??? Dann bleibste ja ungenau, oder?

Wie das? Die Methode konvergiert quadratisch gegen die Wurzel. der Startwert ist nie höher als $2^{32}$ (s. Kommentar), wobei double genau ist bis $2^{53}$ (zumindest bei IEEE).

Ich konnte heute nicht wirklich irgendwas vernünftiges Implementieren, das gute Wetter - ich fühlte mich verpflichtet, den Altweibersommer zu genießen.

volkard

Arcoth schrieb:

Und innendrin nur double??? Dann bleibste ja ungenau, oder?

Wie das? Die Methode konvergiert quadratisch gegen die Wurzel. der Startwert ist nie höher als $2^{32}$ (s. Kommentar), wobei double genau ist bis $2^{53}$ (zumindest bei IEEE).

Ja, aber die Rechnung x = (x + n/x)/2; geschieht dann auf double? also x = (x + n/x)/2;? Dann verlöre n leider 11 Bits. Oder macht das aus mathematischen Gründen da nix aus?

camper

volkard schrieb:

Ja, aber die Rechnung x = (x + n/x)/2; geschieht dann auf double? also x = (x + n/x)/2;? Dann verlöre n leider 11 Bits. Oder macht das aus mathematischen Gründen da nix aus?

Ich bin fast sicher, dass das nichts ausmacht. Vielleicht mal ins Matheforum stellen.

camper

... oder auch nicht
Bruteforce-Suche, mit der Annahme, dass sich isqrt monoton verhält, dann reicht es, nur die Umgebung von Quadratzahlen zu untersuchen.

#include <cstdlib>
#include <cmath>
#include <iostream>

using namespace std;

uint64_t isqrt( uint64_t const n )
{
    double x = n >> (uint64_t)log2(n)/2; // x ist höchstens 2^32
    double old_x;
    for(;;)
    {
        old_x = x;
        x = (x + n/x)/2;

        if( std::abs(x - old_x) < 1 )
            break;
    }

    return x;
}

int main()
{
    for ( uint64_t i = 2; i < (1ull<<32); ++i )
    {
        if ( ( i & ( i - 1 ) ) == 0 )
            cout << i << '\n';
        uint64_t square = i * i;
        if ( isqrt( square - 1 ) != i - 1 )
            cout << "error: " << i << " ^2 - 1\n";
        if ( isqrt( square ) != i )
            cout << "error: " << i << " ^2\n";
        if ( isqrt( square + 1 ) != i )
            cout << "error: " << i << " ^2 + 1\n";
    }
}

findet schnell Fehler, z.B. ideone (auf meinem Rechner bekomme ich andere Zahlen). Ganz nebenbei terminiert isqrt nicht für n==1

volkard

Arcoth schrieb:

Für x86:

Ja, ruhig (in Maßen, um z.B. den einen guten Befehl zu holen) inline-asm und 64-bitter voraussetzen.

Arcoth schrieb:

// Böse geklaut von http://stackoverflow.com/questions/994593/how-to-do-an-integer-log2-in-c
int64 log2(int64 x)
{
	int64 y;
	asm ( "\tbsr %1, %0\n"
		: "=r"(y)
		: "r" (x) );
	return y;
}

int64 isqrt( int64 const n )
{
	double x = n >> log2(n)/2; // x ist höchstens 2^32
	double old_x;
	for(;;)
	{
		old_x = x;
		x = (x + n/x)/2;

		if( std::abs(x - old_x) < 1 )
			break;
	}

	return x;
}

Mein Testprogramm beschwert sich

isqrt(3891049815758399)!=62378279

volkard

camper schrieb:

... oder auch nicht
Bruteforce-Suche, mit der Annahme, dass sich isqrt monoton verhält, dann reicht es, nur die Umgebung von Quadratzahlen zu untersuchen.

#include <cstdlib>
#include <cmath>
#include <iostream>

using namespace std;

uint64_t isqrt( uint64_t const n )
{
    double x = n >> (uint64_t)log2(n)/2; // x ist höchstens 2^32
    double old_x;
    for(;;)
    {
        old_x = x;
        x = (x + n/x)/2;

        if( std::abs(x - old_x) < 1 )
            break;
    }

    return x;
}

int main()
{
    for ( uint64_t i = 2; i < (1ull<<32); ++i )
    {
        if ( ( i & ( i - 1 ) ) == 0 )
            cout << i << '\n';
        uint64_t square = i * i;
        if ( isqrt( square - 1 ) != i - 1 )
            cout << "error: " << i << " ^2 - 1\n";
        if ( isqrt( square ) != i )
            cout << "error: " << i << " ^2\n";
        if ( isqrt( square + 1 ) != i )
            cout << "error: " << i << " ^2 + 1\n";
    }
}

findet schnell Fehler, z.B. ideone (auf meinem Rechner bekomme ich andere Zahlen). Ganz nebenbei terminiert isqrt nicht für n==1

Hast aber geschummelt. Sein log2 packt mehr.
http://ideone.com/Zxt5at

camper

volkard schrieb:

Hast aber geschummelt. Sein log2 packt mehr.
http://ideone.com/Zxt5at

ideone schrieb:

Time limit exceeded

volkard

camper schrieb:

volkard schrieb:

Hast aber geschummelt. Sein log2 packt mehr.
http://ideone.com/Zxt5at

ideone schrieb:

Time limit exceeded

Ja, weil er in den ersten paar Sekunden keinen Fehler findet und nicht in der Zeit selber terminiert.
Upps, der 0-Fehler ist ja noch drin. Moment.

Auch mit
if(n<100) return sqrt(double(n));
time exceeded. muss man wohl lokal testen.

ok, alles repariert. trotzdem fehler bei 557580^2-1

Mit if( std::abs(x - old_x) < 0.5 ) kommt er bis um 6000000^2, aber hübsch ist anders.

camper

Die Genauigkeit des Startwertes sollte ohnehin irrelevant sein.
x=1 ist genauso gut oder schlecht (nur ein bisschen langsamer).

Th69

Bei euch geht es um die ganzzahlige Wurzel, oder?

Diesen Code hatte ich mal vor ca. 10 Jahren geschrieben:

template <class T> T sqrt(T x)
{
	if(x <= 1) return x;

	const T X = 1 << (std::numeric_limits<T>::digits >> 1);
	const T X1 = 1 << (std::numeric_limits<T>::digits >> 2);
	T x1 = (x < X)? 1 : X1;
	T x2 = (x < X)? X1 : X;
	T n;

	while(n = (x1 + x2)/2, x2 > x1 + 1)
	{
		if(n * n > x)
			x2 = n;
		else
			x1 = n;
	}

	return n;
}

Kommt ohne Umwandlung nach double aus. Was haltet ihr davon?

volkard

Th69 schrieb:

Bei euch geht es um die ganzzahlige Wurzel, oder?

Diesen Code hatte ich mal vor ca. 10 Jahren geschrieben:
template <class T> T sqrt(T x)
{
	if(x <= 1) return x;

	const T X = 1 << (std::numeric_limits<T>::digits >> 1);
	const T X1 = 1 << (std::numeric_limits<T>::digits >> 2);
	T x1 = (x < X)? 1 : X1;
	T x2 = (x < X)? X1 : X;
	T n;

	while(n = (x1 + x2)/2, x2 > x1 + 1)
	{
		if(n * n > x)
			x2 = n;
		else
			x1 = n;
	}

	return n;
}
Kommt ohne Umwandlung nach double aus. Was haltet ihr davon?

Kleine Reparatur war nötig.

const T X = T(1) << (std::numeric_limits<T>::digits >> 1);
    const T X1 = T(1) << (std::numeric_limits<T>::digits >> 2);

Ich lasse meinen Tester laufen…
Ist korrekt bis obenhin.
Ist 6-mal langsamer als meine Implementierung.

Columbo

Ist 6-mal langsamer als meine Implementierung.

Was ist denn deine Implementierung? Oder ist die streng geheim?

Mit if( std::abs(x - old_x) < 0.5 ) kommt er bis um 6000000^2, aber hübsch ist anders.

(Vor allem, weil das den Prozess verlangsamt)

Aber Moment. Die Variante, die ich implementiert habe, konvergiert laut Wikipedia quadratisch - und die Bedingung ist immer richtig. Das Problem ist eher, dass double nicht gut genug ist.

Schaue ich mir morgen früh genau an.

Und kann das mal jemand splitten? Oder gehört das hier rein?

volkard

Arcoth schrieb:

Ist 6-mal langsamer als meine Implementierung.

Was ist denn deine Implementierung? Oder ist die streng geheim?

'Türlich ist sie nicht geheim. Und sie ist bestimmt auch noch suboptimal. Ich wollte sie nur genau Dir nicht gleich verraten, damit Du Dich ein wenig mehr in Algos/Mathe reinstupst. Ich geifere nach den Verbesserungen, die Ihr klar sehen werdet, wo ich betriebsblind bin.
6-mal ist natürlich kein Maß. 6-mal bei campers Test über den geamten Zahlenbereich.

volkard

Arcoth schrieb:

Und kann das mal jemand splitten? Oder gehört das hier rein?

Den Thread haben wir erfolgreich entführt. Ist auch kein Problem, der TO ist hier nicht mehr. Allenfalls könnte eine Mod den Threadtitel ändern.

Columbo

#include <iostream>
#include <cmath>

using uint64 = std::uint64_t;

// Zu viel Generizität? Bin mir nicht sicher, du predigst einen mMn. ungesunden Minimalismus.
template<typename T>
T intlog2(T x)
{
	asm ( "bsr %1, %0"
		: "=r"(x)
		: "r" (x) );

	return x;
}

/// Böse geklaut (aber lieb modifiziert)
// von http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
uint64 isqrt(uint64 num)
{
	uint64 bit = 1ull << (intlog2(num) ^ 1);

	uint64 res = 0;
	while( bit )
	{
		if( num >= res + bit )
		{
			num -= res + bit;
			res = (res >> 1) + bit;
		}
		else
			res >>= 1;

		bit >>= 2;
	}
	return res;
}

int main()
{
	std::cout << isqrt(3891049815758399);
}

Leider verstehe ich das Verfahren an sich noch nicht, das muss ich gleich noch bei Wikipedia genau nachlesen. (Schriftliches Wurzelziehen, auau).

Allenfalls könnte eine Mod den Threadtitel ändern.

Hätten wir bloß einen... aber die sind ja bekanntlich nie da wenn man sie braucht

knivil

Ich hoffe, so steht es nicht wirklich in deinem Code. Da ist nämlich UB drin, solange int 32-Bit groß ist.

Dieser "Fehler" ist leicht zu beheben.

volkard

Arcoth schrieb:

// Zu viel Generizität? Bin mir nicht sicher, du predigst einen mMn. ungesunden Minimalismus.
template<typename T>
T intlog2(T x)
{
    asm ( "bsr %1, %0"
        : "=r"(x)
        : "r" (x) );

    return x;
}

Den Minimalismus könnte man auch Feigheit nennen.
Warum Fehler riskieren? Du tust hier für 2 Funktionen intlog2(uint32) und intlog2(uint64) das Templatefaß aufmachen. Danach überlegste, ob es überhaupt schlau ist, x wiederzuverwenden und machst es evtl bei intlog2(uint32) aber bei intlog2(uint64) nicht.
Naja, kannst hier ruhig mal die Sache generisch machen, der Code ist ja erstmal völlig gleich. Und es kann ja zu gar keinem Fehler kommen. Richtig?

Nicht ganz.

double x=16;
	cout<<intlog2(x)<<'\n';//ausgabe 3.06321e-322//kacke

	complex<int> x=16;
	cout<<intlog2(x)<<'\n';//ausgabe 4//ok

	complex<int> x(0,16);
	cout<<intlog2(x)<<'\n';//ausgabe 36//im ernst jetzt?

	char x='a';
	cout<<intlog2(x)<<'\n';//compilerfehler operand type mismatch
//aha, wird die ausgabe sein, für kleinere typen brqauchts 
//anderen code ohne x wiederzuverwenden. 

	char const* x="evil";
	cout<<intlog2(x)<<'\n';//absturz

Das war generischer Unfug.

Columbo

Teste doch wenigstens mal rudimentär bevor Du sowas postest. Ist voll ärgerlich, daß ich Deinen Code nie messen kann, weil er falsch ist.

Wat!?

Stimmt, die Funktion gibt teilweise Blödsinnige Werte (von 18 - 31 ist die "Wurzel" 6). Entweder habe ich den Algo ruiniert, oder die Implementierung ist einfach falsch. In letzterem Fall darf ich mich wohl nicht auf Wikipedias Coder verlassen.
Versuche ich mal zu fixen.

Edit: Nein, ich bin Schuld.

Columbo

Ja, der Fehler ist schnell gefunden, und peinlich. Ich habe natürlich das bescheuerte Bit-Spielchen mit dem XOR falsch gemacht. So ist die Zeile zu korrigieren:

uint64 bit = 1ull << (intlog2(num) & ~1);

Das heißt eigentlich log- log%2 , also das letzte Bit muss auf 0 gesetzt werden, es wird auf 2 runter-gerundet. Habe wohl mal wieder schlecht geschlafen.

Edit: Bezeichner korrigiert.

knivil

template<typename T> 
 T intlog2(T x) 
{ 
     asm ( "bsr %1, %0" 
         : "=r"(x) 
         : "r" (x) ); 

     return x; 
}

Aus dem Bauch heraus: Templates und Assembler nie mischen. Templates sind fuer Metaprogrammierung durch Typen. Assembler kennt keine Typen. C++ kennt keinen Assembler. Der obige Code laeuft nur mit 'nem Integer der in ein Register passt. Warum also Template?