Fakultät größer 21 berechnen

[honey_bunny]

Hallo miteinander,

ich habe folgendes Problem: Mein Programm zur Berechnung der Fakultät einer Zahl liefert nur bis 21! genaue Werte, ab 22! wird an 17. Stelle gerundet.

Info: Ich verwende den Datentyp "long double".

Gibt es einen größeren Datentyp als diesen, oder eine gute Alternative zur Berechnung?

Brauche gringend Hilfe! :confused:

Viele Grüße

Claudia

[nic_power]

Hallo,

ein "long double" zu verwenden ist in diesem Fall generell eine schlechte Idee, da sich die Fakultät mit reiner Integer-Arithmetik bestimmen läßt und Fließkommaberechnungen mit Rechenungenauigkeiten behaftet sind. Allerdings wirst Du auch mit integer-Rechnungen schnell den Gültigkeitsbereich überschreiten. Lösen läßt sich das Problem, wenn Du eine "bignum"-Bibliothek verwendest:

http://www.swox.com/gmp/

Damit lassen sich Rechnungen beliebiger Genauigkeit durchführen. Falls es nur um die Berechnung geht, kannst Du auch ein Algebra-Tool wie "Mupad" einsetzen.

Nic

[kLeiner_HobBes]

Ähm .. ich würde für die Fakultät einen ganzzahligen Datentyp nehmen, z.B. long int. Der hat nämlich bis zu 64 Bit, und das reicht für ein bissel mehr als nur 22! (Was noch locker im 32-Bit-Bereich liegt).

Und aufpassen, dass du durchgehend den gleichen Datentyp verwendest!

HTH

Benjamin

[Deagle--Knight]

da int auch so schön komma zahlen darstellen kann, was bei der fakultät ja besonders wichtig ist ?

[nic_power]

Hallo,

wie so ist das bei der Fakultät wichtig? Ganz im Gegenteil. Bei der Fakultät treten keinerlei Kommastellen auf, da es sich um eine reine Integer-Berechnung handelt.

Nic

[nic_power]

Hallo,

Original geschrieben von kLeiner_HobBes

Ähm .. ich würde für die Fakultät einen ganzzahligen Datentyp nehmen, z.B. long int. Der hat nämlich bis zu 64 Bit, und das reicht für ein bissel mehr als nur 22! (Was noch locker im 32-Bit-Bereich liegt).

Ahem. Die Fakultät von 22 liegt nicht mal mehr ansatzweise im Bereich von 32 Bit, selbst 64 Bit reichen nicht aus, um die Fakultät von 22 aufzunehmen.

Nic

[FinalFantasy]

Ich denke da verwechselt jemand Multiplizieren mit Addieren. *gg*

[kLeiner_HobBes]

uuups!! Sry, hab mich doch tatsächlich vertippt :floet: :floet:

liegt im 70 bit-Bereich (69,...)

*schäm*

[Crush]

Eine Multiplikation auf x-beliebig viele Bits selber zu programmieren ist eigentlich sehr einfach. Binär ist das relativ simpel zu handhaben und Kommastellen stellen auch kein Problem dar (da werden die Bits im Fall von 0,irgendwas nur in die andere Richtung geshiftet, was etwas Mehraufwand darstellt, aber in diesem Fall ja eh nicht gebraucht wird), wenn man z.B. die ersten 64 Bit für die Kommastellen verwendet und die nächsten 256 Bits für den Ganzahlteil. Die Division ist da schon etwas komplizierter, aber wenn man weiß wie, geht das auch sehr einfach mit ein paar Zeilen Code.

[nic_power]

Hallo,

Warum sollte man das Rad zweimal erfinden, wenn vollständige Bignum-Bibliotheken zur Verfügung stehen? Aber poste doch einfach mal Deinen Code zur Multiplikation/Division beliebig langer Integer Zahlen. Das würde dem Fragesteller sicherlich weiterhelfen (sofern er auf eine zusätzliche Bibliothek verzichten möchte).

Nic

[Crush]

