PDA

View Full Version : n! - very slowly.



danbaron
28-08-2011, 11:47
I had an old script to multiply big integers.

It uses the simplest and slowest method of multiplication, like you would do with a pencil and paper.

I decided to try to add a factorial function.

It works correctly, but, it is really slow.

On my computer it took 448 seconds to find 200!.

In function TBMain(), you specify the value, "n", which you want the factorial of.

Then, you run it.

If, for instance, you set, n=200, then, as the program runs, you will see in the console window, it count down from 200 as it calculates, i.e.,

200
199
198
197
.
.

After it reaches 2, it will display the value of n!.

I checked using Racket, and 200! has 375 digits.

The elapsed time to find 200! in DrRacket was, 0 milliseconds.

At least our answers were the same.


' code ----------------------------------------------------------------------------------------------
Uses "console"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%neg = -1
%pos = +1

%intsize = 500

Type genint
digit(%intsize) As Byte
sign As Integer
error As Byte
End Type

'--------------------------------------------------------------

Function TBMain()
Local s As String
Local t1, t2, tt As Integer
Local n As Integer

' Set "n", for n!.
n = 200

t1 = Timer
s = genintfactorial(n)
t2 = Timer
tt = t2 - t1
PrintL
PrintL "elapsed time = ", tt, " seconds"
PrintL
PrintL n, " factorial:"
Console_WriteLine(s)

End Function

'--------------------------------------------------------------

Function genintfactorial(n As Integer) As String
Local num1, num2, num3 As genint
Local s As String
Local i As Integer

s = "1"
genintsetvalue(num1,s)
For i = n To 2 Step -1
PrintL i
s = Format$(i)
genintsetvalue(num2,s)
genintmult(num1,num2,num3)
genintsetequal(num1,num3)
Next
Return genintgetstringrep(num3)
End Function

'--------------------------------------------------------------

Function genintsetequal(gn2 As genint, gn1 As genint)
'Set gn2 = gn1.
Local i As DWord
For i = 1 To %intsize
gn2.digit(i) = gn1.digit(i)
Next
End Function

'--------------------------------------------------------------

Function genintsetvalue(gn As genint, vs As String)
Local i, l As DWord
genintzero(gn)
l = Len(vs)
For i = 1 To l
gn.digit(l - i + 1) = Val(Mid$(vs, i, 1))
Next
End Function

'--------------------------------------------------------------

Function genintzero(gn As genint)
Local i As DWord
For i = 1 To %intsize
gn.digit(i) = 0
Next
End Function

'--------------------------------------------------------------

Function genintsetsign(gn As genint, s As Integer)
gn.sign = s
End Function

'--------------------------------------------------------------

Function genintgetstringrep(gn As genint) As String
Local i, j As DWord
Local vs As String

vs = ""

If gn.sign > -1 Then
vs = vs & "+"
Else
vs = vs & "-"
EndIf

i = 0

Do
If %intsize - i = 0 Then Return "0"
If gn.digit(%intsize - i) > 0 Then Exit Do
i += 1
Loop

For j = i + 1 To %intsize
vs = vs & Chr$(48 + gn.digit(%intsize - j + 1))
Next

Return vs
End Function

'--------------------------------------------------------------

Function genintgetsize(gn As genint) As DWord
Local i As DWord

i = 0

Do
If %intsize - i = 0 Then Return 0
If gn.digit(%intsize - i) > 0 Then Exit Do
i += 1
Loop

Return %intsize - i
End Function

'--------------------------------------------------------------

Function genintseterror(gn As genint)
gn.error = TRUE
End Function

'--------------------------------------------------------------

Function genintshiftleft(gn As genint, nplaces As DWord)
Local sz, i, j As DWord

sz = genintgetsize(gn)

For i = 1 To sz
j = sz - i + 1
gn.digit(j + nplaces) = gn.digit(j)
Next

For i = 1 To nplaces
gn.digit(i) = 0
Next

End Function

'--------------------------------------------------------------

Function genintadd(a As genint, b As genint, c As genint)
Local i As DWord
Local s, v, k As Byte
Local asize, bsize, maxsize
Local temp As genint

asize = genintgetsize(a)
bsize = genintgetsize(b)
maxsize = asize
If bsize > maxsize Then maxsize = bsize

genintzero(temp)

k = 0
For i = 1 To maxsize + 1
s = a.digit(i) + b.digit(i) + k
v = Mod(s, 10)
k = s \ 10
temp.digit(i) = v
Next

For i = 1 To %intsize
c.digit(i) = temp.digit(i)
Next

End Function

'--------------------------------------------------------------

Function genintmult(a As genint, b As genint, c As genint)
Local i, j As DWord
Local p, v, k As Byte
Local asize, bsize, tsize As DWord
Local temp As genint

genintzero(c)
genintzero(temp)

asize = genintgetsize(a)
bsize = genintgetsize(b)
tsize = asize + bsize
If tsize > %intsize Then
genintseterror(c)
Return
EndIf

