Jump to content
IGNORED

Pi Day - Calculate Pi on the TI


Vorticon

Recommended Posts

2 hours ago, twoodland said:

I would like to see the max # of digits possible using the SAMS board!

 

That would not be trivial given that the SAMS memory is accessible only in 4K chunks. Theoretically, you should be able to calculate over 166,000 decimals with it using Picalc's algorithm.?

  • Like 3
Link to comment
Share on other sites

Here's an interesting way to approximate Pi taken from A.K. Dewdney's book "The Armchair Universe" from 1988:

Assume a square field with a circle within it which touches all 4 sides of the square representing a pond in a cartesian field. The origin of the field is (0,0) and is the center of the pond, with (1,1) being the northeastern corner of the field. Fire a cannon randomly at the field n times, count the number of shots that land in the pond np, and calculate API=(n/np)*4. The result should be a very rough approximation of Pi. To get anywhere close though, n should be at least 100,000. However, the number of iterations required to get successive decimals of pi grows exponentially, and you probably will need something around 4 billion tries to get the first 4 decimals right!

Here's the very short program:

 

10 CALL CLEAR
20 RANDOMIZE
30 INPUT "NUMBER OF TRIES? ":N
40 FOR I=1 TO N
50 X=RND::Y=RND
60 D=SQR(X^2+Y^2)
70 IF D<=1 THEN TOTAL=TOTAL+1
80 PRINT "ITERATION # ";I
90 NEXT I
100 API=(TOTAL/N)*4
110 PRINT "APPROXIMATE PI=";API
120 END

I set up the program to run 200000 iterations on my TI before I left for work today and we'll see what I have when I get back :)

  • Like 5
Link to comment
Share on other sites

So it took about 14.5 hrs to complete 200,000 iterations, and I got Pi = 3.17228, whereas Dewdney got 3.14137 with 134,000 iterations. Hmmmm... My result is much worse than his even with a larger number of iterations, and I wonder if the randomizing function is to blame here... There is an important fact to mention here though: Dewdney had asked readers to send in runs of 1000 and he got 134 runs returned, and he averaged the combined result. Most likely each reader used a different machine with a different randomizing function although I can't imagine it should make much of a difference...

I'm going to do another run to match Dewdney's 134,000 iterations and see what I get.

  • Like 2
Link to comment
Share on other sites

Pi = 3.17322 with 134000 iterations.... So it appears that a larger number of iterations is going to be needed :) How about 1 million?

By the way, if anyone wants to give it a shot as well, it is best to delete line 80 because printing to the screen slows things down significantly.

More to come!

Link to comment
Share on other sites

On 3/16/2020 at 12:50 PM, Vorticon said:

 



50 X=RND::Y=RND
60 D=SQR(X^2+Y^2)
70 IF D<=1 THEN TOTAL=TOTAL+1

 

You can speed up things here.

1) SQR is not needed because if SQR(x)<1 then it means that x<1.

2) you can directly write D=RND^2+RND^2

3) or even better replace those three lines with:

50 IF RND^2+RND^2<=1 THEN TOTAL=TOTAL+1

 

I'm not sure if a condition returns -1 when true... Else, everything can be packed as:

 

50 TOTAL=TOTAL-( RND^2 + RND ^2 <=1)

 

Tell me if it runs faster..

 

Guillaume.

  • Like 1
Link to comment
Share on other sites

1 hour ago, moulinaie said:

You can speed up things here.

1) SQR is not needed because if SQR(x)<1 then it means that x<1.

2) you can directly write D=RND^2+RND^2

3) or even better replace those three lines with:

50 IF RND^2+RND^2<=1 THEN TOTAL=TOTAL+1

 

I'm not sure if a condition returns -1 when true... Else, everything can be packed as:

 

50 TOTAL=TOTAL-( RND^2 + RND ^2 <=1)

 

Tell me if it runs faster..

 

Guillaume.

Very nice Guillaume! While there is no execution speed difference between the packed and unpacked statements, the original listing is slower by 15 seconds over 200 iterations, which translates into about 4.2hrs over 200,000 iterations!

I'm going to use this for a 1 million iterations trial :)

Link to comment
Share on other sites

