## Program that calculates factorials

Discuss software for the Apple 1/replica 1

### Program that calculates factorials

I ported this program from Altair 680 VTL-2 code. It calculates factorials beginning with 1 and reaches 328 before it halts with ">32767" error. Still, the factorial of 328 is over 680 digits long! I need to find a way to keep the array elements below the 32767 limit. Apparently VTL-2 did not represent negative numbers so the array elements could be as high as 65535. If anybody has any ideas, please let me know.

For those of you who don't remember, a factorial is a number multiplied by lesser numbers until zero is reached. For example, the factorial of 5 (expressed as 5!) would be 5*4*3*2*1=120 or 6!=6*5*4*3*2*1=720 and so on....

Here's the listing:

1 REM CALCULATION OF FACTORIALS
2 REM PORTED FROM ALTAIR VTL-2
3 REM BY MIKE PENZO DEC 9, 2007
4 DIM A(500)
5 FOR I=1 TO 500
6 LET A(I)=0
8 NEXT I
10 LET P=1
20 LET L=2
30 LET A(1)=1
40 LET I=2
50 LET A(I)=0
60 LET I=I+1
70 IF L>I THEN 50
80 PRINT
90 PRINT
100 PRINT P;
110 PRINT "!="
120 REM
130 LET I=L+1
140 LET I=I-1
150 IF A(I)=0 THEN 140
160 PRINT A(I);
170 LET I=I-1
180 IF I=0 THEN 220
190 LET V=A(I)/10
200 PRINT V;
205 PRINT A(I)-(V*10);
210 GOTO 170
220 LET P=P+1
230 LET I=1
240 LET C=0
250 LET X=A(I)
260 LET A(I)=P*X
270 IF A(I)<X THEN 320
280 LET A(I)=A(I)+C
290 LET C=A(I)/100
300 LET Q=C
305 LET A(I)=A(I)-(Q*100)
310 LET I=I+1
320 IF L>I THEN 250
330 IF C=0 THEN 80
340 LET L=L+1
350 REM
360 LET A(I)=C
370 GOTO 290
380 END
mike_p

Posts: 14
Joined: Dec Sat 08, 2007 9:51 pm

This is really quite impressive.

Do you have a reference to the actual algorithm being employed? It will not be easy to reverse engineer the BASIC code and work out how to overcome the 32K limitation.

Ken
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 6:42 pm
Location: Sydney

OK - I've worked out what's going on here.

Try this variant:

Code: Select all
1 REM CALCULATION OF FACTORIALS
2 REM PORTED FROM ALTAIR VTL-2
3 REM BY MIKE PENZO DEC 9, 2007
4 DIM A(500)
5 FOR I=1 TO 500
6 LET A(I)=0
8 NEXT I
10 LET P=1
20 LET L=2
30 LET A(1)=1
40 LET I=2
50 LET A(I)=0
60 LET I=I+1
70 IF L>I THEN 50
80 PRINT
90 PRINT
100 PRINT P;
110 PRINT "!="
120 REM
130 LET I=L+1
140 LET I=I-1
150 IF A(I)=0 THEN 140
160 PRINT A(I);
170 LET I=I-1
180 IF I=0 THEN 220
190 LET V=A(I)/10
200 PRINT V;
205 PRINT A(I)-(V*10);
210 GOTO 170
220 LET P=P+1
222 REM P = Z * 100 + Y
225 Z = P / 100
228 Y = P - Z * 100
230 LET I=1
240 LET C=0
250 LET X=A(I)
260 REM X = H * 100 + R
265 H = X / 100
270 R = X - H * 100
275 REM A(I) = P * (H * 100 + R) + C
277 LET Q=C
280 C = P * H + Z * R + (Y * R + C) / 100
290 A(I) = (P * H - C + R * Z) * 100 + Y * R + Q
310 LET I=I+1
320 IF L>I THEN 250
330 IF C=0 THEN 80
340 LET L=L+1
350 REM
360 LET A(I)=C
363 LET C=A(I)/100
366 LET A(I)=A(I)-(C*100)
370 GOTO 310
380 END

