## 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 10: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 7:42 pm
Location: Sydney

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

Try this variant:

Code: Select all
`1 REM CALCULATION OF FACTORIALS2 REM PORTED FROM ALTAIR VTL-23 REM BY MIKE PENZO DEC 9, 20074 DIM A(500)5 FOR I=1 TO 5006 LET A(I)=08 NEXT I10 LET P=120 LET L=230 LET A(1)=140 LET I=250 LET A(I)=060 LET I=I+170 IF L>I THEN 5080 PRINT90 PRINT100 PRINT P;110 PRINT "!="120 REM130 LET I=L+1140 LET I=I-1150 IF A(I)=0 THEN 140160 PRINT A(I);170 LET I=I-1180 IF I=0 THEN 220190 LET V=A(I)/10200 PRINT V;205 PRINT A(I)-(V*10);210 GOTO 170220 LET P=P+1222 REM P = Z * 100 + Y225 Z = P / 100228 Y = P - Z * 100230 LET I=1240 LET C=0250 LET X=A(I)260 REM X = H * 100 + R265 H = X / 100270 R = X - H * 100275 REM A(I) = P * (H * 100 + R) + C277 LET Q=C280 C = P * H + Z * R + (Y * R + C) / 100290 A(I) = (P * H - C + R * Z) * 100 + Y * R + Q310 LET I=I+1320 IF L>I THEN 250330 IF C=0 THEN 80340 LET L=L+1350 REM360 LET A(I)=C363 LET C=A(I)/100366 LET A(I)=A(I)-(C*100)370 GOTO 310380 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 7: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!=5100864472103711080930193283927293363034467198288537603946858243798113070849885229213136372954649899020351011761059778806133788936493577434500361683415188982472401346338884861141422740713138880590833811293694271365431210201323723370857438895402320411929506183969598863382952818389369235832500911641281417808599643869688726988699026312600541139182490611593239476045067552375315738712391517997520501623882216426473863468587538818060875218240504105205706947492365188106425431862641010474387420418554164289649410408629854298787786204180850445367165509339805260836237886088058511128061771315622325534263276929061018352747735327201401241600000000000000000000000000000000000000000000000000000000000000000000000000000000000`

Ken
Kallikak

Posts: 172
Joined: Jan Sun 29, 2006 7: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 10: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 10: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 10: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 7: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 10: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 10: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 7: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 10: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 7: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 3: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 5:42 pm, edited 1 time in total.
Kallikak

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

Next 