So that's interesting. In my previous trials, I had been using RXB. On a hunch, I tried it with standard XB for 200,000 iterations and I got Pi = 3.1482, a much better result. Is the RND function of RXB different from XB?

 

Also regarding Picalc, I dug up my documentation on it and I used Machin's formula which is Pi = 16 * arctan(1/5) - 4 * arctan(1/239)

 

Fred Kaal did some speed comparisons for 100 decimals in August 2007 which he posted on the old Yahoo listserve:

Geneve: ExBasic = 7min 32sec, C99 = 30 seconds

TI 99/4A: ExBasic = 13min 51sec, C99 = 1min 16sec

I'm seriously thinking about creating versions for Pascal and assembly just for kicks :)

 

In case you want to look over the XB or C99 code, here's the disk which is also found on the tigameshelf.net site under Edutainment

 

 

Picalc.dsk

  • Like 1
Link to comment
Share on other sites

So that's interesting. In my previous trials, I had been using RXB. On a hunch, I tried it with standard XB for 200,000 iterations and I got Pi = 3.1482, a much better result. Is the RND function of RXB different from XB?

 

Yes, RXB and XB G.E.M. use the random number generator from TI BASIC. This is more than 5x faster than the XB random number generator, but it appears the more complicated random number generator was not crazy after all, but was used to fix some deficiencies in the older TI BASIC version. It might be possible to use a flag to tell XB which RND generator to use. i.e. CALL LOAD(-1,1) to use the slower XB version or CALL LOAD(-1,0) to use the fast one.

 

Link to comment
Share on other sites

Hello,

 

I have made a version of the PI ransom calculation with MLC. Here is how it looks:

 

A square table is calculated during 30 seconds.

Then, it only requires 15 seconds to test 25.000 points into the square.

 

So, calculating up to 1.000.000 points will only require 10min30s (more or less).

 

The MLC Source Code:

; this program computes an approximation of PI
; it gets a random point into a square of vertice 1
; and counts how many belong to the quarter of circle inside the square
; the ratio (into circle) / (total points) is the same
; as the ratio of surfaces, so PI/4.


100 CALL CLEAR
$MLC F 110 10 3000
870 CALL CLEAR::RANDOMIZE::T=INT(RND*1000)
880 CALL CHAR(64,"007EA42424244400")::CALL CHAR(33,"0060920C60920C00")
890 PRINT "Building Square Table..."::PRINT "(about 30 sec...)"
900 CALL LINK("INITSQ",T)::PRINT
910 PRINT "Every !15sec, 25000 pts"::PRINT
920 T=0::C=0::ST=25000
1000 T=T+ST::PRINT "#Pts";T;
1010 CALL LINK("PIRAND",ST,M)::C=C+M
1020 PRINT "@!";4*C/T
1070 GOTO 1000
$INITSQ 0
	FARRAY 0
	DIMTABLE A 4096	; the square table
	FLOAT 512
	FMOVE 0 2		; r2 divisor = 512
	FLOAT -3
	FMOVE 0 3		; r3 increment = -3 (first dum value for loop)
	FLOAT 2
	FMOVE 0 4		; r4 increment of increment = 2
	FLOAT 1
	FMOVE 0 5		; r5 current square=1 (first dum value for loop)
	
	FOR I 0 4095
		FMOVE 3 0		; old increment
		FMOVE 4 1		; 2
		FADD			; new increment
		FMOVE 0 3		; saved
		FMOVE 5 1		; old square
		FADD			; new square
		FMOVE 0 5		; saved
		FMOVE 0 1		; for division r1=dividend
		FMOVE 2 0		; 512 and r0 divisor
		FDIV			; square/512
		INTEGER S		; to int
		PUTTABLE A(I) S	; one more square in the table
	NEXT
	GETPARAM 1 I		; get T value
	NDO I I				; and get T random values
		RND				; to start at a random position of the generator
	NLOOP
$$
$PIRAND 0
	GETPARAM 1 N		; total points
	CLEAR M				; reset the number of points into the 1/4 circle
	NDO N N
		RND					; random from 0000 to 0FFF
		GETTABLE A(Z) X		; its square reduced to 0000-7FFF
		RND					; another
		GETTABLE A(Z) Y		; another
		ADD X Y				; sum from 0000-FFFF
		IF>=				; if high bit not set, then...
			INC M			; one more into the 1/4 circle
		ENDIF
	NLOOP
	PUTPARAM 2 M			; returns how many into the 1/4 circle
