n über k - berechnen!?

[Plague]

Hallo wie kann ich das mathematische "n über k" berechnen?
Kennt da jemand den Quelltext?

Gruß
Thomas

[bttb930]

n über k = n! / (k! (n-k)!)

das ist leicht zu programmieren, oder?

[Kryoko]

Hi !

Wenn ich mich nicht irre meinst du mit "n über k" die Binomialkoeffizienten (irre ich mich ??)
Berechnet werden die mit
n! / ((n - k)! * k!)

Delphi-Quellcode:

function fakultaet(N: Int64): Int64;
var 
  i: Integer;
begin
  Result := 1;
  for i := 1 to N do
    Result := Result * i
end;

function nueberk(n, k: Integer): Int64;
begin
  Result := fakultaet(n) / (fakultaet(n - k) * fakultaet(k))
end;

[xineohp]

jungs, denkt doch mal ein bisschen nach!

schon 21! sprengt mit 5,09*10^19 die int64-Grenzen (2^63 = 9,22*10^18) !!

[Kryoko]

ist mir auch klar...
der Quellcode sollte niemals 100% funktionstüchtig sein...sonst hät' ich da noch was mit kürzen etc. geschrieben.
er wollte nur Quellcode. :evil:

Oder alles mit Extended FLoats geschrieben...

[xineohp]

hm, da hast du wohl Recht ... aber gerade das ist doch das Problem an der Sache, so ist die Lösung doch witzlos, da sie in 50% der Fälle nicht geeignt ist :zwinker:

[Dani]

Delphi-Quellcode:

function Binomial(n, k: Integer): Extended; 
 begin
  If (k < 0) or (n < 0) then raise {Fantasie-Excpetion};
  If (k = 0) or (k=n) then Result := 1
   else Result := n/k * Binomial(n-1, k-1);
 end;

bzw.

Delphi-Quellcode:

function Binomial(n, k: Cardinal): Extended; 
 begin
  If (k = 0) or (k=n) then Result := 1
   else Result := n/k * Binomial(n-1, k-1);
 end;

[xineohp]

womit aber das Problem mit den zu großen Zahlen bestehen bleibt.

[Dani]

Die Grenze des Berechenbaren liegt aber viel höher als bei Int46. Für die Praxis reicht es bestimmt.

[xineohp]

eben nicht!

die Standardaufgaben zu Bernoulliketten arbeiten mit n=10 bis n=100. Und ab 21 ist int64 zu klein!

[Dani]

Zitat:

die Standardaufgaben zu Bernoulliketten arbeiten mit n=10 bis n=100. Und ab 21 ist int64 zu klein!

Der Ansatz, den ich gepostet habe kommt auch mit n = 5000 klar (wird aber irgendwann ungenau). Die Fakultät wird wegen der Rekursion nie als Ganzes berechnet.

[JasonDX]

probier das mal:

Delphi-Quellcode:

function nuberk(n, k: integer): extended;
var
  i: integer;
begin
  result := 1;
  for i := n downto 2 do
    if (i in [k + 1..n]) then
      result := result * i
    else if(i in [2..n-k]) then
      result := result / i;
end;

ich habs nicht probiert und kann auch gar nicht garantieren, dass der code funktioniert, es is eine (un)überlegte "optimierung" eines montag morgens ;)

[dizzy]

Das mit dem Mengenoperator wird schwierig sobald die Intervallgrenzen ausserhalb von 0..255 liegen. Da solltest du besser mit 2 Vergleichen dran gehen.

Gruss,
Fabian
(der auch nicht die Richtigkeit überprüft hat ^^)

[xineohp]

@Dani: Da hast du Recht ... ich hatte nicht daran gedacht, dass die Gleitkommatypen soviel größer sind, sorry.

[ibp]

desweiteren werfe ich mal das kürzen des quotienten in die runde!

Code:

/5\    1*2*3*4*5      4*5
| | = ------------ = ------
\3/    1*2*3*(1*2)    (1*2)