For i = 1 To asize
k = 0
For j = 1 To bsize + 1
p = a.digit(i) * b.digit(j) + k
v = Mod(p, 10)
k = p \ 10
temp.digit(j) = v
Next

genintshiftleft(temp, i - 1)
genintadd(c, temp, c)
genintzero(temp)
Next

End Function

' output --------------------------------------------------------------------------------------------
.
.
5
4
3
2

elapsed time = 448 seconds

200 factorial:
+7886578673647905035523632139321850622951359776871732632947425332443594499634033
42920304284011984623904177212138919638830257642790242637105061926624952829931113
46285727076331723739698894392244562145166424025403329186413122742829485327752424
24075739032403212574055795686602260319041703240623517008587961789222227896237038
97374720000000000000000000000000000000000000000000000000

C:\Users\root\Desktop\NEW THIN>

ErosOlmi
28-08-2011, 12:33
thinBasic and Racket language cannot be compared in term of speed!

Why?

thinBasic is a just in time 100% interpreted language
Source code is continuously parsed and interpreted by thinBasic software
Execution speed will always be dependent by software
racket is a just in time compiler: http://docs.racket-lang.org/guide/performance.html
Source code is loaded, compiled into machine code, executed inside hardware
Once compiled execution speed is only hardware dependent
If you want to compare something:
you can to compare Racket with Oxygen (they both are just in time compilers) or with true compilers like PowerBasic or Freebasic.

The only way to compare thinbasic with a language that creates compiled applications (or JIT compilers like racket) is to create in thinBasic native (compiled) functions.
That's why thinBasic has modules in which script ideas can get full hardware acceleration by compiling them into native functionality.

Ciao
Eros

danbaron
28-08-2011, 13:17
I didn't mean that.

I meant that there are ways to multiply (FFT - Fast Fourier Transform) that are unbelievably faster than the way I did the multiplication.

I'm almost positive that Racket uses FFT.

There are lots of tricks for fast multiplication, division, square roots, etc., of big integers.

I think it took mathematicians many years to develop those techniques.

If you don't know those tricks, then what you make (what I made), can look terrible compared to what others have made that do use them.

So far, it's just a fact, that, for instance, Racket, and much much more so, Mathematica, are light years ahead of what I can dream of doing.

I wasn't saying one bad word about thinBasic.

I was saying that my creation is relatively pitiful.

It wouldn't have mattered if it was compiled in C.

Maybe if I can learn how to do FFT multiplication, things will improve.

Then, maybe I could try to implement arbitrary precision math functions, sine, cosine, exp, etc.

On the other hand, the most important thing is not speed, it is getting the correct answer.

I wouldn't show a script in thinBasic in order to indicate how bad thinBasic is.

Since I like to post on the thinBasic site, I guess that wouldn't make much sense.

I still like to post my programs, even when they are relatively not very good, I think, that's how you learn.

And, to me, something doesn't have to be great in order to show it, I don't really care that much what people think of what I post.

If they don't like it, they don't have to look at it.

How many people could make an interpreter or a compiler?, you don't need to be defending thinBasic to me.

Why would I write one line of code in thinBasic, if I didn't like it?

I don't know why.

ErosOlmi
28-08-2011, 13:30
Sorry Dan. :oops:

Maybe I overwhelmed because speed comparison is one of the common tasks people do without knowing what compilers and interpreters are (not your case).
I'm quite badly sensible on this area.


Anyway, just to demonstrate how bad I'm sensible on speed comparison, I wrote down a thinBasic module compiled with Power Basic 9.05 using your source code.
With a module, getting 200! in any case takes 18 seconds, so much slower than native racket.

In case someone would like to have a look, attacked full module source and script example.

Ciao and sorry again
Eros



Calculating factorial of 200


elapsed time = 18 seconds


200 factorial:
+7886578673647905035523632139321850622951359776871732632947425332443594499634033
42920304284011984623904177212138919638830257642790242637105061926624952829931113
46285727076331723739698894392244562145166424025403329186413122742829485327752424
24075739032403212574055795686602260319041703240623517008587961789222227896237038
97374720000000000000000000000000000000000000000000000000

Charles Pegge
28-08-2011, 13:37
Here is a Binary-Coded-Decimal multiplier I posted a long time ago. One of your maths challenges Dan. If you like, add a factorial function to it and see how it performs :)

Charles

danbaron
28-08-2011, 14:01
18 seconds seems pretty good to me.

But, that's what I meant.

Somehow, they know how to do it really fast.

Why shouldn't we know too?

------------------------------------------------------

I downloaded yours, Charles.

But, I don't know how to do anything with Oxygen.

Maybe I could call the Oxygen multiply function from the thinBasic factorial function I make?

Or, I guess it wouldn't kill me to read the Oxygen documentation.

Charles Pegge
28-08-2011, 15:06
Here it is with the factorial iteration at the bottom. It won't keep you waiting for long :)



