Silnia na przykładzie szybkiego mnożenia

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