Container verschnellern....

volkard

Manchmal verwende ich Templates, um die automatischen Casts loszuwerden.
Wer intlog2 mit was anderem als einem unsigned auruft, denkt, die Funktion würde was anderes machen, als sie macht. Daher

template<typename T>
int intlog2(T)=delete;

template<>
int intlog2<uint32>(uint32 x)
{
	uint32 r;
    asm ( "bsrl %1, %0"
        : "=r"(r)
        : "r" (x) );
    return r;
}

template<>
int intlog2<uint64>(uint64 x)
{
	uint64 r;
    asm ( "bsrq %1, %0"
        : "=r"(r)
        : "r" (x) );
    return r;
}

volkard

Arcoth schrieb:

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.

Damit klappt Dein Code.
886 Sekunden.
Th69 hatte 1000.
Ich habe 240.
Soll ich ihn zeigen oder willste noch ein wenig basteln?

volkard

Newton ganzzahlig. 860s.

UInt64 sqrt(UInt64 x){
	if(x==0) return 0;
	UInt64 sn=x>>(intlog2(x)/2);
	UInt64 so;
	do{
		so=sn;
		sn=(so+x/so)/2;
	}while(sn<so);
	if(sn*sn>x) --sn;
	return sn;
}

knivil

Nur so als Zwischenfrage da ich auch mal ueber optimales isqrt bei Primzahlsieb nachgedacht habe: Wie gross ist der Anteil von isqrt im Vergleich zum Rest des Problems?

volkard

knivil schrieb:

Nur so als Zwischenfrage da ich auch mal ueber optimales isqrt bei Primzahlsieb nachgedacht habe: Wie gross ist der Anteil von isqrt im Vergleich zum Rest des Problems?

Allerhöchstens bei häufigen Aufrufen mit sehr kleinen Primzahlkandidaten messbar, würde ich sagen. Aber da schlägt eh if(n ist klein) return sqrt(double(n)) zu. Man braucht halt ein isqrt, was bis 2^64 genau ist.

knivil

Mein Problem ist der Ganzzahl-Newton. Wie du ihn formuliert hast, funktioniert er fuer alles bis 2^32 (grad getestet). Leider weiss ich nicht, wie ich ihn fuer 2^64 verifizieren kann. Ansonsten: Sehr nice!

volkard

knivil schrieb:

Mein Problem ist der Ganzzahl-Newton. Wie du ihn formuliert hast, funktioniert er fuer alles bis 2^32 (grad getestet). Leider weiss ich nicht, wie ich ihn fuer 2^64 verifizieren kann. Ansonsten: Sehr nice!

Eigentlich teste ich mit campers Code.
Bis 2^64 lief er bei mir durch.

camper

volkard schrieb:

Newton ganzzahlig. 860s.

UInt64 sqrt(UInt64 x){
	if(x==0) return 0;
	UInt64 sn=x>>(intlog2(x)/2);
	UInt64 so;
	do{
		so=sn;
		sn=(so+x/so)/2;
	}while(sn<so);
	if(sn*sn>x) --sn;
	return sn;
}

wie sieht es aus, wenn du

UInt64 sn=sqrt((double)x);

verwendest?

volkard

camper schrieb:

wie sieht es aus, wenn du
UInt64 sn=sqrt((double)x);
verwendest?

Viel viel viel besser.

Geht runter auf 320s oder so. Und ermöglicht die Folgenden Verbesserungen:

Nach dem guten double-Schätzwert braucht man nie mehr als eine Newton-Runde.

UInt64 sqrt(UInt64 x){//240s
	UInt64 s=std::sqrt(double(x));
	if(x>=4503599761588224ull)//2^52+2^27
		s=(s+x/s)/2;
	return s;
}

Nichtmal das! Die Wurzel passt ja noch fein in den double.

UInt64 sqrt(UInt64 x){//103s-107s
	if(x==0) return 0;
	UInt64 s=sqrt(double(x));
	if(x<4503599761588224ull)//2^52+2^27
		return s;
	if(s*s>x)
		--s;
	return s;
}

UInt64 sqrt(UInt64 x){//evtl ein halbes Sekundchen langsamer, zu ungenaue Messungen
	UInt64 s=sqrt(double(x));
	if(s*s>x)
		--s;
	return s;
}

Das werde ich produktiv einsetzen, denn wenn ich sqrt(UInt64) anwerfen, dann schleiche ich eh bei großen Zahlen rum, wo der Nachtest gebraucht wird.
Hübsch, wenn die Jagd nach Geschwindigkeit zu einfachem Code führt.

knivil

Bis 2^64 lief er bei mir durch.

Whooot? Achso, es lauft auch nur bis 2^32.

