Silnia powyżej 13 - procesor 32b

Siema, widziałem podobne tematy, ale nie wiem, czy mnie satysfakcjonują. W asmie jestem dość świeży.
To mój kod(pisany w środowisku Delphi, dlatego tak wygląda) - na początku w EAX mam liczbę, dla której liczymy silnie:

mov ECX, EAX  {bo w EAX pojawia się wartość X, od teraz będzie w ECX}
mov EAX, 1    {Ustawiam EAX na 1 - pierwszy czynnik}
mov EBX, 1    {i drugi czynnik - licznik pęrli}
xor EDX, EDX

@loop:
  mul EBX     {EDX:EAX = EAX * EBX}
  inc EBX     {EBX = EBX + 1}
  cmp EBX, ECX  {czy EBX = X}
  jle @loop    {if EBX<=x goto loop}

OK i niby wszystko jest ok, ale do liczby 13 włącznie. A to dlatego, że kolejne liczby zajmują już więcej niż 4 bajty. I problem polega na tym, że mul EBX mnoży EBX przez EAX, a powinien przez EDX:EAX. Czy da się to jakoś zrobić?

Zaznaczam, że procesor jest 32 bitowy, tak więc nie posiada 64 bitowych rejestrów

Musisz wykonać kilka mnożeń.

Domyślam się ;) Możesz rozwinąć albo rzucić jakimś hasłem wg którego mam szukać?

Hasło: własna arytmetyka

a wlasna arytmetyke w asm implementuje sie latwiej niz w jezykach wyzszego poziomu i dziala ona sprawnie.

o co mniej wiecej chodzi?
mnozenie liczb nb (naturalny binarny)

materialy:
wyklady profesora pwr
http://www.zak.ict.pwr.wroc.pl/materialy/architektura/wyklad%20AK1/AK1-3-11-%20Dodawanie%20i%20mnozenie.pdf
slajd MUL-7 to obrazuje
tworzysz iloczyny cząstkowe a nastepnie je dodajesz na odpowiednich pozycjach.
to wlasna arytmetyka

wynik w tablicy, długość tablicy (w słowach), a możliwość policzenia silni

rozm   n! max
1       12
2       20
3       27
4       34
5       40
6       46
7       51
8       57
9       62
10      67
11      73
12      78
13      83
14      88
15      93
16      98
17      102
18      107
19      112 

typedef unsigned int u32;

#define MUL(hi, lo) { \
	unsigned long long int a=(unsigned long long int)hi * lo; \
	hi = a>>32; lo = a; }

int main(){
	const int M=17; // wystarczy dla 100!
	u32 t[M], p, ax, dx, i, n;
	char * w, CY;
	n=100; //scanf("%d", %n);
	//*************************************	
	w=(char*)t; 
	i=M;
	do{
		*(u32*)w = 0; 
		w += 4;
		i--;
	} while (i>0); 
	t[0]=1;
	while( n != 0 ){
		i=M;
		p=0;
		w=(char*)t;
		do {
			ax = n;
			dx = *(u32*)w;
			MUL(dx, ax)
			ax += p;
			CY = (ax<p); // samo się zrobi :)
			dx = dx + 0 + CY; 
			*(u32*)w = ax;
			w += 4;
			p=dx;
			i--;
		} while( i != 0 );
		n--; }
	//*******************************************************
	int k;
	for( k=M-1; t[k]==0; k-- ) ;
	while( k>=0 )
		printf("%08x", t[k--]);
	puts("\n00001b30964ec395dc24069528d54bbda40d16e966ef9a70eb21b5b2943a321cdf
               "10391745570cca9420c6ecb3b72ed2ee8b02ea2735c61a000000000000000000000000");} 

dla porównania wynik wg WolframAlpha
powinno dać się przełożyć jeden do jeden http://ideone.com/MLEY6

test30 napisał(a):

a wlasna arytmetyke w asm implementuje sie latwiej niz w jezykach wyzszego poziomu i dziala ona sprawnie.

Jak to łatwiej :D ?

Np. w Python'ie rekurencyjnie zasuwa super szybko silnia np. z 5000 ;)

import sys
sys.setrecursionlimit(100000)

silnia = (lambda n:
	1 if n==0 else
	n*silnia(n-1))
	
print silnia(5000L)

Bez rekurencji kod też jest krótki, a działa zapewne szybciej

def silnia(n):
    s=1
    for i in range(2,n+1):
        s*=i
    return s
