n über k - berechnen!?

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