Mein Testprogramm braucht 1 Minute 29 Sekunden:

uint64_t isqrt(uint64_t x)
{
    if(x==0) return 0;
    uint64_t sn=x>>(intlog2(x)/2);
    uint64_t so;
    do{
        so=sn;
        sn=(so+x/so)/2;
    }while(sn<so);
    if(sn*sn>x) --sn;
    return sn;
}

int main()
{
    std::cerr << isqrt(1) << '\n';
    std::cerr << isqrt(2) << '\n';
    std::cerr << isqrt(3) << '\n';
    std::cerr << isqrt(4) << '\n';
    std::cerr << isqrt(5) << '\n';

    for(auto i = 6u; i < (1u<<31); ++i)
        if(isqrt(i) != (uint64_t)std::sqrt(i))
        //if((uint64_t)std::sqrt(i) != (uint64_t)std::sqrt(i))
            std::cerr << i << '\n';
}

volkard

knivil schrieb:

Bis 2^64 lief er bei mir durch.

Whooot? Achso, es lauft auch nur bis 2^32.

Mein Testprogramm braucht 1 Minute 29 Sekunden:

uint64_t isqrt(uint64_t x)
{
    if(x==0) return 0;
    uint64_t sn=x>>(intlog2(x)/2);
    uint64_t so;
    do{
        so=sn;
        sn=(so+x/so)/2;
    }while(sn<so);
    if(sn*sn>x) --sn;
    return sn;
}

int main()
{
    std::cerr << isqrt(1) << '\n';
    std::cerr << isqrt(2) << '\n';
    std::cerr << isqrt(3) << '\n';
    std::cerr << isqrt(4) << '\n';
    std::cerr << isqrt(5) << '\n';

    for(auto i = 6u; i < (1u<<31); ++i)
        if(isqrt(i) != (uint64_t)std::sqrt(i))
        //if((uint64_t)std::sqrt(i) != (uint64_t)std::sqrt(i))
            std::cerr << i << '\n';
}

campers Programm von 22 Okt 2013 17:36

knivil

Aha, ich sehe. Ich bin noch nicht ganz ueberzeugt, dass es ausreicht, nur vor und nach der Quadratzahl zu testen. Es spricht aber nix dagegen. Gut finde ich, dass isqrt auf std::sqrt reduziert wurde.

Columbo

Soll ich ihn zeigen oder willste noch ein wenig basteln?

Letzteres.

Ich weiß genau, wie deines so schnell ist: Lookup-Tables.

Gut finde ich, dass isqrt auf std::sqrt reduziert wurde.

Das finde ich hingegen ein Verbrechen. sqrt berechnet viel genauer als nötig. Wenn std::sqrt schnell ist, dann kann isqrt mit einem ähnlichen Verfahren viel schneller sein. Und das verschenkt man.

volkard

Arcoth schrieb:

Soll ich ihn zeigen oder willste noch ein wenig basteln?

Letzteres.

Ups, habs schon gezeigt.

Arcoth schrieb:

Ich weiß genau, wie deines so schnell ist: Lookup-Tables.

Jo, auch. Ein paar Milliönchen Transistoren um die ersten Bits nachzuschauen, wird man sich schon gönnen, damit ja die Framerate bei 3D-Spielen ordentlich ist und die Prozessoren von den Gamern hübsch gekauft werden.

Columbo

volkard schrieb:

Ups, habs schon gezeigt.

Wo denn das? Du hast Newton gezeigt, aber der dauert doch vier mal länger als deine Methode (800 statt 200 sek.)?

Edit: Ach, das! Na, das ist doch keine Kunst, du Betrüger!

volkard

Arcoth schrieb:

Edit: Ach, das! Na, das ist doch keine Kunst, du Betrüger!

Wärst Du drauf gekommen?

Columbo

volkard schrieb:

Arcoth schrieb:

Edit: Ach, das! Na, das ist doch keine Kunst, du Betrüger!

Wärst Du drauf gekommen?

Darauf gekommen, einfach sqrt zu nutzen, und den Fehler auszubügeln? Ich hätte das - selbst wenn es mir in den Sinn gekommen wäre - nicht umgesetzt. Ich halte std::sqrt in diesem Fall nach wie vor für unpassend, da es einfach Rechenzeit verschwendet.

Gibt es keinen wunderhübschen Algo der 100 Sekunden braucht? Und wirklich nur die gefloorte Wurzel berechnet?

knivil

Gibt es keinen wunderhübschen Algo der 100 Sekunden braucht? Und wirklich nur die gefloorte Wurzel berechnet?

Sicher, aber nicht auf deiner Wunschhardware. Aehnliche Ansaetze sind hier: http://www.azillionmonkeys.com/qed/sqroot.html