print silnia(5000)

a możesz sprawdzić jakie byłyby czasy dla Pytho'nowej wersji takiej silni?

f( a, b )={ 
	my(c);
	if( b == a, return(a));
	if( b-a > 1,
		c=(b + a) >> 1;
		return(f(a, c) * f(c+1, b))
	);
	return( a * b );
}
 
fact(n) = f(1, n) 

Musisz wykonać 8 mnożeń, poprzesuwać/porozdzielać bity i pododawać.
Dla przykładu, jeśli słowo maszynowe miałoby 4 bity i chcemy pomnożyć 8-bitową liczbę przez 4-bitową:

(litery reprezentują bity)

AABBCCDD * EEFF = 
= (AA000000 + BB0000 + CC00 + DD) * (EE00 + FF) = 
=
    AA000000 * EE00 +
    BB0000 * EE00 +
    CC00 * EE00 +
    DD * EE00 +
    AA000000 * FF +
    BB0000 * FF +
    CC00 * FF +
    DD * FF
=
    AA * EE * 100000000 +    <--- to będzie zawsze obcięte więc pomijamy

    BB * EE * 1000000 +      <--- połowa bitów trafi do wyższego słowa wyniku, druga połowa będzie obcięta
    AA * FF * 1000000 +      

    CC * EE * 10000 +        <--- wszystkie bity trafią do wyższego słowa wyniku, bez przesuwania
    BB * FF * 10000 +

    DD * EE * 100 +          <--- połowa bitów leci do wyższego słowa wyniku, druga połowa do niższego
    CC * FF * 100 +

    DD * FF                  <--- wszystkie bity trafią do niższego słowa wyniku, bez przesuwania

Dodając kolejne składniki do niższego słowa wyniku musisz pamiętać aby również przenosić bit przepełnienia do wyższego słowa.

Przy większej ilości bitów będzie podobnie, zmieniają się tylko długości:

AAAAAAAABBBBBBBBCCCCCCCCDDDDDDDD * EEEEEEEEFFFFFFFF = 
    AAAAAAAA * EEEEEEEE * 100000000000000000000000000000000 +
    BBBBBBBB * EEEEEEEE * 1000000000000000000000000 +
    AAAAAAAA * FFFFFFFF * 1000000000000000000000000 +
    CCCCCCCC * EEEEEEEE * 10000000000000000 +
    BBBBBBBB * FFFFFFFF * 10000000000000000 +
    DDDDDDDD * EEEEEEEE * 100000000 +
    CCCCCCCC * FFFFFFFF * 100000000 +
    DDDDDDDD * FFFFFFFF

Dla 5000 między około 11 ms, dla 100'000 około 1,7 s.

znaczy się potrafi mnożyć

a tu wersja silni na dwa słowa (64 bity)

#define MUL(k) { \
        unsigned long long int a=(unsigned long long int)ax * k; \
        dx = a>>32; ax = a; }
#define PUSH(k) stos=k;
#define POP(k)  k=stos;

main(){
	unsigned int ax,dx,n,lo,hi,stos;
	
	n=20;
	hi=0; lo=1;
	while( n > 1 ){
		ax=n; 
		MUL(lo)
		lo = ax;
		PUSH(dx)
		ax=n;
		MUL(hi)
		POP(dx)
		hi = ax + dx;
		n--; }
	printf("%08x%08x", hi, lo);} 

OK, dzięki za odpowiedzi :)

Xitami napisał(a):

a możesz sprawdzić jakie byłyby czasy dla Pytho'nowej wersji takiej silni?

f( a, b )={ 
	my(c);
	if( b == a, return(a));
	if( b-a > 1,
		c=(b + a) >> 1;
		return(f(a, c) * f(c+1, b))
	);
	return( a * b );
}
 
fact(n) = f(1, n) 

0:00:00.002795

def silnia( a, b ):
	c=0
	if b == a: return a
	if b-a > 1:
		c=(b + a) >> 1;
		return (silnia(a, c) * silnia(c+1, b))
	return a * b

Dla porównania silnia(5000):
iteracyjna: 0:00:00.007711
rekurencyjna: 0:00:00.011260
z przesunięciami: 0:00:00.002795

edit --------------------

Dla porównania silnia(100 000):

iteracyjna: 0:00:03.986534

rekurencyjna: Segmentation fault :D

z przesunięciami: 0:00:00.288030