'
'===================
'MONSTROUS FACTORIAL
'===================

Uses "oxygen"
Dim src As String
src="


basic


'------------------------------------------------------------------
function multiply(byval ia as string, byval ib as string) as string
'==================================================================


dim as string a,b,c,d
dim as long pa,pb,pc,pd,la,lb,lc,ld
dim as long nd,sh,qa


a=ia
b=ib
la=len a
lb=len b
lc=la+lb'+10
ld=lc'+20
c=nuls lc 'LINE ACCUMULATOR
d=nuls ld 'BLOCK ACCUMULATOR
pa=*a
pb=*b
pc=*c
pd=*d




pushad






'SETUP POINTERS
'==============


mov esi,pa : add esi,la
mov edi,pb : add edi,lb
mov edx,pc : add edx,lc
mov ebx,pa




mov qa,esi 'RIGHT START POSITION FOR NUMBER A
mov nd,edi 'SETUP NEXT DIGIT POINTER (B NUMBER)
mov sh,edx 'SETUP POSITION SHIFT POINTER






'CONVERT FROM ASCII TO BINARY CODED DECIMAL
'==========================================




mov edi,pa
mov ecx,la
(
dec ecx
jl exit
sub byte [edi],48
inc edi
repeat
)
mov edi,pb
mov ecx,lb
(
dec ecx : jl exit
sub byte [edi],48
inc edi
repeat
)






nextline:
'========


'MULTIPLY BY ONE DIGIT
'WORKING FROM RIGHT TO LEFT


dec edi
mov cl,[edi]
mov ch,0
(
dec esi
cmp esi,ebx : jl exit
mov al,[esi]
mov ah,0
mul cl
add al,ch 'ADD CARRY VALUE
mov ch,0 'CLEAR CARRY VALUE
(
cmp al,10
jl exit 'NO CARRY
mov ch,10 'DIVISOR
div ch '
mov ch,al 'CARRY VAL IN CH
mov al,ah 'REMAINDER NOW IN AL
)
dec edx
mov [edx],al
repeat
)
'FINAL CARRY
(
cmp ch,0
jz exit
dec edx
mov [edx],ch
)


'ADD TO BLOCK ACCUMULATOR
'========================


mov esi,pc : add esi,lc
mov edi,pd : add edi,ld
mov ah,0
mov ebx,pc




'BCD ADDITION
'
'WORKING FROM RIGHT TO LEFT


(
dec esi
cmp esi,ebx : jl exit
dec edi
mov al,0
xchg al,[esi] 'LOAD AND THEN CLEAR LINE DIGIT
mov cl,[edi]
add al,ah 'PREVIOUS CARRY
add al,cl 'OPERAND
(
mov ah,0
cmp al,10 : jl exit
sub al,10
inc ah
)
mov [edi],al
repeat
)




mov ebx,pa


mov esi,qa 'START POSITION FOR NUMBER A


mov edi,nd 'NEXT DIGIT IN NUMBER B
dec edi
mov nd,edi


cmp edi,pb : jle fwd done


'SHIFT OUTPUT TO LINE ACCUM


mov edx,sh
dec edx
mov sh,edx




jmp long nextline






done:






'CONVERT FROM BCD TO ASCII
'=========================




mov edi,pd
mov ecx,ld
add ecx,edi
(
cmp edi,ecx : jge exit
add byte [edi],48 : inc edi
repeat
)






'TRIM LEADING ZEROS
'==================




mov edi,pd
mov ecx,ld
add ecx,edi
(
cmp edi,ecx : jge exit
mov al,[edi]
inc edi
cmp al,48 : jg exit
repeat
)
sub edi,pd
mov nd,edi




popad


function=mid(d,nd,ld)


end function




'----
'MAIN
'====




string a=""
string b="1"


sys factorial=200


for i=1 to factorial
b=multiply(b,str(i))
next


print b




"


'MsgBox 0,O2_PREP src ': Stop
O2_ASMO src
If Len(O2_ERROR) Then
MsgBox 0, O2_ERROR
Else
O2_EXEC
End If

ErosOlmi
28-08-2011, 16:12
Charles,

your O2 code on my computer it takes 0.002 seconds (two milliseconds) for 200! and 0.063 seconds for 1000!

Ciao
Eros

Charles Pegge
28-08-2011, 17:24
Better than I expected Eros. No time for a cup of tea. Long Division is going to be hard work though.

Charles

danbaron
28-08-2011, 22:06
Wow, that's good!

And the entire code is not very long.

It makes me think that what Racket and Mathematica are doing is less like magic.

Charles Pegge
29-08-2011, 00:12
I have now exposed the multiply and factorial functions so they can be called from thinBasic. . I have later versions of Oxygen but this program works with the o2 included with thinBasic V1.8.8.0.

Charles

danbaron
29-08-2011, 06:09
For 10,000!,

I got,

35,660 digits,

in,

36.182 seconds.

(No complaints at all.)