Basically the problem arose because of a P*X multiplication going over 32K. But the result is subsequently divided by 100 (since the a[i]'s actually hold the values mod 100), so by a little re-ordering of the calculation the problem can be avoided.

I am testing the code above now and have just calculated 175! correctly. Will watch and see how far it gets. I think the only limitation with this version will be available RAM (and hence the DIM A() statement)

I also tested the code on a modern computer and correctly generated the 9991 digit factorial of 3246 without needing an intermediate value over 32767.

Ken
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 6:42 pm
Location: Sydney

Well it's made it to factorial 340 without problem now. (I verified the result with Mathematica)

I think the only limits now are P < 32767 and maybe the 960 characters to display on the screen

Code: Select all
340!=
5100864472103711080930193283927293363034
4671982885376039468582437981130708498852
2921313637295464989902035101176105977880
6133788936493577434500361683415188982472
4013463388848611414227407131388805908338
1129369427136543121020132372337085743889
5402320411929506183969598863382952818389
3692358325009116412814178085996438696887
2698869902631260054113918249061159323947
6045067552375315738712391517997520501623
8822164264738634685875388180608752182405
0410520570694749236518810642543186264101
0474387420418554164289649410408629854298
7877862041808504453671655093398052608362
3788608805851112806177131562232553426327
6929061018352747735327201401241600000000
0000000000000000000000000000000000000000
00000000000000000000000000000000000

Ken
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 6:42 pm
Location: Sydney

Kallikak wrote:Well it's made it to factorial 340 without problem now. (I verified the result with Mathematica)

I think the only limits now are P < 32767 and maybe the 960 characters to display on the screen

I ran it under an Apple II emulator (for acceleration purposes) and it ran fine up to 449!, where it died with a "RANGE ERR" when the variable "I" hit 501. It can run further by making the array A() larger. With DIM A(1000), it ran up to 807!.

Tom
cowgod

Posts: 31
Joined: Sep Fri 14, 2007 9:38 pm
Location: New Jersey

Ken,
Thanks for the workaround! I kept trying to single step through the code and identify the problem but the large number of iterations frustrated me quickly.

The original program was limited by the small amount of Altair 680 memory (usually 1K). I'm looking forward to running your mod on my Replica with the largest array 32K of memory will allow.
Mike
mike_p

Posts: 14
Joined: Dec Sat 08, 2007 9:51 pm

I set LOMEM = 512, HIMEM = 32767 and dimensioned the array for 15,500 elements. It's off and running, I'll let you know how high it gets!
mike_p

Posts: 14
Joined: Dec Sat 08, 2007 9:51 pm

I hope you are patient.

You could rem out the value printing code to speed things up a great deal.

Ken
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 6:42 pm
Location: Sydney

Holy toledo! After 25 hours and 20 minutes of processing time, my Replica is at 1429! and the result is one huge number. This thing might process for quite some time before erroring out. It should be interesting to determine the highest factorial it will calculate given 32K of RAM. This Replica is quite the supercomputer!

Ken, thanks again for the code fix.
mike_p

Posts: 14
Joined: Dec Sat 08, 2007 9:51 pm

I'm at 2249! and still running with no end in sight. Here's the final listing:

LOMEM = 512
HIMEM = 32767

1 REM CALCULATION OF FACTORIALS
2 REM PORTED FROM ALTAIR VTL-2
3 REM BY MIKE PENZO
4 REM AND KEN WESSEN DEC 2007
5 PRINT "INITIALIZING ARRAY, PLEASE WAIT"
6 DIM A(15500)
8 FOR I=1 TO 15500
10 LET A(I)=0
15 NEXT I
20 LET P=1
25 LET L=2
30 LET A(1)=1
40 LET I=2
50 LET A(I)=0
60 LET I=I+1
70 IF L>I THEN 50
80 PRINT
90 PRINT
100 PRINT P;
110 PRINT "!="
130 LET I=L+1
140 LET I=I-1
150 IF A(I)=0 THEN 140
160 PRINT A(I);
170 LET I=I-1
180 IF I=0 THEN 220
190 LET V=A(I)/10
200 PRINT V;
205 PRINT A(I)-(V*10);
210 GOTO 170
220 LET P=P+1
225 LET Z=P/100
228 LET Y=P-Z*100
230 LET I=1
240 LET C=0
250 LET X=A(I)
265 LET H=X/100
270 LET R=X-H*100
277 LET Q=C
280 LET C=P*H+Z*R+(Y*R+C)/100
290 LET A(I)=(P*H-C+R*Z)*100+Y*R+Q
310 LET I=I+1
320 IF L>I THEN 250
330 IF C=0 THEN 80
340 LET L=L+1
360 LET A(I)=C
363 LET C=A(I)/100
366 LET A(I)=A(I)-(C*100)
370 GOTO 310
380 END
mike_p

Posts: 14
Joined: Dec Sat 08, 2007 9:51 pm

Factorial 4776 has 15500 digits exactly, so looks like you are only half way through!
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 6:42 pm
Location: Sydney

Well I finally stopped the run of this program so I could perform a ROM dump for Ken. My Replica was at 3715!. The result compared exactly to the result from a factorial calculator I found on the web . The number is 11,652 digits long so I won't post it here!
mike_p

Posts: 14
Joined: Dec Sat 08, 2007 9:51 pm

I wonder how many people have run a 2 week long calculation on an Apple 1 before...

Ken
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 6:42 pm
Location: Sydney

If I only had the time I would try to convert it to assembler and see how much of a difference it would make...

However the code is pretty cryptic, does anyone have a version with some comments on what is really going on ?

fsafstrom

Posts: 154
Joined: Dec Tue 26, 2006 2:57 pm
Location: San Antonio, Texas

Yeah - Basic is just awful for structure.

To do the conversion that made it able to run so high I didn't need to work out the algorithm in full detail, but it seems to me it is pretty much just doing what you would do if you were calculating each successive factorial from the last manually - i.e. multiplying it out digit by digit, but using base 100. The 'digits' are being stored in A, and C is the carry.

Ken
Last edited by Kallikak on Jan Thu 03, 2008 4:42 pm, edited 1 time in total.
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 6:42 pm
Location: Sydney

Next