Silnia - szybkie DIV i MOD

---

Czas wykonania instrukcji arytmetycznych (w mikro sekundach).
Jak widać zawsze całkowitoliczbowe dzielenie dział najsłabiej.

Ma może ktoś dane o bardziej współczesnych procesorach?

Przy okazji obliczeń wielokrotnej precyzji, często pojawia się fragment

      // Metoda 1
      c:=  c + t[i]*X;
      t[i]:=c mod 10000;
      c:=c div 10000;

Łatwo można pozbyć się jednej operacji dzielenia

      // Metoda 2
      c:=  c + t[i]*X;
      r:=c div 10000;
      t[i]:=c-r*10000;
      c:=r;

Mało tego, dzielenia można pozbyć się w ogóle

      // Metoda 3
      c:=  c + t[i]*X;
      asm
        movl  c,%eax
//        incl  %eax            ;konieczne gdy frac(T)<0.5
        movl  $-776530087,%edx
        mull  %edx
        shrl  $13, %edx
        movl  %edx,r
      end ['EAX','EDX'];
      t[i]:=c-r*10000;
      c:=r; 

Są kompilatory, które potrafią zastąpić dzielenie przez stałą (nie koniecznie każdą) mnożeniem i przesunięciem. Warto sprawdzić swój kompilator.
Dzielenie C przez 10'000, wykonane jest następująco: (C * $D1B71759) shr (13+32).
Napisałem 13+32, tak naprawdę to wykonywane jest przesunięcie w prawo o 13 bitów starszego słowa 64 bitowego wyniku mnożenia, i właśnie to jest zwracaną wartością.
Jak obliczyć magiczną stałą $D1B71759?
Szukamy największej naturalnej potęgi dwójki, mniejszej od dzielnika: DZIELNIK > K = 2 ^ N.
W naszym przypadku K=8192, N = 13.
Obliczamy T= K/DZIELNIK * (2^32).
Czyli T=8'192/10'0002^32=0.8192 4'294'967'296= 3'518'437'208.8832.
Zaokrąglona wartość T jest właśnie magiczną stałą użytą w mnożeniu. O ile bitów przesunąć starsze słowo wyniku mówi nam N.
Należy sprawdzić, czy część ułamkowa wartości T jest większa od 0.5. Jeżeli tak nie jest to przed mnożeniem należy dzielną zwiększyć o jeden.
Na przykład dzielenie przez 100. K=64, N=6, T= 2'748'779'069.44, frac(T)==0.44 < 0.5, konieczna jest poprawka (instrukcja incl, przed mull).

Załączam program liczący 26'550! (silnia), PIII 996MHz, czas obliczeń:
1 - 29.66 sekundy
2 - 17.15
3 - 10.10
Fragment ten jest powtarzany ponad 310 milionów razy.

 uses sysutils;

const
      METODA=     3;
      radix=  10000;
      max=    33434;
      n=      26550;                 // policzymy N!
      log= round(ln(radix)/ln(10));  // log10(radix)
var
      t:array[0..max-1] of cardinal;
      tz,i,j,c,r:cardinal;
      kk,k:integer;
      s:string;
      mc:qword;
      d:double;                      // czas
begin
      d:=time; mc:=0;

      write(n/1000:0:3,'! Radix ',
          radix/1000:0:3,' (',log,'), ');

      // no to liczymy
      t[0]:=1;
      k:=0;              // indeks najstarszej "cyfry"
      tz:=0;             // zera na koncu omijamy
      for i:=2 to n do begin
          c:=0;
          while t[tz]=0 do
              tz+=1;
          for j:=tz to k do begin  mc+=1;
              c+= t[j]*i;   //n*radix+radix < 2^32
              if METODA=1 then begin

                  t[j]:=c mod radix;
                  c:=c div radix;

              end else if METODA=2 then begin

                  r:=c div radix;
                  t[j]:=c - r*radix;
                  c:=r;

              end else begin

                  asm
                    movl  c,%eax
                    movl $-776530087,%edx
                    mull  %edx
                    shrl  $13, %edx
                    movl  %edx,r
                  end ['EAX','EDX'];
                  t[j]:=c - r*radix;
                  c:=r;

              end
          end;
          while c<>0 do begin
              k+=1;
              r:=c div radix;
              t[k]:=c - r*radix;
              c:=r
          end;
      end;
      // policzone
      d:=time-d;

      // pokazemy maksymalnie 300*log poczatkowych cyfr
      writeln;
      kk:=k-299;
      if kk<0 then kk:=0;
      //kk:=0;          //jezeli chcesz zobeczyc calosc
      write(t[k]:log, #32);
      for j:=k-1 downto kk do begin
          if (k-j) mod 10=0 then writeln;
          str(t[j]:log,s);
          for i:=1 to length(s) do
              if s[i]=#32 then write(0)
              else write(s[i]);
          write(#32);
      end;
      if kk>0 then write('...');
      writeln; writeln;

      // najstarsza "cyfra" może mieć mniej niż log cyfr
      c:=t[k]; j:=0;
      while c<>0 do begin
          c:=c div 10; j+=1;
      end;
      writeln(k*log+j, ' cyfr, ',
              mc/1e6:0:6, ' mln. mnozen, ',
              d*24*60*60:0:2, ' sekund,  ',k+1,' limb.');
      writeln
end.

W komórkach tablicy znajdują się cyfry wyniku, pod indeksem zero najmłodsza, pod indeksem K najstarsza.
Na początek wstawiamy do tablicy jedynkę (t[0]:=1; k:=1). Następnie mnożymy naszą superliczbę prze kolejne liczby, od 2 do N (pamiętając o ewentualnym przeniesieniu).
Silnia w prawie jednej trzeciej składa się z zer, warto je ominąć (zmienna TZ).
Dokładnie, zer tych będzie:

   q:=5; zera:=0;
   repeat
     zera:= zera + n div q;
     q:= q*5;
   until q > n;

Można w komórkach tablicy przechowywać jedną cyfrę, lecz można też 4 (LOG), wtedy liczymy jakby w układzie o podstawie równej 10'000 (RADIX). Dzięki temu wykonamy czterokrotnie mniej obliczeń.
Jak wielkie może być N?. Oczywiście musimy mieć miejsce w tablicy na wynik, 26'550! to 105'932 cyfr dziesiętnych. Potrzebujemy więc tablicy z miejscem na 21'187 supercyfr (MAX).
Liczbę cyfr dziesiętnych silni poda funkcja

function liczbacyfr(n: cardinal): cardinal;
begin
      result:=  trunc(1 + (ln(2*pi*n)/2 + n*ln(n) - n) / ln(10));
end;

Posługujemy się 32 bitowymi liczbami bez znaku, ważne jest by wartość "c:=c + t[j]*i" nie przekroczyła 2^32-1. Osiągnąć może ona wartość (z pewnym przeszacowaniem) równą RADIX * N + RADIX. Trzeba zdecydować się na pewien kompromis, albo zwiększamy RADIX (np. zmiana na 100'000, zmniejsza czas z 10.10 do 8.65), a wtedy maksymalna silnia jaką możemy obliczyć się zmniejsza, albo odwrotnie (kosztem czasu).
-----------------------------------------------------------------------------------------------------------[ Xitami ]-