Wie Du willst: Hier erstmal die Multiplikation auf 64 Bit (eine Anpassung auf x-beliebig viele Bits ist ja beschrieben):

	// Sollen die Register größer sein, muß man nur die Rx-Register um entsprechend viele zusätzliche Register erweitern

	// und das Rollen und die Addition um diese ergänzen

	unsigned int RA,RB1,RB2,RC1,RC2,TMP;

	RA = 10;

	RB1=0x20000001;

	RB2= RC1= RC2= TMP = 0;

	while ( RA )

		{

		// Bei 1er Bits RB zu RC addieren

			if ( RA & 1)

			{

				TMP = RC1;

				RC1 += RB1;

				RC2 += RB2;

				// C kann leider nicht mit Übertrag addieren, deshalb die TMP-Variable (wäre nur 1 Assembler-Befehl)

				if ( RC1 < TMP && RB1 != 0)

				RC2++;

			}

			RA >>= 1;

			// solange noch Bits da sind Multiplikator * 2 (=1x nach links rollen)

			if ( RA )

			{

				// C kann leider nicht rollen (wäre auch nur 1 Assembler-Befehl)

				RB2 = RB2 << 1 | RB1 >> 31;

				RB1 <<= 1;

			}

		}

	CString Erg;

	Erg.Format("Ergebnis (hexadezimal): %p%p dez hi: %i (*4294967296) + dez lo: %i", RC2, RC1, RC2, RC1);

	MessageBox(Erg);

[/PHP]

Von der Division hab ich leider nur noch die Assembler-Version einigermaßen funktionstüchtig (ist glaub noch irgendwo der Wurm drin), aber sollte ich die wesentlich kürzere C-Version evtl. noch korrigiert bekommen, kann ichs - falls gewünscht - noch nachreichen:

[PHP]
typedef struct {unsigned long hi; unsigned long lo;} i64;
// 64-Bit Division (64-Bit-Zahl/64-Bit Divisor/64-Bit Rest...) Rest steht im Hilfsregister
i64 divident,ergebnis,divisor,hilfsregister;
divident.lo=0x23462;divident.hi=0x0;
// divisor.lo=0x80000000;divisor.hi=0;
divisor.lo=0x3;divisor.hi=0x0;
ergebnis.lo=1;ergebnis.hi=0;
hilfsregister.lo=0;hilfsregister.hi=0;
_asm
{
push eax // divisor
push ebx // divident
push ecx // hilfsregister
push edx // ergebnis

mov eax,divisor.lo
mov ebx,divisor.hi

xor ecx,ecx
bsr cx,divident.hi
jne clp1
mov ecx,divident.lo
mov divident.hi,ecx
xor cx,cx
mov ergebnis.hi,1
jmp noclp
clp1:
rol divident.lo,1
rcl divident.hi,1
rol ergebnis.lo,1
rcl ergebnis.hi,1
loop clp1
noclp:
xor ecx,ecx
bsr ecx,divident.hi
jne noclp2
clp2:
rol divident.lo,1
rcl divident.hi,1
rol ergebnis.lo,1
rcl ergebnis.hi,1
loop clp2
noclp2:



xor ecx,ecx

mov ecx,hilfsregister.lo
mov edx,hilfsregister.hi
lp:
rol divident.lo,1
rcl divident.hi,1
rcl ecx,1
rcl edx,1
cmp edx,ebx
jc h1
cmp ecx,eax
jc h1
sub edx,ebx
sbb ecx,eax
clc
h1:
cmc
rcl ergebnis.lo,1
rcl ergebnis.hi,1
jnc lp
mov hilfsregister.lo,ecx
mov hilfsregister.hi,edx
pop edx
pop ecx
pop ebx
pop eax
}
TRACE("Ergebnis: %p Rest: %p",ergebnis.lo,hilfsregister.lo);

Wer allerdings meinen Code verwendet soll mich irgendwo im Sourcecode erwähnen =8-)

[nic_power]

Hallo,

naja, wirklich hilfreich ist das nicht. Zum einen reichen für das oben angesprochene Problem 64 Bit nicht aus und zum anderen benötigt man für die Division/Multiplikation von 64 Bit Zahlen keine eigene Implementierung, da dies mit den Standarddatentypen in C/C++ wesentlich einfacher (und portabler/schneller) geht, was spricht gegen die Verwendung von "long long" bzw. "int" (je nach Architektur)? Gängige Architekturen erledigen die Multiplikation/Division mit einer Instruktion in einem Takt.

Der Schritt zu einer Arithmetik mit beliebiger Länge ist da deutlich größer, da Du um die Verwendung von Arrays nicht herum kommst.