$$
$END

 

You'll find attached a ZIP file with MLC and the PIRANDOM program for to be used with CLASSIC99.

 

Put the files into your DSK1 folder and type:

 

RUN DSK1.PIRANDOM

 

That's all !

 

Guillaume.

pirandom.png

DISK1.zip

  • Like 1
Link to comment
Share on other sites

36 minutes ago, Vorticon said:

While clearly you are using the same concept as the XB one, I'm not able to quite understand the algorithm which is impressively fast because of my unfamiliarity with MLC. 

Would you mind expanding on it a bit?

To avoid the calculation of the squares, I build a table with the squares of every integer from 0 to 4095 (because the RND function in MLC returns a value from 0 to 4095).

 

When you look at the square list:

0  1  4  9  16  25

You can see that the offsets are the odd numbers

0   (+1)   1   (+3)   4   (+5)   9   (+7)   16   (+9)   25   and so on... So, no need to multiply !

 

The maximum is 4095^2= 16 769 025 , too large for an integer!

That's why I divide every result by 512, so the maximum is 32752.

 

Then, hard to explain rapidely MLC but I tried to work as close as the ROM does to limit the number of extra code around every float instruction.

The fact is that the basic operations are performed between FAC and ARG that I call register 0 and 1 in MLC.

So, I have to fill every time R0 and R1 with the arguments of the operation an get the result from R0  to store it into R2-R3-R4-R5.

 

That's the explanation for "INITSQ", Am I clear enough????

(not so sure...)

 

Guillaume.

 

 

  • Like 2
Link to comment
Share on other sites

44 minutes ago, moulinaie said:

To avoid the calculation of the squares, I build a table with the squares of every integer from 0 to 4095 (because the RND function in MLC returns a value from 0 to 4095).

 

When you look at the square list:

0  1  4  9  16  25

You can see that the offsets are the odd numbers

0   (+1)   1   (+3)   4   (+5)   9   (+7)   16   (+9)   25   and so on... So, no need to multiply !

 

The maximum is 4095^2= 16 769 025 , too large for an integer!

That's why I divide every result by 512, so the maximum is 32752.

 

Then, hard to explain rapidely MLC but I tried to work as close as the ROM does to limit the number of extra code around every float instruction.

The fact is that the basic operations are performed between FAC and ARG that I call register 0 and 1 in MLC.

So, I have to fill every time R0 and R1 with the arguments of the operation an get the result from R0  to store it into R2-R3-R4-R5.

 

That's the explanation for "INITSQ", Am I clear enough????

(not so sure...)

 

Guillaume.

 

 

Yes I get it. Nicely done! So MLC's RND function only returns integers 0-4095? That's an interesting choice I have not seen before. Why is that?

Unfortunately, your method cannot be converted to XB without some contortions given the differences in RND returns,  which will defeat the potential speed benefit...

Perhaps you should consider converting Picalc to MLC! 

Link to comment
Share on other sites

1 minute ago, Vorticon said:

Yes I get it. Nicely done! So MLC's RND function only returns integers 0-4095? That's an interesting choice I have not seen before. Why is that?

Unfortunately, your method cannot be converted to XB without some contortions given the differences in RND returns,  which will defeat the potential speed benefit...

Perhaps you should consider converting Picalc to MLC! 

Why this choice? I found a very simple method to generate a random number that only uses integers. One simple 16-bits addition for each call ! Then I remove some of the 16 bits that I think are not so "at random", so I keep 12-bits out of 16.

 

Then in the PIRAND subroutine, only integers are used. That's why this is so fast.

 

Have you got the source code of PICALC ?

 

Link to comment
Share on other sites

5 hours ago, senior_falcon said:

So that's interesting. In my previous trials, I had been using RXB. On a hunch, I tried it with standard XB for 200,000 iterations and I got Pi = 3.1482, a much better result. Is the RND function of RXB different from XB?

 

Yes, RXB and XB G.E.M. use the random number generator from TI BASIC. This is more than 5x faster than the XB random number generator, but it appears the more complicated random number generator was not crazy after all, but was used to fix some deficiencies in the older TI BASIC version. It might be possible to use a flag to tell XB which RND generator to use. i.e. CALL LOAD(-1,1) to use the slower XB version or CALL LOAD(-1,0) to use the fast one.

 

