quotient aus fakultäten (im datentypbereich belassen)



  • hallöchen,
    habe ein programm geschrieben was mir die mikrokanonischen zustandssummen für alle möglichen energien bei vorgegebener teilchen zahl geschrieben (statistische physik)
    doch das ist nicht das problem 🙂
    das programm funktioniert hervorragend, allerdings nur bis zu N=170, ab 171 gibt er statt der gewünschten zahlen nur "1.#INF" aus, da bei 170 die größte potenz 79 ist und die möglichkeiten von double meines wissens nach weit darüber liegen, kanns doch eigentlich nicht direkt an der größe von double liegen oder ??
    leider fehlt mir die erfahrung und so musste ich nach einiger suche doch die segel streichen und wende mich an euch
    ich arbeite mit ms vcpp
    hier der code

    #include <iostream.h>
    
    int betrag(int m)
    {
    	if(m<0)
    		return -m;
    	else
    		return m;
    }
    double fak(double c)
    {
    	double a=c;
    	if(a<2)
    		return 1;
    	for(double i=a-1;i>1;i--)
    		a=a*i;
    	return a;
    }
    
    double zustaende(int a,int b)
    {
    	double m = a;
    	int N = b;
    	double summe=0;
    	double x=0,y=0,z=0;				//x=Anzahl -1, y=Anzahl 0, z=Anzahl 1
    
    	if(m==0)y=N;
    
    	if(m<0)
    	{
    		x=betrag(m);
    		y=N+m;
    	}
    	if(m>0)
    	{
    		z=m;
    		y=N-m;
    	}
    
    	if(m>=-N&&m<=N)
    	{
    		for(double i=0;i<((N-betrag(m))/2)+1;i++)
    		{
    			summe+=(fak(N)/(fak(x)*fak(y)*fak(z)));
    			x+=1;y-=2;z+=1;
    		}
    	}
    	else return 0;
    	return summe;
    }
    
    void main()
    {
    	int N=1;
    	while(N!=0)
    	{			//Teilchenzahl
    	cout<<"\n"<<"Teilchenzahl eingeben:";
    	cin>>N;
    
    	for(int i=-N;i<=N;i++)
    	{
    		cout<<"m = "<<i<<":		"<<zustaende(i,N)<<"\n";
    	}
    	}
    
    }
    

    viele dank für eure hilfe
    grüße
    dr_ceran



  • Ne direkte Antwort hab ich darauf auch nicht, aber schau dir mal Folgendes an:

    http://www.datasource.de/programmierung/tab33_cppdatentypen.html

    Kannst ja mal den 80-bit-Datentyp bei 32-bit-Systemen versuchen. Hab aber keine Ahnung ob dass klappt.



  • sorry, aber des funzt leider auch net, was aber nicht wirklich verwunderlich ist, da ja bei N=170 die höchste potenz 79 in den berechnungen ist, bei 171 müsste es dann bei höchstens 80 in der potenz liegen...
    und 80 ist ja für double noch kein problem, denn dieser datentyp geht ja eigentlich bis zur potenz 308 ....
    in einem testprogramm in dem ich mir einige double zahlen mal habe ausgeben lassen, kam ich auch bis zu potenz 308....
    wenn jemand eine idee hat wo der fehler im code zu suchen ist ??
    denn es scheint mittlerweile ganz eindeutig am code und nich am double zu liegen



  • Warum verwendest du überhaupt double? soweit ich das sehe rechnest du doch mit Natürlichen Zahlen.

    versuch doch mal

    long int
    

    für Zahlen von (-2 147 483 648 .. 2 147 483 647)



  • dr_ceran schrieb:

    int N=10.2;
    	while(N!=0)
    

    Also für einen int-Wert eine 10,2 einzugeben finde ich schonmal ein bisschen seltsam.


  • Mod

    das erwartete ergebnis mag vielleicht im darstellbaren bereich von double liegen, aber das gilt ja nicht unbedingt für alle zwischenergebnisse. was wäre denn fak(171) ? ... sobald du eine unendlichkeit an irgendeiner stelle hast, bekommst du sie nicht mehr durch einfache rechenoperationen weg (das ist ja auch der sinn dabei). wenn du also - ohne auf grössere datentypen auszuweichen - mit gtrösseren N rechnen willst, müsstest du darüber nachdenken, evtl. statt jede fakultät einzeln zu berechnen, diese bereits mit der folgenden division zu kombinieren. eine andere möglichkeit wäre, für fak(N) nur für die zwischenrechnung einen vorfaktor von 2^-1023 zu wählen, und ihn zum schluss wieder eliminieren. dann erfolgt der überlauf erst etwas später.



  • Was erhoffst du dir davon:

    return 000
    

    Die Zahl 0 oder den String "000"?


  • Mod

    eine weitere möglichkeit wäre, alle zahlen auf einen bereich von 1bis 2 zu normieren und die normierungsfaktoren (2er potenzen) separat zu verechnen. etwa so:

    double fak(int c, int& norm)
    {
        double f = 1;
        norm = 0;
        for( int i = 2; i <= c; ++i )
        {
            f *= i;
            while( f >= 2 )
            {
                f /= 2;
                ++norm;
            }
        }
        return f;
    }
    

    die gesuchte zahl wäre hier also f*2^norm. mit norm kann man dann ganz normal (entsprechen potenzgesetzen) weiter rechnen; erst ganz zum schluss ist es sinnvoll, es mit der mantisse f zu kombinieren. andere faktoren als 2 sind denkbar um den rechenaufwand zu reduzieren. hier nur zur veranschaulichung der idee.



  • @camper:
    genau das ist es fak(171)= infinity...
    dein lösungsansatz sieht vielversprechend aus...
    nochmals vielen dank für die hilfe
    sorry wegen der seltsamen fehler im code, aber es war spät und ich hab zich was probiert und nich alles wieder verändert 😞



  • Gast221212 schrieb:

    Was erhoffst du dir davon:

    return 000
    

    Wenigstens ein oktales Literal im Code? 😃



  • @camper:
    hab jetzt mit deiner idee bzgl. der normierung ein wenig rumprobiert, und wenns nur um eine fakultät geht is das natürlich gold wert auch wenns um einen quotienten zweier fakultäten geht...
    das problem is aber, dass ja in meiner fkt folgender term auftaucht...

    (fak(N))/(fak(x)*fak(y)*fak(z))
    

    d.h. hier wird die normierung wieder zum fallstrick...
    die idee ist daher denk ich vielleicht im "dynamischen" berechnen dieses quotienten zu finden, also quasi die einzelnen fakultäten wachsen lassen und nach jedem wachstumsschritt das produkt speichern, danach nur noch die hinzukommenden faktoren (die also bei jeder einzelnen fakultät dazukämen) wieder als quotient errechnen und dem alten produkt auffaktorisieren....
    dann is man denk ich ständig im "grünen" bereich...
    falls jemand zu dieser idee (da die ja vielleicht nich neu ist) den code schon hat dann bitte hier posten...
    ansonsten werd ich mich jetzt mal selbst versuchen


  • Mod

    ist doch simpel (Potenzgesetze): multiplikation wird zur addition der exponenten, division zu subtraktion etc.; also
    N!x!y!z!=N!normiert2norm_Nx!_normiert2norm_xy!_normiert2norm_yz!_normiert2normz\frac{N!}{x!y!z!}=\frac{N!_{normiert}2^{norm\_N}}{x!\_{normiert}2^{norm\_x}y!\_{normiert}2^{norm\_y}z!\_{normiert}2^{norm_z}}

    =N!normiertx!normierty!normiertz!normiert2norm_Nnorm_xnorm_ynorm_z=\frac{N!_{normiert}}{x!_{normiert}y!_{normiert}z!_{normiert}}2^{norm\_N-norm\_x-norm\_y-norm\_z}



  • @camper:
    genau so geht`s, so leid es mir tut, eben nich...=)
    wenn du das mal probierst so kommt bei dem quotient von zwei fakultäten (mit ausreichend hohen N) noch was ordentliches raus...logisch, da sich dann die normierungen gerade rauskürzen...
    aber im nenner hast du nochmal die normierung^2 und das überschreitet den double bereich gerade wieder....
    aber meine idee sieht momentan ganz vielversprechend aus, wenns fertig is werd ich es posten...
    solltest du lust haben, kannst du ja deinen ansatz mal ausbauen, vielleicht hab ich nur was nich verstanden... wenn deine lösung bei N=300 immer noch nich infinity ausgibt, dann hast dus und ich verneige mich in erfurcht


  • Mod

    #include <iostream>
    
    typedef double real;
    
    real fak(int c, int& norm)
    {
    	real f = 1.0;
    	norm = 0;
    	for( int i = 2; i <= c; ++i )
    	{
    		f *= i;
    		while( f >= 4294967296.0 )
    		{
    			++norm;
    			f /= 4294967296.0;
    		}
    	}
        return f;
    }
    
    real combine(real f, int norm)
    {
    	if( norm < 0 )
    	{
    		while( norm++ ) f /= 4294967296.0;
    	}
    	else
    	{
    		while( norm-- ) f *= 4294967296.0;
    	}
    	return f;
    }
    
    real zustaende( int m, int N )
    {
        int x = 0, y = 0, z = 0;
        if( m == 0 ) y = N;
        if( m < 0 )
        {
    		x = -m;
    	    y = N + m;
        }
        if( m > 0 )
        {
            z = m;
    	    y = N - m;
        }
    	real summe = 0.0;
        if( m >= -N && m <= N )
        {
            for( int i = 0; i < ( ( N - abs( m ) ) / 2 ) + 1; ++i )
            {
    			int normN, normx, normy, normz;
    			real fakN = fak( N, normN );
    			real fakX = fak( x, normx );
    			real fakY = fak( y, normy );
    			real fakZ = fak( z, normz );
    			summe += combine( fakN / ( fakX * fakY * fakZ ), normN - ( normx + normy + normz ) );
                ++x;
    			y -= 2;
    			++z;
            }
        }
        return summe;
    }
    
    int main()
    {
    	int N = 1;
    	while( N != 0 )
        {
    		std::cout<< "Teilchenzahl eingeben:";
        	std::cin >> N;
    	    for( int i = -N; i <= N; ++i )
        	{
            	std::cout << "m = " << i << "  :\t" << zustaende( i, N ) << std::endl;
        	}
        }
        return EXIT_SUCCESS;
    }
    

    funktioniert bis einschliesslich N=649. wenn du g++'s long double benutzt, geht es noch ein ganzes stück weiter - allerdings ist die normierungsfunktion recht ineffizient und das macht sich dann bemerkbar 😉 habe hier bereits eine grössere potenz genommen, um es nicht bei kleinen zahlen schon zum problem werden zu lassen. es sollten aber unbedingt 2er potenzen sein, alles andere würde zusätzliche rundungsfehler produzieren.



  • @camper:
    kopfschüttelnd betrachtete ich deinen code, und fragte mich, wie ich es mir erlauben konnte an deiner idee zu zweifeln 🙂
    in meiner beschränktheit übersah ich doch glatt die möglichkeit, die normierungen klein zu wählen, dafür häufiger auszuführen und zu zählen, um sie dann ganz ähnlich immer im datentypbereich bleibend wieder zu entfernen....
    wunderschön kann ich da nur sagen...
    der code gefällt mir sehr gut...
    hab wieder ne menge gelernt...
    und darum nochmals danke für deine rege anteilnahme an meinem problem (und vor allem dessen lösung 😃 )



  • Was willst Du eigentlich genau ausrechnen? Wenn Dein Endergebnis platz in nem unsigned long long hat, dann kannst Du auch mit denen Rechnen, indem Du modulo einer ausreichend großen Primzahl rechnest. Das dürfte um einiges Effizienter sein als mit den doubles.

    MfG Jester


  • Mod

    @Jester: die zahlen werden recht gross für grosse N, da reicht auch long long recht schnell nicht mehr aus. natürlich muss man sich überlegen, wieweit man das ausbauen will, bis man erkennt, das man mit einer echten bignum-bibliothek besser bedient ist.
    P.S. wieso muss es eigentlich eine primzahl sein?

    ich tippe mal das es um atomphysik geht. aber das ist geraten.



  • Ich weiß, daß die Fakultät wächst, aber das ist nicht so wichtig. Mir geht es um das Endergebnis. 100 über 5 ist zum Beispiel im Vergleich zu der riesigen 100! garnicht so groß. Wie gesagt es kommt nur auf die Größe des Endergebnisses an, nicht auf die der Zwischenschritte. Allein davon hängt es ab was nötig ist.

    MfG Jester


  • Mod

    ich bezog mich auch auf das endergebnis. bereits für N=43 passt es nicht mehr für alle m in den unsigned long long bereich. bei N=650 steigt er wie gesgt bei double aus. bei long double treten ab N=10343 überläufe auf.



  • @camper:
    es handelt sich um theoretische physik (genauer statistik/thermodynamik)
    ich deutete es schon im ersten post an...
    man hat eine bestimmte anzahl an teilchen im magnetfeld, welche nicht wechselwirken... diese können einen spin von -1,0,1 aufweisen (genauer gesagt die projektion eines spins auf die z-achse, wie es in der th.physik üblich ist)
    nun wird eben mit dem programm berechnet wie groß die mikrokanonische(kein energie-und teilchenaustausch) zustandssumme zu einem bestimmten energiewert wird... das m was im code vorkommt ist hier der entscheidende teil der energie, der dann noch mit zwei festen größen multipliziert wird um zu energie zu gelangen
    dieses m ist für N teilchen die summe aus deren spinprojektionen, kann also bei N=1 auch nur 1,0,-1 sein....
    d.h. es wird also die zahl der möglichkeiten gesucht um mit einer festen anzahl an teilchen auf den jeweiligen energiewert bzw m zu kommen...
    viele grüße


Anmelden zum Antworten