Container verschnellern....

[It0101]

volkard schrieb:

und wechle mal nach uint64_t

Ist Teil des C-Standards, daher nicht in C++ zu verwenden. *pfeif*

[knivil]

It0101 schrieb:

Ist Teil des C-Standards, daher nicht in C++ zu verwenden. *pfeif*

Pfeife: http://en.cppreference.com/w/cpp/types/integer

[It0101]

Gut. Dann sind die inzwischen dabei. Ich nehme alles zurück und behaupte das Gegenteil

[volkard]

knivil schrieb:

Aber 5896847556 ist auch kein guter Testcase, da nur Primzahlen bis 80000 gebraucht werden, also etwa 8000 Stueck.

Vielleicht magste die eher:
2^61-1: vieele Teiler
2^61-7: mehrfacher Primfaktor und große Primfaktoren
2^61-87: drei große Primfaktoren

[Columbo]

Ich glaube, ich weiß, was du meinst: Du checkst die "Standard-Teiler" ab, sprich, alle unter ~100000 (ist das ein gutes Maß?). Solange die Zahl dadurch nicht komplett zerlegt wird, zerlegt man sie in zwei kleine Teiler, die man wiederum in ihre Primfaktoren zerlegt.

Meine beginnt so:

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

Für x86:

// 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;
}

Die Approximation kann man natürlich schlechter machen.

Die Approximation kann auch dem GCC überlassen werden, der hat eine schöne Funktion, genannt** __builtin_clzll **.

[volkard]

Arcoth schrieb:

Ich glaube, ich weiß, was du meinst: Du checkst die "Standard-Teiler" ab, sprich, alle unter ~100000 (ist das ein gutes Maß?).

Würde mal damit rechnen, daß vielleicht auch 256 gut ist. Lassen wir es offen und messen es später aus.

Arcoth schrieb:

Solange die Zahl dadurch nicht komplett zerlegt wird, zerlegt man sie in zwei kleine Teiler, die man wiederum in ihre Primfaktoren zerlegt.

Jo.

[volkard]

Arcoth schrieb:

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;
}

Hoffe, uint64 ist gemeint.
Und innendrin nur double??? Dann bleibste ja ungenau, oder?

Arcoth schrieb:

Die Approximation kann man natürlich schlechter machen.

Und alles andere. Da geht noch einiges. Von knivils Trick gehen wir eh aus, hoffe ich.

Arcoth schrieb:

Die Approximation kann auch dem GCC überlassen werden, der hat eine schöne Funktion, genannt** __builtin_clzll **.

Die rechnet Dir nur log2(n) blitzschnell aus.

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