Interesting... So the TI Basic RND function is sub-optimal... Learn something every day.

Link to comment
Share on other sites

7 minutes ago, moulinaie said:

Why this choice? I found a very simple method to generate a random number that only uses integers. One simple 16-bits addition for each call ! Then I remove some of the 16 bits that I think are not so "at random", so I keep 12-bits out of 16.

 

Then in the PIRAND subroutine, only integers are used. That's why this is so fast.

 

Have you got the source code of PICALC ?

 

I would love to see that RND algorithm.

Here's Picalc's XB source:

 

10 REM  PI CALCULATOR
20 REM  BY WALID MAALOULI
30 REM  AUGUST 2007
40 REM
100 CALL CLEAR
110 OPTION BASE 1
120 DIM SUM(652),SUM1(652),TERM(652),TEMP(652)
130 DISPLAY AT(1,8):"Pi Calculator" :: DISPLAY AT(3,6):"By Walid Maalouli" :: DISPLAY AT(5,9):"August 2007"
140 DISPLAY AT(10,2)BEEP:"# of decimals (mult. of 8)" :: DISPLAY AT(11,2):"(Maximum of 5200 decimals)" :: ACCEPT AT(13,12)VALIDATE(DIGIT):D
150 IF D>5200 THEN 140
160 ITR=INT(D/1.4) :: D=INT(D/8)+2
170 SUM(1)=3 :: SUM(2)=20000000 :: TERM(1)=0 :: TERM(2)=20000000 :: S=0 :: DENOM1=3 :: TBASE=25 :: MULT=16
180 FOR N=1 TO ITR+1
190 IF FLAG=0 THEN DISPLAY AT(20,5):"Term 1 iteration #";N;" " ELSE DISPLAY AT(20,5):"Term 2 iteration #";N;" "
200 IF N=1 THEN 220
210 FOR I=1 TO D :: TERM(I)=TEMP(I) :: NEXT I
220 DENOM=TBASE :: REMAINDER=0
230 GOSUB 1050
240 FOR I=1 TO D :: TEMP(I)=TERM(I) :: NEXT I
250 DENOM=DENOM1 :: REMAINDER=0
260 GOSUB 1050
270 FOR I=1 TO D
280 IF S=0 THEN 300
290 SUM(I)=SUM(I)+MULT*TERM(I) :: GOTO 310
300 SUM(I)=SUM(I)-MULT*TERM(I)
310 NEXT I
320 IF S=0 THEN S=1 ELSE S=0
330 FOR I=D TO 2 STEP-1
340 IF SUM(I)>=100000000 THEN 350 ELSE IF SUM(I)<0 THEN 390 ELSE 420
350 QUOTIENT=INT(SUM(I)/100000000)
360 SUM(I)=SUM(I)-QUOTIENT*100000000
370 SUM(I-1)=SUM(I-1)+QUOTIENT
380 GOTO 420
390 QUOTIENT=INT(SUM(I)/100000000)+1
400 SUM(I)=SUM(I)-(QUOTIENT-1)*100000000
410 SUM(I-1)=SUM(I-1)+QUOTIENT-1
420 NEXT I
430 IF FLAG=2 THEN 680
440 DENOM1=DENOM1+2
450 NEXT N
460 IF FLAG=1 THEN 620
470 TBASE=57121 :: DENOM1=3 :: MULT=4
480 FOR I=1 TO D
490 SUM1(I)=SUM(I)
500 SUM(I)=0 :: TERM(I)=0
510 NEXT I
520 TERM(1)=4
530 DENOM=239 :: REMAINDER=0
540 GOSUB 1050
550 FOR I=1 TO D
560 SUM(I)=TERM(I) :: TERM(I)=0
570 NEXT I
580 DENOM=239 :: TERM(1)=1 :: REMAINDER=0
590 GOSUB 1050
600 FLAG=1 :: S=0
610 GOTO 180
620 PRINT :: PRINT "Finalizing calculations..." :: PRINT
630 FOR I=1 TO D
640 SUM1(I)=SUM1(I)-SUM(I) :: SUM(I)=SUM1(I)
650 NEXT I
660 FLAG=2
670 GOTO 330
680 CALL CLEAR :: DISPLAY AT(2,3)BEEP:"Calculations complete!"
690 DISPLAY AT(5,3):"Send results to:" :: DISPLAY AT(8,5):"1- Screen" :: DISPLAY AT(10,5):"2- Printer (PIO)" :: DISPLAY AT(12,5):"3- File"
700 CALL KEY(0,K,ST) :: IF ST=0 THEN 700
710 IF K<49 OR K>51 THEN 700
720 ON K-48 GOTO 730,860,950
730 CALL CLEAR :: PRINT "Pi=3."
740 FOR I=2 TO D-1
750 SUM$=STR$(SUM(I))
760 IF LEN(SUM$)=8 THEN 780
770 SUM$="0"&SUM$ :: GOTO 760
780 PRINT SUM$;" ";
790 L=L+1 :: IF L=64 THEN 800 ELSE 830
800 PRINT :: PRINT :: DISPLAY AT(24,1)BEEP:"Press any key to continue" :: L=0
810 CALL KEY(0,K,ST) :: IF ST=0 THEN 810
820 CALL CLEAR
830 NEXT I
840 PRINT :: PRINT :: DISPLAY AT(24,1)BEEP:"End. Press any key to exit"
850 CALL KEY(0,K,ST) :: IF ST=0 THEN 850 ELSE STOP
860 OPEN #1:"PIO"
870 PRINT #1:"Pi=3."
880 FOR I=2 TO D-1
890 SUM$=STR$(SUM(I))
900 IF LEN(SUM$)=8 THEN 920
910 SUM$="0"&SUM$ :: GOTO 900
920 PRINT #1:SUM$;" ";
930 NEXT I
940 CLOSE #1 :: END
950 DISPLAY AT(22,1)BEEP:"Enter path.filename:" :: ACCEPT AT(24,1):FL$
960 OPEN #1:FL$,OUTPUT,VARIABLE 80
970 PRINT #1:"Pi=3."
980 FOR I=2 TO D-1
990 SUM$=STR$(SUM(I))
1000 IF LEN(SUM$)=8 THEN 1020
1010 SUM$="0"&SUM$ :: GOTO 1000
1020 PRINT #1:SUM$;" ";
1030 NEXT I
1040 CLOSE #1 :: END
1050 REM  DIVIDE SUBROUTINE
1060 FOR I=1 TO D
1070 DIVIDEND=REMAINDER*100000000+TERM(I)
1080 TERM(I)=INT(DIVIDEND/DENOM)
1090 REMAINDER=DIVIDEND-TERM(I)*DENOM
1100 NEXT I
1110 RETURN

 

  • Like 1
