+Vorticon Posted March 15, 2020 Share Posted March 15, 2020 Happy birthday Pi! 9 Quote Link to comment Share on other sites More sharing options...
+dhe Posted March 15, 2020 Share Posted March 15, 2020 Maybe next year we can have a contest of sorts. Well make sure everyone uses a standard algorithm and people can pick a language to calculate PI in? Quote Link to comment Share on other sites More sharing options...
twoodland Posted March 16, 2020 Share Posted March 16, 2020 I would like to see the max # of digits possible using the SAMS board! Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 16, 2020 Author Share Posted March 16, 2020 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.? 3 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 16, 2020 Author Share Posted March 16, 2020 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 5 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 17, 2020 Author Share Posted March 17, 2020 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. 2 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 17, 2020 Author Share Posted March 17, 2020 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! Quote Link to comment Share on other sites More sharing options...
+mizapf Posted March 17, 2020 Share Posted March 17, 2020 Compared to just calculating p = 22/7, this is not really encouraging. 2 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 17, 2020 Author Share Posted March 17, 2020 54 minutes ago, mizapf said: Compared to just calculating p = 22/7, this is not really encouraging. It's the journey my friend 2 Quote Link to comment Share on other sites More sharing options...
+TheBF Posted March 17, 2020 Share Posted March 17, 2020 1 hour ago, mizapf said: Compared to just calculating p = 22/7, this is not really encouraging. Or 335/113 Quote Link to comment Share on other sites More sharing options...
+dhe Posted March 18, 2020 Share Posted March 18, 2020 once we settle on an algorithm, I'll give a shot at doing it in fortran99. 2 Quote Link to comment Share on other sites More sharing options...
moulinaie Posted March 18, 2020 Share Posted March 18, 2020 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. 1 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 18, 2020 Author Share Posted March 18, 2020 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 Quote Link to comment Share on other sites More sharing options...
senior_falcon Posted March 18, 2020 Share Posted March 18, 2020 13 hours ago, TheBF said: Or 335/113 Should be 355/113 1 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 18, 2020 Author Share Posted March 18, 2020 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 1 Quote Link to comment Share on other sites More sharing options...
senior_falcon Posted March 18, 2020 Share Posted March 18, 2020 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. Quote Link to comment Share on other sites More sharing options...
moulinaie Posted March 18, 2020 Share Posted March 18, 2020 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. DISK1.zip 1 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 18, 2020 Author Share Posted March 18, 2020 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? Quote Link to comment Share on other sites More sharing options...
moulinaie Posted March 18, 2020 Share Posted March 18, 2020 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. 2 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 18, 2020 Author Share Posted March 18, 2020 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! Quote Link to comment Share on other sites More sharing options...
moulinaie Posted March 18, 2020 Share Posted March 18, 2020 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 ? Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 18, 2020 Author Share Posted March 18, 2020 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. Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 18, 2020 Author Share Posted March 18, 2020 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 1 Quote Link to comment Share on other sites More sharing options...
moulinaie Posted March 18, 2020 Share Posted March 18, 2020 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. 3 Quote Link to comment Share on other sites More sharing options...
+Vorticon Posted March 20, 2020 Author Share Posted March 20, 2020 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). 2 Quote Link to comment Share on other sites More sharing options...
Recommended Posts
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.