[Wolfgang Mix]

Ich möchte jetzt 'nen alten Thread wiederbeleben, der für mich unbefriedigend
endete. Die Wahrscheinlichkeit, einen 6er im Lotto zu landen, ist
6 über 49, also ca. 1/14 Mio. Mit meinem Aldi-Taschenrechner für 3,99
löse ich das problemlos. Da bei Int64-Werten bei 21! Schluß ist,
habe ich alle Variablen mit dem Typ Extended rechnen lassen.


Laut Embarcadero ist ...

Zitat:

The range for the Extended type
is between 3.6 x 10^-4951 and 1.1 x 10^4932. The size in bytes of an Extended value is 10.

Mit diesem Code prüfe ich das:

Delphi-Quellcode:

function fakultaet(N: Extended): Extended;
var
  i: Integer;
begin
  Result := 1;
  for i := 1 to trunc(N) do
    Result := Result * i
end;

function nueberk(n, k: Extended): Extended;
begin
  Result := fakultaet(n) / (fakultaet(n - k) * fakultaet(k))
end;

procedure TForm1.Button1Click(Sender: TObject);
begin
  Edit1.Text:=FloatToStr(nueberk(49,6));
end;

Ergebnis:

6 über 49 ---> 1,18366178998794E-60 //Falsch

49 über 6 ---> 13983816 //Richtig

Da die eine Lösung lediglich der Kehrwert der anderen Lösung ist, wundert mich das.

Frage: Was läuft hier falsch?

Grüß

Wolfgang

[jfheins]

Die Formel gilt nur für n >= k

Die Wahrscheinlichkeit für 6 Richtige (ohne Superzahl) ist also 49 über 6 = 13983816

6 über 49 gibt es nicht / ist nicht definiert oder ist gleich 0. (Je nach Geschmack/Definition)

Btw.: Was für ein Kehrwert? 14 Millionen ist nicht der Kehrwert von 1E-60 - und n mit k zu vertauschen ist kein Kehrwert sondern einfach nicht sinnvoll ;)

[Wolfgang Mix]

@jfheins

Zitat:

Die Wahrscheinlichkeit für 6 Richtige (ohne Superzahl) ist also 49 über 6 = 13983816

Das ist leider nicht richtig!

Die Wahrscheinlichkeit für 6 Richtige (ohne Superzahl) ist 6 über 49 = 1/13983816,

sonst hättest Du ja 13983816 Gewinne pro Spiel, und 6 über 49 gibt es und das ist der
Kehrwert von 49 über 6!

[jfheins]

Okay, die Wahrscheinlichkeit für 6 Richtige ist

1/(49 über 6) = 1/13983816 = 7,151124 * 10^-8

Aber 1/(49 über 6) ist nicht gleich (6 über 49) :!:

[patti]

Zitat:

Zitat von Wolfgang Mix

und 6 über 49 gibt es und das ist der
Kehrwert von 49 über 6!

Nein! 6 über 49 gibt es einfach nicht. Die Wahrscheinlichkeit eines Lotogewinns ist (49 über 6)^(-1). Siehe dazu auch Wikipedia:

Zitat:

Ist dabei k > n, so ist (n über k) = 0

(Wikipedia-Artikel "Binomialkoeffizient")

mfg

[Wolfgang Mix]

@jfheins Danke, da war mein Knick in der Denke,
war schon zu lange her, daß ich das 'mal brauchte.
Der Code ist also okay :)

@patti: hat sich geklärt, danke

Gruß

Wolfgang

[Uwe Raabe]

Zitat:

Zitat von Wolfgang Mix

Die Wahrscheinlichkeit für 6 Richtige (ohne Superzahl) ist 6 über 49 = 1/13983816

jfheins hat Recht: 6 über 49 ist nicht definiert, da bei (n über k) k <= n sein muss. Das ergibt sich allein dadurch, daß n! für negative Zahlen nicht definiert ist.

