Silnia na przykładzie szybkiego mnożenia

sapero

Obliczanie silnii dowolonej liczby
Prezentuję algorytm szybkiego mnożenia długich liczb, a jako przykład użycia - wykożystałem obliczenie silnii. Algorytm narodził się w czasach ZX Spectrum+ - wtedy pisało się programy tylko w assemblerze i to tak, aby kod był jaknajszybszy.

Składnia programu - język basic - ostatni krzyk mody - ibasic. Pomijam jego zalety bo to nie dotyczy algorytmu.

na początku trzeba założyć dwie tabele o pewnej ilości bajtów
char - określa zmienną jednobajtową, MaxDigits - liczba znaków w tabelach. Dla 65536 znaków można na spokojnie obliczać silnię 10000
tabela A to wynik, a tabela B to coś w stylu reszty. Jeden bajt w tabeli to jedna cyfra 0-9

const MaxDigits=65536
char a[MaxDigits],b[MaxDigits]

następnym krokiem jest podanie liczby n! - tu będzie nią N
trzeba sprawdzić czy n=0 | 1
jeśli n=0 albo 1 - wynik =n
ewentualnie czyścimy tabelę A:RtlZeroMemory(&A,MaxDigits)
wpisujemy liczbę 2 na początek wyniku (zeby nie mnożyć przez 1 :) )

a[0]=2
LenA=1

LenA - zmienna przechowująca ilość cyfr wyniku. Jest potrzebna aby przyspieszyć obliczanie
czyścimy obie tabele

wchodzimy w pętlę mnożenia dla liczb 3->N

timestart=GetTickCount()
for BC=3 to N
  mul=BC /* mul będzie modyfikowane */

  RtlZeroMemory(&B,LenA) /* czyścimy 'resztę' */

  do                 /* mnożenie długiej liczby a[] przez mul */
    if (mul&1)    /* jeśli mul jest nieparzyste */
      AddAtoB() /* dodaj wynik do 'reszty' */
      mul--        /* zmniejsz mul o 1 */
    else
      RLCA()      /* pomnóż wynik przez 2 */
      mul/=2      /* podziel mul przez 2 */
    endif
  until mul<2    /* powtarzaj if N>1 */

  AddBtoA()    /* dodaj 'resztę' do wyniku */
next BC

time=GetTickCount()-timestart

zmienna time ma ilość ms wykonywania się obliczeń - podziej ją przez 1000, wyświetl wynik do dwóch miejsc po przecinku
print "Time:"+using("#.##",time/1000.0),"s"
no i teraz trzeba wyświetlić wynik lub wpisać do pliku. Zaczynamy od ostatniego bajtu (na pozycji LenA-1) kończąc na pozycji 0

for x=LenA - 1 to 0 step -1
  print (albo write do pliku) chr$(a[x]+48)
next x

a teraz funkcje:
Mnożenie wyniku przez 2

sub RLCA
	Carry  = 0
	for HL=0 to LenA-1
		a[HL] = a[HL]*2 + Carry
		Carry = (a[HL]>9)
		if Carry
			a[HL] = a[HL]-10
		endif
	next HL	
	a[HL]+=Carry
	LenA+=Carry

	return
endsub

dodanie wyniku do reszty

sub AddAtoB
	Carry=0	
	for HL=0 to LenA-1
		b[HL]= b[HL] + a[HL] + Carry
		Carry = (b[HL]>9)
		if Carry
			b[HL] = b[HL]-10
		endif
	next HL
	b[HL]+=Carry
	LenA+=Carry	
	return

endsub

dodanie reszty do wyniku

sub AddBtoA
	Carry  = 0
	for HL=0 to LenA-1
		a[HL]= a[HL] + b[HL] + Carry
		Carry = (a[HL]>9)
		if Carry
			a[HL] = a[HL]-10
		endif
	next HL
	a[HL]+=Carry
	LenA+=Carry	
	return
endsub

Carry to zmienna przechowująca przeniesienie; Carry=True jeśli ostatnia operacja dodawania/mnożenia dała wynik>9
HL i BC - zwykłe zmienne, ich nazwy wywodzą się z Z80 :)

gotowca można ściągnąć

for HL=0 to LenA-1 działa tak: zmienna HL przyjmuje wartość 0. Wykonują się kolejne rozkazy aż program napotka rozkaz NEXT HL, wtedy HL przyjmuje kolejną wartość (tutaj o 1 większą) i dopuki HL<(LenA-1) - pętla powraca do następnego rozkazu po FOR

Carry = (a[HL]>9) działa tak: jeśli wyrażenie po znaku równości będzie TRUE to zmienna CARRY przyjmie wartość 1 (TRUE), w przeciwnym razie przyjmie wartość 0 (FALSE)