Link to comment
Share on other sites

27 minutes ago, Vorticon said:

I would love to see that RND algorithm.

 

The idea is this  you have an array of  seventeen 16bits integers:

from r(0) to r(16):

1966,23,3,31415,2718,7,666,17881,48,57,28030,991,8,7412,69,884,11357

 

You start with two pointers : i=4 and j=16

 

to get a random number:

r(j) = r(j)+r(i) (you don't care about overflow, just keep the 16 lower bits)

RND = r(j) shifted 4 bits to the right

i=i-1 :: if i<0 then i=16

j=j-1 :: if j<0 then j=16

 

That's all.

 

 

27 minutes ago, Vorticon said:

Here's Picalc's XB source:

 


10 REM  PI CALCULATOR
2

 

I will study it slowly.

I made an assembly adaptation of SUPER PI for the Motorola MC68030 (that runs on the Atari Falcon/TT).

It was easy because of the 32 bits of integers and the ability to multiply 32*32 = 64 bits.

The challenge will be different with a simple TI...

 

Guillaume.

  • Like 3
Link to comment
Share on other sites

I let the MLC version run for a while, and it seems that the Pi value is getting stuck around the 1.427 value +/- 0.005 or so after over 119 million iterations. I wonder if this is a limitation of the RND algorithm.

The 1 million iterations trial with the XB version is still running, and I'll be curious what kind of result I will get (should be complete within the next 48hrs).

 

 

 

PI.jpg

  • Like 2
Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...
  • Recently Browsing   0 members

    • No registered users viewing this page.
×
×
  • Create New...