Nic

[kLeiner_HobBes]

Ich würde auch die bignum-Bibliothek verwenden. Vor allem bist du da, etwas überspitzt ausgedrückt, wirklich nur noch auf die Größe deines Arbeitsspeichers begrenzt. Und ich denke, dass du, wenn du so etwas selber implementieren wolltest, es kaum so performant hinkriegen würdest wie GMP.

Viel Spaß beim rumprobieren

Benjamin

[Crush]

Ist das jetzt ausreichend? (96*128 Bit-Multiplikation->128 Bit Ergebnis)

Darum gehts mir ja grad: Diese Methode ist performanter als ein arraybehaftetes Bignum-Verfahren. In optimiertem Assembler läuft alles übrigens nochmal mindestens um den Faktor 10 schneller! Ich werd glaub mal Speed-Tests machen, was schneller ist, diese Methode oder GMP mit 128 Bit.

// Sollen die Register größer sein, muß man nur die Rx-Register um entsprechend viele zusätzliche Register erweitern

// und das Rollen und die Addition um diese ergänzen

	unsigned int RA1,RA2,RA3,RB1,RB2,RB3,RB4,RC1,RC2,RC3,RC4,TMP1,TMP2,TMP3;

	RA1=RA2=RA3=RB1=RB2=RB3=RB4=RC1=RC2=RC3=RC4=TMP1=TMP2=TMP3=0;

	RA1 = 0x1234;

	RB4 = 0x1234;