a[HL] - jeśli HL=0 to chodzi o pierwszy bajt z tabeli A - HL jest jakby wskaźnikiem/pointerem do kolejnych bajtów w tabeli. Pierwszy bajt w obu tabelach jest bajtem (cyfrą wyniku lub reszty) najmniej znaczącym (jedności)

if Carry - działa taksamo jak if Carry=1

LenA+=Carry - jeśli wynik ma jedną cyfrę, załóżmy że chwilowy wynik=5 - pomnożenie cyfry 5 przez 2 da liczbę 10 (czyli 2 cyfry). i zmienna Carry zgodnie z założeniem będzie równa 1 [patrz linijka Carry = (a[HL]>9)]. Każdy rozkaz NEXT w basicu ZAWSZE najpierw zmienia zmienną pętlową a potem sprawdza warunek podany w FOR, więc po wyjściu z FOR-NEXT zmienna HL wskazuje na pierwszy wolny bajt w tabeli wyniku A (z lewej strony, tutaj na drógą cyfrę, której zresztą jeszcze nie ma). Pierwsza cyfra wyniku=0 ponieważ warunek (liczba>9) został spełniony - wpisano do tabeli 10 a potem odjęto 10 i ustawiono Carry dla następnego mnożenia przez 2 (popularne JEDEN DALEJ). Ale pętla przebiegała dla liczb 0-1 więc wykonała się tylko raz - teraz wystarczy dopisać jedynkę z lewej strony wyniku (b[HL]+=Carry) albo b[HL]=Carry i zwiększyć ważną liczbę cyfr wyniku o jeden: LenA+=Carry (jeśli Carry=0 to LenA nie zwiększy się)

czyli zmienna Carry pełni dwuznaczną rolę - ustawia się gdy dodawanie/mnożenie przekroczyło 9, a zarazem po wyjściu z pętli wiadomo czy trzeba dopisać cyfrę do wyniku

wyniki:
10000! - pentium4 2.4GHz - 91s, athlonXP 2.1GHz - 77s

9 komentarzy

kurka , utrudniłem se życie :P <lol>

Dla kompletnych alików, przedstawiam kod, w ktorym smai decydujecie do jakiego n silna ma byc obliczona, znaczy sie nadajecie taka możliwość programowi, poprzez dopisanie wiekszych wartość , liczby elementów naturalnych do tablicy.

Oto kod:

program Silnia_Rekurencyjna;

{$APPTYPE CONSOLE}

uses
SysUtils;

{tablica kolejnych liczb naturalnych od 1 do 15}

const A: Array[1..15] of Longint=(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);

var
i ,n : Integer; //okreslenie zmiennych
silnia : Int64; //dla ogromnych wartości silni zmiencie typ na REAL

begin
silnia:=1; //nadanie wartosci poczatkowej
Writeln('Podaj liczbe naturalna: ');
Readln(n); //pobranie zmiennej n

for i:=1 to n do //wzrastanie potegowania [obliczania silni]
silnia:= silnia * A[i];

Writeln('Silnia: ', silnia); //wyswietlenie wyniku
Readln(n);

end.

Kodzik prosty, przyjemny i nie ma żadnych problemów z 0! :P . Pozdro for all !

Chyba dalo sie to prosciej zrobic :)
Ok zrobilem wlasnna silnie ponizej kod
--Cut Here--
function Silnia(liczba : integer; i : integer; n : integer):integer;
begin
liczba := 1;
for i := 1 to n do
begin
liczba := liczba * i;
end;
result := liczba;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
d : integer;
begin
d := Silnia(1,1,StrToInt(Edit1.Text));
Memo1.Lines.Add(IntToStr(d))

end;
--Cut Here--
To caly kod :P oczywiscie zeby silnia dzialala trzeba wstawic na forme komponent memo, jeden edit i button :)

" (..) jeśli n=0 albo 1 - wynik =n (...) "

Otóż nie, bo jeśli dla 0!=1 i dla 1!=1...

Wielkie dzieki za maila;]

Byłbym wdzięzny za algorytm i jesli moglbys-> matys359@wp.pl

niestety ostatnio pisałem w C jakieś 6 lat temu na spectrumie, C tylko znam ze strony tłumaczenia z C na basica. Operacje na tablicach cyfr raczej są proste (tu jest tylko dodawannie, a mnożenie przez 2 to też dodawanie)

ale jeśli komuś ciężko się czyta taki kod, to mogę to przerobić na czysty algorytm.

A to że dział nie ten - też już wiem :) ale nie widzę możliwości przeniesienia tematu

Mam prosbe moglbys tez pokazac implementacje w C++?? Bo o delphi nie mam pojecia a chetnie kod bym przeczytal;]

to chyba nie ten dział; to pownno być w "z pogranicza"