Übrigens gibt (49 über 6) nicht die Wahrscheinlichkeit eines Gewinns mit 6 Richtigen an, sondern die Anzahl der Kombinationen die man in 6 aus 49 bilden kann. Eine dieser Kombinationen ist dann der Gewinn. Deshalb ist die Wahrscheinlichkeit 1/(49 über 6). Wenn du also (49 über 6) unterschiedliche Tips abgibst, ist sicher ein 6er dabei.

[patti]

Zitat:

Zitat von Uwe Raabe

Übrigens gibt (49 über 6) nicht die Wahrscheinlichkeit eines Gewinns mit 6 Richtigen an, sondern die Anzahl der Kombinationen die man in 6 aus 49 bilden kann. Eine dieser Kombinationen ist dann der Gewinn. Deshalb ist die Wahrscheinlichkeit 1/(49 über 6)

Dies wiederum folgt direkt aus der Formel zur Berechnung einer Laplace-Wahrscheinlichkeit:

Code:

Wahrscheinlichkeit = (Anzahl der günstigen Fälle) / (Anzahl der möglichen Fälle)

da alle Kombinationen gleich wahrscheinlich sind. Weil es nur einen günstigen Fall gibt, ist die Wahrscheinlichkeit also (49 über 6)^(-1).

mfg

[Uwe Raabe]

Zitat:

Zitat von jfheins

Aber 1/(49 über 6) ist nicht gleich (6 über 49) :!:

Vollkommen richtig!

[Wolfgang Mix]

Habe den Code noch leicht abgeändert.
Die Funktion fakultät bzw nueberk funktioniert bis
n(max) = 1754! = 1,97926189010501E4930.
Das sollte für Experimente reichen.

Delphi-Quellcode:

function fakultaet(N: integer): Extended;
var  i: Integer;
begin 
  Result := 1; 
  for i := 1 to N do
    Result := Result * i
end; 

function nueberk(n, k: integer): Extended;
begin 
  Result := fakultaet(n) / (fakultaet(n - k) * fakultaet(k)) 
end;

procedure TForm4.Button1Click(Sender: TObject);
begin
  Edit1.Text:=FloatToStr(nueberk(1754,600));
  Edit2.Text:=FloatToStr(fakultaet(1754));
  //n!(max) = 1754!=1,97926189010501E4930
end;

Gruß

Wolfgang

[stoxx]

http://www.matheplanet.com/default3....ne%26spell%3D1


oder hier (auf die harte Tour)

http://www.delphi-library.de/topic_B...n_35050,0.html

[Wolfgang Mix]

Danke für die Links,

werde ich mir 'reinziehen.

Soll ich noch etwas ändern oder ist etwas noch faul?

Grüß

Wolfgang

[ibp]

oben habe ich die Möglichkeit zum Kürzen des Bruches beschrieben..

ansonsten wäre das noch möglich:

http://upload.wikimedia.org/math/5/1...5b1b818978.png

siehe hier: Wikipedia Binomialkoeffizient

[Corpsman]

Zu dem thema kann ich auch noch was beisteuern :


Binomial Source

[Wolfgang Mix]

Danke für die Beiträge.
Eine weitere Optimierung bringt ja eigentlich nur etwas,
wenn k in die Größenordnung von n kommt.
Dann dürte man allerdings nicht n! zuerst ausrechnen lassen,
und der Code wird um einiges länger.
Werde mich mit dem Thema weiter befassen.

Gruß

Wolfgang

[Corpsman]

Schau dir meine Samples von #29 an, da sind die Varianten drin, die ohne größere Zwischenergebnisse auskommen.

[Wolfgang Mix]

Danke, mache ich.

Vielleicht habt Ihr ja auch Lust, an der Optimierung mitzuarbeiten ? :)

[Amateurprofi]

Ich hab das vor langer Zeit mal so gelöst, würde das heute aber etwas anders machen.

Delphi-Quellcode:

FUNCTION Combinations(const n,k:extended):extended;
const maxk=32767;
var r1,r2:extended; i,kk,nn,len:integer; ok:boolean;
    list:array of integer;