//	RB1=0x20000001;

	while (RA1|RA2|RA3)

		{

		// Bei 1er Bits RB zu RC addieren

			if ( RA1 & 1)

			{

				TMP1 = RC1;

				TMP2 = RC2;

				TMP3 = RC3;


				RC1 += RB1;

				// C kann leider nicht mit Übertrag addieren, deshalb die TMP-Variable (wäre nur 1 Assembler-Befehl)

				if ( RC1 < TMP1 && RB1)				RC2++;

				if ( RC2 < TMP2 && RB2)				RC3++;

				if ( RC3 < TMP3 && RB3)				RC4++;

				TMP2 = RC2;

				TMP3 = RC3;


				RC2 += RB2;

				if ( RC2 < TMP2 && RB2)				RC3++;

				if ( RC3 < TMP3 && RB3)				RC4++;

				TMP3 = RC3;


				RC3 += RB3;

				if ( RC3 < TMP3 && RB3)				RC4++;

				RC4 += RB4;

			}

			RA1 = RA1>>1 | RA2 << 31;

			RA2 = RA2>>1 | RA3 << 31;

			RA3 >>= 1;

			// solange noch Bits da sind Multiplikator * 2 (=1x nach links rollen)

			if ( RA1 | RA2 | RA3 )

			{

				// C kann leider nicht rollen (wäre auch nur 1 Assembler-Befehl)

				RB4 = RB4 << 1 | RB3 >> 31;

				RB3 = RB3 << 1 | RB2 >> 31;

				RB2 = RB2 << 1 | RB1 >> 31;

				RB1 <<= 1;

			}

		}

	CString Erg;

	Erg.Format("Ergebnis (hexadezimal): %p-%p-%p-%p dez hi3: %i (*4294967296hoch3) +dez hi2: %i \

(*4294967296hoch2) +dez hi1: %i (*4294967296) + dez lo: %i", RC4,RC3, RC2, RC1, RC4,RC3, RC2, RC1);

	MessageBox(Erg);

[/PHP]

...außerdem glaube ich, daß der Code manchen zu Lehrzwecken auch nützlich sein kann

So, jetzt ist aber Schluß für heute mit Binärarithmetik

[nic_power]

Hallo,

Original geschrieben von Crush

Ist das jetzt ausreichend? (96*128 Bit-Multiplikation->128 Bit Ergebnis)

Darum gehts mir ja grad: Diese Methode ist performanter als ein arraybehaftetes Bignum-Verfahren. In optimiertem Assembler läuft alles übrigens nochmal mindestens um den Faktor 10 schneller! Ich werd glaub mal Speed-Tests machen, was schneller ist, diese Methode oder GMP mit 128 Bit.

Wenn Du keine Arrays verwendest, wie übergibst Du die Datentypen an die Funktionen und wie werden die Ergebnisse an den Aufrufer zurückgeliefert? Irgendeinen Datentyp musst Du ja verwenden.

An den Benchmark-Ergebnissen im Vergleich zu Bignum wäre ich interessiert (insbesondere für große Zahlen >1024 Bit).

Nic

[Crush]

Wenn Du keine Arrays verwendest, wie übergibst Du die Datentypen an die Funktionen und wie werden die Ergebnisse an den Aufrufer zurückgeliefert? Irgendeinen Datentyp musst Du ja verwenden.

Wie man sieht hab ich hier alles einfach in Rx1-Rx3 oder Rx4 gesplittet, was man natürlich auch als Array laufen lassen kann. Wie gesagt ist alles in Assembler wesentlich einfacher zu realisieren und ich werd das wohl irgendwann mal auch noch rein interessehalber auf 1024 Bit ummodeln, das bringt natürlich dann nur dann benchmarkmäßig etwas, wenn man auch tatsächlich mit Zahlen multipliziert, die sich in den oberen 1000 Bit-Bereichen aufhalten, weil ja nur von rechts an bis zum letzten Bit im Divisor hochgearbeitet wird.

[nic_power]

Hallo,

Darum gehts mir ja grad: Diese Methode ist performanter als ein arraybehaftetes Bignum-Verfahren. In optimiertem Assembler läuft alles übrigens nochmal mindestens um den Faktor 10 schneller!

Wie hast Du denn die Performance Deiner Lösung gegenüber GMP gemessen, ich kann das nicht mal ansaztweise nachvollziehen? Ich habe gerade mal einen kleinen Benchmark laufen lassen. GMP ist fast 10x schneller (und baut seinen Vorsprung deutlich aus, wenn die Zahlen länger werden)! Ich habe dafür die Beispielrechnung aus Deiner Implementierung verwendet (also 0x1234*0x1234) mit 100 Mio Durchläufen.

Bench_BN , 100000000 Loops: 3530000 Ticks Time: 3530 ms (0x1234 x 0x1234)

Bench_BN , 100000000 Loops: 9160000 Ticks Time: 9160 ms (0x12341234123412341234123412341234 x 0x12341234123412341234123412341234)

Bench_BN_FAC , 100000000 Loops: 113700000 Ticks Time: 113700 ms

Bench_INT32 , 100000000 Loops: 100000 Ticks Time: 100 ms

Bench_INT , 100000000 Loops: 32150000 Ticks Time: 32150 ms (0x1234 x 0x1234)

Bench_INT , 100000000 Loops: 215060000 Ticks Time: 215060 ms (0x123412341234123412341234 x 0x1234)

(BN=GNU GMP, INT32=Native 32 Bit Multiplikation, INT=Deine Implementierung).

Da GMP (siehe unten) eine Optimierung bzgl. der benötigten Genauigkeit automatisch durchführt, habe ich sicherheitshalber auch noch einem Benchmark mit einer 128x128Bit Multiplikation durchgeführt. Ich denke die Ergebnisse sprechen für sich (insbesondere was die Skalierbarkeit betrifft wenn die Zahlen größer werden). Den Code für den Benchmark gibts hier: http://www.nleymann.de/misc/fac.c (wäre eventuell etwas viel, um ihn hier komplett zu posten). Wenn man das "#define PRINT_OUT 0" weg kommentiert sieht man übrigens sehr schön warum es so schwer ist vernünftige Benchmarks zu implementieren, wenn optimierende Compiler eingesetzt werden (aber das ist eine andere Diskussion).

Compiliert wurde der Code mit dem gcc auf einem Opteron-System unter Linux (allerdings nur 32 Bit Arithmetik): gcc -m32 -O2 -o fac fac.c -lgmp Auf einer SPARC unter Solaris sehen die Ergebnisse identisch aus.

Wie gesagt ist alles in Assembler wesentlich einfacher zu realisieren und ich werd das wohl irgendwann mal auch noch rein interessehalber auf 1024 Bit ummodeln, das bringt natürlich dann nur dann benchmarkmäßig etwas, wenn man auch tatsächlich mit Zahlen multipliziert, die sich in den oberen 1000 Bit-Bereichen aufhalten, weil ja nur von rechts an bis zum letzten Bit im Divisor hochgearbeitet wird.

Ist schon klar. GMP erledigt dies von Hause aus und arbeitet nur mit der maximal notwendigen Genauigkeit die bei Bedarf angepasst wird. Es wäre ja auch nicht sinnvoll, eine 32 Bit Addition mit 1024 Bit Genauigkeit durchzuführen. Sowas sollte man generell bei einer Implementierung mit berücksichtigen.

@honey_bunny: (um nicht völlig OT zu sein, vielleicht sollten wir ja einen extra Thread öffnen)

Wenn Du die gmp Bibliothek verwendest, kannst Du die Fakultäten mit beliebiger Genauigkeit sehr einfach ausrechnen lassen:

#include <gmp.h>


main (int argc, char **argv)

{

  mpz_t p;

  unsigned int counter;


  mpz_init_set_str(p,"1",10);

  for(counter=1; counter<=atoi(argv[1]); counter++) {

    mpz_mul_ui (p, p, counter);

  }

  gmp_printf("%Zd\n", p);

}

[/php]

Compilieren läßt sich das ganze mit "gcc -m32 -O2 -o fac fac.c -lgmp. Aufruf über "./fac 22" (zum Beispiel).

Nic

[Shadax]

Original geschrieben von kLeiner_HobBes

Ähm .. ich würde für die Fakultät einen ganzzahligen Datentyp nehmen, z.B. long int. Der hat nämlich bis zu 64 Bit [...]

Wer sagt das?

[Shadax]

Original geschrieben von nic_power

unsigned int counter;

// [...]

for(counter=1; counter<=atoi(argv[1]); counter++) {

// [ ...]

}

Du vergleichst hier einen vorzeichenbehatfteten Datentyp mit einem vorzeichenlosen... möge der Herr im Himmel gnädig sein... ;-)