//------------------------------------------------------------------------------
FUNCTION CheckDivisible(var n:integer; k:integer):boolean; register;
asm
      push  ebx
      mov   ebx,eax
      mov   eax,[ebx]
      mov   ecx,edx
      xor   edx,edx
      div   ecx
      or    edx,edx
      jne   @1
      mov   [ebx],eax
@1:   sete  al
      pop   ebx
end;
//------------------------------------------------------------------------------
begin
   if (k>n) then begin
      result:=0;
   end else if (k=n) or (k=0) then begin
      result:=1;
   end else if (k>maxk) then begin
      result:=NaN;
   end else begin
      nn:=Trunc(n);
      len:=Trunc(k);
      SetLength(list,len);
      r1:=1;
      r2:=1;
      // alle n in Liste
      for i:=0 to len-1 do begin
         list[i]:=nn;
         Dec(nn);
      end;
      ok:=true; // damit der Compiler nicht meckert
      // alle k gegen n aus liste kürzen, ggfs. in r1 aufmultiplizieren
      for kk:=len downto 2 do begin
         for i:=0 to len-1 do begin
            ok:=CheckDivisible(list[i],kk);
            if ok then break;
         end;
         if not ok then r1:=r1*kk;
      end;
      // alle n aus Liste in r2 aufmultiplizieren
      for i:=0 to len-1 do r2:=r2*list[i];
      // resultat zurückgeben
      result:=int(r2/r1);
   end;
end;

[Wolfgang Mix]

Sieht gut aus!
Und wie würdest Du das heute machen?

[Amateurprofi]

Vor allem schneller!
Ich werde mir das im Laufe des Wochenendes mal ansehen und (so hoffe ich) eine schnellere Version posten.
Vor allem die Prüfung auf Teilbarkeit läßt sich sicherlich besser gestalten.

[rollstuhlfahrer]

Ich hab in den einfachen Code das mit dem Kürzen mal integriert. Die Zwischenergebnisse sollten hier wesentlich kleiner sein.

Delphi-Quellcode:

function fakultaet2(start, ende: integer): extended;
var
  i: integer;
begin
  if ende > start then
    raise EMathError.Create('Start kleiner als Ende');
  Result := 1;
  for i := start to ende do
  begin
    Result := Result * i;
  end;
end;

function fakultaet(N: integer): Extended;
var  i: Integer;
begin
  Result := 1;
  for i := 1 to trunc(N) do
    Result := Result * i
end;

function nueberk(n, k: integer): Extended;
begin
  Result := fakultaet2(k, n) / (fakultaet(n - k))
end;

Bernhard


Der Code ist fehlerbehaftet und errechnet falsche Ergebnisse. Eine korrigierte Version befindet sich in Beitrag #56.

[edit=Admin]Anmerkung des Autors eingefügt. Mfg, Daniel[/edit]

[Amateurprofi]

Zitat:

Zitat von Wolfgang Mix

Und wie würdest Du das heute machen?

Ich hab das mal komplett in ASM hingestolpert.
Wesentliche Änderungen:
1) wenn N-K < K ist wird NK(N,K) als NK(N,N-K) ausgeführt.
2) die Routine die versucht Faktoren durch Kürzen zu eleminieren, arbeitet vollständig in der FPU

Besonders bei großen N und großen K ergeben sich durch (1) erhebliche Verbesserungen gegenüber der ursprünglichen Version.

Code:

NK(50,10)       1000000 x ausgeführt = 390 ms (ursprünglich 1200 ms)
NK(25000,24999)       1 x ausgeführt =   0 ms (ursprünglich 1800 ms)

Bei Fehlern wird als Ergebnis NaN zurückgegeben.
Andere Fehler wie Überschreitung des Rechenbereiches oder Out of Memory werden nicht abgefangen.

Delphi-Quellcode:

FUNCTION NK(n,k:integer):extended;
asm
               test     eax,eax
               js       @ReturnNan           // n<0, Fehler
               test     edx,edx
               js       @ReturnNan           // k<0, Fehler
               je       @Return1             // k=0, Result=1
               mov      ecx,eax
               sub      ecx,edx              // d:=n-k
               js       @Return0             // k>n, Result=0
               je       @Return1             // k=n, Result=1
               cmp      ecx,edx              // if d<k then k:=d
               jae      @1
               mov      edx,ecx
@1:            cmp      edx,1
               jne      @2                   // k<>1
@ReturnN:      push     eax
               fild     dword [esp]
               pop      eax
               jmp      @End
@ReturnNan:    lea      eax,@Nan.Pointer[eax]
               fld      tbyte [eax]
               jmp      @End
@NaN:          dw       0,0,0,$C000,$FFFF
@Return0:      fldz
               jmp      @End
@Return1:      fld1
               jmp      @End
               //---------------------------------------------------------------
@2:            push     ebx                  // EBX retten
               // Platz für k Integers auf Stack reservieren
               lea      ecx,[edx*4]
               sub      esp,ecx
               // n, n-1, n-2, ... n-k+1 auf Stack legen
               sub      ecx,4
@Loop1:        mov      [esp+ecx],eax
               sub      eax,1
               sub      ecx,4
               jnc      @Loop1
               //---------------------------------------------------------------
               // alle k gegen n aus liste kürzen, ggfs. in st aufmultiplizieren
               //---------------------------------------------------------------
               push     edx                  // k
               fld1
               fld1
               fild     dword [esp]          // st0=k, st1=Const 1, st2=Produkt(k)
               lea      eax,[edx-1]
@Loop2:        mov      ecx,edx              // Index Listenende
@Loop3:        // prüfen, ob ein Wert aus der Liste durch k teilbar ist
               fild     dword [esp+ecx*4]    // st0=n
               fdiv     st,st(1)             // st0=n/k
               fld      st                   // dto
               frndint                       // st0=Int(n/k)
               fcomip   st,st(1)
               je       @Divisible           // n durch k teilbar
               ffree    st(0)
               fincstp
               sub      ecx,1
               jne      @Loop3
               fmul     st(2),st             // k zu st(1) multiplizieren
@NextK:        fsub     st,st(1)             // k:=k-1
               sub      eax,1
               jne      @Loop2               // weiter bis k=1
               jmp      @Finish
@Divisible:    // Ist teilbar. wenn quotient=1 dann n aus Liste entfernen,
               // sonst Quotient in Liste
               fcomi    st,st(2)
               je       @Quotient1           // Quotient ist 1
               fistp    dword [esp+ecx*4]    // Quotient in Liste
               jmp      @NextK
@Quotient1:    // Qoutient=1, n aus Liste entfernen
               ffree    st(0)
               fincstp
               mov      ebx,[esp+edx*4]
               mov      [esp+ecx*4],ebx
               sub      edx,1
               jmp      @NextK
@Finish:       // Alle n aus Liste aufmultiplizieren
               fcomp                        // k vom FPU-Stack nehmen
@Loop4:        fimul    dword [esp+edx*4]
               sub      edx,1
               jne      @Loop4
               // Und durch Produkt(k) dividieren
               fdivrp
               // Liste vom Stack nehmen
               pop      edx
               lea      esp,[esp+edx*4]
               pop      ebx
@End:
end;

[Corpsman]

@Amateurprofi
Warum liefert deine funktion einen Extended ?

n über k ist immer ein Integer, oder Int64 = Ganze Zahl ?

[Wolfgang Mix]

@Corpsman:

Int64 kann nur 21! Mit Extended schaffe ich immerhin 1754!

[Amateurprofi]

Zitat:

Zitat von Corpsman

@Amateurprofi
Warum liefert deine funktion einen Extended ?

n über k ist immer ein Integer, oder Int64 = Ganze Zahl ?

Weil man vielleicht auch mal größere Werte berechnen möchte.