[nic_power]

Hallo,

da muss der Herr im Himmel nicht gnädig sein, da die C-Erfinder dankenswerterweise an implizite Typumwandlungen gedacht haben. Es gibt in C eine klare Regelung, wie mit unterschiedlichen Integer-Datentypen bei Zuweisungen oder Vergleichen umgegangen wird.

Problematisch wird das ganze nur, wenn der Zieldatentyp bei der Konvertierung überlaufen kann (bei guten Compilern kann man entsprechende Warnings erzeugen lassen: -Wsign-compare beim gcc beispielsweise).

Nic

[Shadax]

Original geschrieben von nic_power

Problematisch wird das ganze nur, wenn der Zieldatentyp bei der Konvertierung überlaufen kann .

Genau -deswegen- sollte ja der Herr im Himmel ja gnädig sein (damit der User keine Werte < 0 übergibt) ;-)

[nic_power]

Hallo,

Original geschrieben von Shadax

Genau -deswegen- sollte ja der Herr im Himmel ja gnädig sein (damit der User keine Werte < 0 übergibt) ;-)

Das hat aber nichts mit dem Vergleich zu tun sondern mit der generellen Behandlung von Fehleingaben durch den Nutzer. Die Fakultät ist mathematisch nur für natürliche Zahlen definiert (>=0, mit 0!=1).

Nic

[honey_bunny]

Hallo,

ich bin glaub ich einfach zu blöd: Ich hab mir dieses gmp-Dingsda runtergeladen, aber das ist ein riesiger Ordner und ich weiß nicht was ich damit machen soll...

Kann mir jemand beschreiben wie ich als Dummy da vorgehen muss? (Benutze Visual Studio 6.0 pro)

Vielen Dank für die super Antworten

[nic_power]

Hallo,

prinzipiell muss Du gmp selbst auf der Ziel Plattform übersetzten. Alternativ kannst Du Dir auch mal die folgende Seite ansehen, dort gibt es vorkompilierte Binaries für verschiedene Windows-Compiler:

http://www.cs.nyu.edu/exact/core/gmp/

Nic

[honey_bunny]

Also, irgendwie klappt das nicht...

Ich hab mir jetzt die gmp.h in das Visual Studio reingepackt, aber welchen Datentyp muss ich denn da nehmen?

mpz_t kann es schonmal nicht sein, weil das Teil inkompatibel mit long double ist. Ich muss mit den Werten tierisch viel rumrechnen (Hypergeometrische Verteilung)!

Würde ja gern den Quelltext posten, aber der ist mächtig lang ;-)

Kann mir nochmal jemand helfen, ich krieg hier noch die totale Krise

Danke schön