Jump to content
IGNORED

Pi Day - Calculate Pi on the TI


Vorticon

Recommended Posts

19 hours ago, Vorticon said:

Well it seems that the random shot method has its limits. After 14.3 billion iterations, PI has converged to 3.141076. I guess we're not going to get that 4th decimal after all. Interesting...

 

PI-MLC.jpg.0eaebe97fff293b74ff05c1a6eaaa5e9.jpg

Hi !

 

The RND function used has a precision of 1/4096 = 0.00024.

The best we can expect is a value in PI +/- 0.00024

3.1413 and 3.1418

 

Guillaume.

  • Like 2
Link to comment
Share on other sites

Hi,

 

Supposing that every value of the RND generator appears as many times as the others, here is a calculation of the "Best PI' you can expect according to the number of bits of the RND value.

 

Note: for MLC, the number of bits is 12.

 

50 CALL CLEAR
100 INPUT "# of bits in RND : ":B
105 S=0::R=2^B-1::R2=R*R
110 FOR I=0 TO R
120 S=S+INT(SQR(R2-I*I))
130 NEXT I
140 PRINT "Best PI with";B;"bits"
150 PRINT 4*S/(2^(2*B))
160 PRINT
170 GOTO 100


Take care of the calculation time. Use the CPU overdrive mode...

 

Guillaume.

  • Like 1
Link to comment
Share on other sites

5 hours ago, moulinaie said:

Hi !

 

The RND function used has a precision of 1/4096 = 0.00024.

The best we can expect is a value in PI +/- 0.00024

3.1413 and 3.1418

 

Guillaume.

Good to know! The problem is that even after 14 billion+ iterations, that 4th decimal was still stuck at 0, outside of the theoretical range...

I would assume then that the XB RND function has the full TI precision and should give a better result, although it would take months to get a decent number of iterations...

Link to comment
Share on other sites

1 hour ago, Vorticon said:

Good to know! The problem is that even after 14 billion+ iterations, that 4th decimal was still stuck at 0, outside of the theoretical range...

I would assume then that the XB RND function has the full TI precision and should give a better result, although it would take months to get a decent number of iterations...

Hi,

 

Sure ! You're right.

But MLC can do FLOAT operations. The lack is that I can see no call to the XB RND function from assembler.

Else, I could try to program a faster version with MLC.

 

Guillaume.

 

Link to comment
Share on other sites

17 minutes ago, moulinaie said:

Hi,

 

Sure ! You're right.

But MLC can do FLOAT operations. The lack is that I can see no call to the XB RND function from assembler.

Else, I could try to program a faster version with MLC.

 

Guillaume.

Well XB RND in GROM is here and it shows the Assembly CALLS used:

<0218> A295 42,03,23 X2SEED BYTE >42,>03,>23,>15,>00 * =   33521, X2 INITIAL VAL
       A298 15,00
<0219> A29A 43,02,3E X1SEED BYTE >43,>02,>3E,>2A,>17 * = 2624223, X1 INITIAL VAL
       A29D 2A,17
<0220>               ***********************************************************
<0221>               *           PSEUDO-RANDOM NUMBER GENERATOR                 
<0222>               *      X(N+1) = (A*X(N)+C) MOD M;  RND = X/M               

99/4 GPL-ASSEMBLER (Pass 3) correct                                   PAGE 0016 
EQUATES EXEC-359
<0223>               *    WHERE:                 X = X2 * 1E7 + X1              
<0224>               *                           A = A2 * 1E7 + A1              
<0225>               *                           C = C2 * 1E7 + C1              
<0226>               *                           M = 1E14                       
<0227>               * ASSUMPTIONS:                                             
<0228>               *  (1) All numbers are integers; fractional parts are      
<0229>               *      truncated                                           
<0230>               *  (2) If the variables listed below start in the ranges   
<0231>               *     specified. They will also end in the ranges specified
<0232>               *                                                          
<0233>               * CONSTANTS: 0 <= A2 < 5E6 ; 0 <= C2 < 1E7                 
<0234>               *            0 <= A1 < 5E6 ; 0 <= C1 < 1E7                 
<0235>               * VARIABLES: 0 <= X2 < 1E7 ; 0 <= T1 <= 1E14 ; 0 <= T2 < 1E
<0236>               *            0 <= X1 < 1E7 ; 0 <= T3 <= 1E14 ; 0 <= T4 < 1E
<0237>               *                                                          
<0238>               *        STACK USAGE:                                      
<0239>               *            CONSTANT REFS      CONTANT REFS    CONTANT REF
<0240>               * +---------+      IN/OUT            IN/OUT          IN/OUT
<0241>               * | STACK+4 | X2*A1(F)(H)       --    ----      --    ---- 
<0242>               * +---------+                                              
<0243>               * | STACK+3 |   T2 (C)(J)       --    ----      --    ---- 
<0244>               * +---------+                                              
<0245>               * | STACK+2 |   T1 (B)(D)   new X1   (E)(N)     --    ---- 
<0246>               * +---------+                                              
<0247>               * | STACK+1 |old X1(A)(G)       T3   (K)(L) new X2   (M)(P)
<0248>               * +---------+                                              
<0249>               ***********************************************************
<0250>               * COMPUTE NEW VALUE FOR X1, SAVE IT IN V@RNDX1
<0251>               *                             STACK
<0252>               *                               SREFS   FAC CONTENTS
<0253> A29F 35,00,05 NRND   MOVE 5,V@RNDX1,@FAC        FAC = X1
       A2A2 4A,A3,A5
<0254> A2A5 86,4F           CLR  @FAC5                 FAC = CLR
<0255> A2A7 87,50           DCLR @FAC6                 FAC = CLR
<0256> A2A9 0F,77           XML  VPUSH          (A)    FAC = X1
<0257> A2AB 31,00,08        MOVE 8,G@RNDA1,@ARG        ARG = A1
       A2AE 5C,A3,57
<0258> A2B1 0F,08           XML  FMUL                  FAC = X1*A1
<0259> A2B3 31,00,08        MOVE 8,G@RNDC1,@ARG        ARG = C1
       A2B6 5C,A3,67
<0260> A2B9 0F,06           XML  FADD               T1=FAC = X1*A1+C1
<0261> A2BB 0F,77           XML  VPUSH          (B)    FAC = T1
<0262> A2BD 31,00,08        MOVE 8,G@RNDEM,@ARG        ARG = 1/1E7
       A2C0 5C,A3,77
<0263> A2C3 0F,08           XML  FMUL                  FAC = T1/1E7
<0264> A2C5 06,00,22        CALL GRINT              T2=FAC = INT(T1/1E7)
<0265> A2C8 0F,77           XML  VPUSH          (C)    FAC = T2
<0266> A2CA 31,00,08        MOVE 8,G@RNDEP,@ARG        ARG = 1E7
       A2CD 5C,A3,6F
<0267> A2D0 0F,08           XML  FMUL                  FAC = T2*1E7
<0268> A2D2 A7,6E,00        DSUB 8,@VSPTR
       A2D5 08
<0269> A2D6 0F,0C           XML  SSUB           (D) X1=FAC = T1-T2*1E7
<0270> A2D8 35,00,05        MOVE 5,@FAC,V@RNDX1        FAC = X1 (new)
       A2DB A3,A5,4A
<0271> A2DE 0F,77           XML  VPUSH          (E)    FAC = X1
<0272>               * COMPUTE NEW VALUE FOR X2, SAVE IT IN V@RNDX2
<0273> A2E0 35,00,05        MOVE 5,V@RNDX2,@FAC        FAC = X2
       A2E3 4A,A3,A0
<0274> A2E6 86,4F           CLR  @FAC5                 FAC = CLR
<0275> A2E8 87,50           DCLR @FAC6                 FAC = CLR
<0276> A2EA 31,00,08        MOVE 8,G@RNDA1,@ARG        ARG = A1
       A2ED 5C,A3,57
<0277> A2F0 0F,08           XML  FMUL                  FAC = X2*A1

99/4 GPL-ASSEMBLER (Pass 3) correct                                   PAGE 0017 
EQUATES EXEC-359
<0278> A2F2 A3,6E,00        DADD 8,@VSPTR
       A2F5 08
<0279> A2F6 0F,77           XML  VPUSH          (F)    FAC = X2*A1
<0280> A2F8 A7,6E,00        DSUB 24,@VSPTR
       A2FB 18
<0281> A2FC 0F,78           XML  VPOP           (G)    FAC = X1
<0282> A2FE A3,6E,00        DADD 32,@VSPTR
       A301 20
<0283> A302 31,00,08        MOVE 8,G@RNDA2,@ARG        ARG = A2
       A305 5C,A3,4F
<0284> A308 0F,08           XML  FMUL                  FAC = X1*A2
<0285> A30A 0F,0B           XML  SADD           (H)    FAC = X2*A1+X1*A2
<0286> A30C 31,00,08        MOVE 8,G@RNDC2,@ARG        ARG = C2
       A30F 5C,A3,5F
<0287> A312 0F,06           XML  FADD                  FAC = X2*A1+X1*A2
<0288> A314 0F,0B           XML  SADD           (J) T3=FAC = X2*A1+X1*A2
<0289> A316 A7,6E,00        DSUB 16,@VSPTR
       A319 10
<0290> A31A 0F,77           XML  VPUSH          (K)    FAC = T3
<0291> A31C 31,00,08        MOVE 8,G@RNDEM,@ARG        ARG = 1/1E7
       A31F 5C,A3,77
<0292> A322 0F,08           XML  FMUL                  FAC = T3/1E7
<0293> A324 06,00,22        CALL GRINT              T4=FAC = INT(T3/1E7)
<0294> A327 31,00,08        MOVE 8,G@RNDEP,@ARG        ARG = 1E7
       A32A 5C,A3,6F
<0295> A32D 0F,08           XML  FMUL                  FAC = T4*1E7
<0296> A32F 0F,0C           XML  SSUB           (L) X2=FAC = T3-T4*1E7
<0297> A331 35,00,05        MOVE 5,@FAC,V@RNDX2        FAC = X2 (new)
       A334 A3,A0,4A
<0298>               * COMPUTE NEW VALUE FOR RND, LEAVE IT IN FAC
<0299> A337 31,00,08        MOVE 8,G@RNDEM,@ARG        ARG = 1/1E7
       A33A 5C,A3,77
<0300> A33D 0F,08           XML  FMUL                  FAC = X2/1E7
<0301> A33F 0F,77           XML  VPUSH          (M)    FAC = X2/1E7
<0302> A341 A3,6E,00        DADD 8,@VSPTR
       A344 08
<0303> A345 0F,78           XML  VPOP           (N)    FAC = X1
<0304> A347 0F,08           XML  FMUL                  FAC = X1/1E7
<0305> A349 0F,08           XML  FMUL                  FAC = X1/1E14
<0306> A34B 0F,0B           XML  SADD           (P)RND=FAC = (X2/1E7)+(X1/1E14)
<0307> A34D 0F,75           XML  CONT
[0008]               ***********************************************************
[0009]                      COPY 'DSK5.MYXB5-D'
<0001>               * CONSTANTS FOR THE RANDOM NUMBER ROUTINE
<0002> A34F 43,01,2B RNDA2  BYTE >43,>01,>2B,>59,>52,>00,>00,>00 * = 1438982
       A352 59,52,00
       A355 00,00
<0003> A357 42,2A,08 RNDA1  BYTE >42,>2A,>08,>15,>00,>00,>00,>00 * = 0420821
       A35A 15,00,00
       A35D 00,00
<0004> A35F 43,02,0B RNDC2  BYTE >43,>02,>0B,>20,>30,>00,>00,>00 * = 2113248
       A362 20,30,00
       A365 00,00
<0005> A367 43,06,36 RNDC1  BYTE >43,>06,>36,>05,>13,>00,>00,>00 * = 6540519
       A36A 05,13,00
       A36D 00,00
<0006> A36F 43,0A,00 RNDEP  BYTE >43,>0A,>00,>00,>00,>00,>00,>00 * = 1E7
       A372 00,00,00
       A375 00,00
<0007> A377 3C,0A,00 RNDEM  BYTE >3C,>0A,>00,>00,>00,>00,>00,>00 * = 1/1E7
       A37A 00,00,00
       A37D 00,00
<0008>               ***********************************************************
<0009>               *                   RANDOMIZE STATEMENT                    

99/4 GPL-ASSEMBLER (Pass 3) correct                                   PAGE 0018 
EQUATES EXEC-359
<0010>               ***********************************************************
<0011> A37F 06,6A,78 NRNDMZ CALL CHKEND            Seed provider?
<0012> A382 63,C0           BS   RNDM1             No
<0013>               * RANDOMIZE given a see value
<0014>               * (99,000,000,000,001 possible starting positions)
<0015>               * (Place-value is ignored in the input number)
<0016> A384 0F,74           XML  PARSE             Parse the seed
<0017> A386 83              BYTE TREMZ           * Up to end of statement
<0018> A387 06,A3,ED        CALL CKSTNM
<0019> A38A 8F,4A           DCZ  @FAC              Check FAC for zero
<0020> A38C 63,B6           BS   GA3B6
<0021> A38E BE,4A,46        ST   >46,@FAC          0 < FAC < 1E14
<0022> A391 0F,77           XML  VPUSH             Let FAC = X2*1E7+X1
<0023> A393 31,00,08        MOVE 8,G@RNDEM,@ARG        ARG = 1/1E7
       A396 5C,A3,77
<0024> A399 0F,08           XML  FMUL                  FAC = X2+X1/1E7
<0025> A39B 06,00,22        CALL GRINT                 FAC = X2
<0026> A39E 35,00,05        MOVE 5,@FAC,V@RNDX2        FAC = X2
       A3A1 A3,A0,4A
<0027> A3A4 31,00,08        MOVE 8,G@RNDEP,@ARG        ARG = 1E7
       A3A7 5C,A3,6F
<0028> A3AA 0F,08           XML  FMUL                  FAC = X2*1E7
<0029> A3AC 0F,0C           XML  SSUB                  FAC = X1
<0030> A3AE 35,00,05        MOVE 5,@FAC,V@RNDX1        FAC = X1
       A3B1 A3,A5,4A
<0031> A3B4 0F,75           XML  CONT                  FAC = X1
<0032> A3B6 BD,A3,A0 GA3B6  DST  @FAC,V@RNDX2          FAC = 0
       A3B9 4A
<0033> A3BA BD,A3,A5        DST  @FAC,V@RNDX1          FAC = 0
       A3BD 4A
<0034> A3BE 0F,75           XML  CONT
<0035>               * RANDOMIZE given number seed value (use GPL RAND function)
<0036>               * (16K possible starting positions)
<0037> A3C0 BF,4A,42 RNDM1  DST  >4201,@FAC            FAC = >4201
       A3C3 01
<0038> A3C4 86,4E           CLR  @FAC4                 FAC4= >00
<0039> A3C6 06,A3,D2        CALL RNDMZ
<0040> A3C9 03,A5           DATA RNDX1
<0041> A3CB 06,A3,D2        CALL RNDMZ             Set up seed
<0042> A3CE 03,A0           DATA RNDX2
<0043> A3D0 0F,75           XML  CONT              Continue on
<0044> A3D2 88,52    RNDMZ  FETCH @FAC8            Fetch address of seed (high b
<0045> A3D4 88,53           FETCH @FAC9            Fetch address of seed (low by
<0046> A3D6 02,63           RAND 99                GPL Randomize
<0047> A3D8 BC,4C,78        ST   @RANDOM,@FAC2     >00<=FAC+2<=FF
<0048> A3DB E6,4C,02        SRL  2,@FAC2           >00<=FAC+2<=3F
<0049> A3DE 02,63           RAND 99                GPL Randomize
<0050> A3E0 BC,4D,78        ST   @RANDOM,@FAC3     >00<=FAC+3<=FF
<0051> A3E3 E6,4D,02        SRL  2,@FAC3           >00<=FAC+3<=3F
<0052> A3E6 35,00,05        MOVE 5,@FAC,V*FAC8     Put in seed
       A3E9 B0,52,4A
<0053> A3EC 00              RTN
<0054> A3ED D6,4C,65 CKSTNM CEQ  >65,@FAC2
<0055> A3F0 6D,BE           BS   ERRSNM
<0056> A3F2 00              RTN
<0057> A3F3 40,01,00 FLT1   BYTE >40,>01,>00,>00,>00,>00,>00,>00
       A3F6 00,00,00
       A3F9 00,00
<0058>               ***********************************************************
17 minutes ago, moulinaie said:

 

 

Link to comment
Share on other sites

10 minutes ago, RXB said:

Well XB RND in GROM is here and it shows the Assembly CALLS used:



<0220>               ***********************************************************
<0221>               *           PSEUDO-RANDOM NUMBER GENERATOR                 
<0222>               *      X(N+1) = (A*X(N)+C) MOD M;  RND = X/M               

 

 

Ok thanks!

 

I got the method... I can simulate it in assembly.

 

But, isn't there a way to call this routine from assembler without rewriting it?

 

Guillaume.

  • Like 1
Link to comment
Share on other sites

2 hours ago, moulinaie said:

 

Ok thanks!

 

I got the method... I can simulate it in assembly.

 

But, isn't there a way to call this routine from assembler without rewriting it?

 

Guillaume.

Most of the work of making a Random Number is done is Assembly in XB ROMs.

 XML  VPUSH

 

 XML  VPOP

 

 XML  FMUL

 

 XML  FADD

 

 XML  SSUB
Link to comment
Share on other sites

Hi all,

 

I have been able to recreate with MLC (using its inline assembler) the Extended Basic RND routine.

 

snap.png.7ef8337890e145b9f305b0b98a68c71d.png

 

You can see on the left the value from MLC and on the right the XB RND call.

They start with the same seed (so, you must reboot your TI before), and the list generated is the same.

 

Then...

 

If you try form the command line

FOR I=1 TO 1000::A=RND::NEXT I

The XB function takes 72 seconds.

 

And if you try:

FOR I=1 TO 1000::CALL LINK("RANDOM",A)::NEXT I

MLC function takes 46 seconds.

And for the same result with the same precision !

 

Download the ZIP file, put everything in DISK1 folder and load DSK1.XBRND to test it.

DISK1.zip

 

Next work... Calculate PI with it!

 

Guillaume.

 

  • Like 3
Link to comment
Share on other sites

1 hour ago, moulinaie said:

They start with the same seed (so, you must reboot your TI before), and the list generated is the same.

 

Or...you can start each program with the same seed by using the same number (seed) with RANDOMIZE as the first statement.

10  RANDOMIZE 111

 

...lee

  • Like 1
Link to comment
Share on other sites

3 minutes ago, Lee Stewart said:

 

Or...you can start each program with the same seed by using the same number (seed) with RANDOMIZE as the first statement.


10  RANDOMIZE 111

 

...lee

Yes you're right!

I meant that with the same seed, the two lists are equal. That was to prove that MLC simulates exactly what XB does.

And... as MLC is not concerned by RANDOMIZE, but uses the boot seed of the TI, you must reboot the TI to be sure that the seeds ar the same.

 

I don't know if I'm clear... ?

 

Guillaume.

 

Link to comment
Share on other sites

2 hours ago, moulinaie said:

Hi all,

 

I have been able to recreate with MLC (using its inline assembler) the Extended Basic RND routine.

 

snap.png.7ef8337890e145b9f305b0b98a68c71d.png

 

You can see on the left the value from MLC and on the right the XB RND call.

They start with the same seed (so, you must reboot your TI before), and the list generated is the same.

 

Then...

 

If you try form the command line


FOR I=1 TO 1000::A=RND::NEXT I

The XB function takes 72 seconds.

 

And if you try:


FOR I=1 TO 1000::CALL LINK("RANDOM",A)::NEXT I

MLC function takes 46 seconds.

And for the same result with the same precision !

 

Download the ZIP file, put everything in DISK1 folder and load DSK1.XBRND to test it.

DISK1.zip 9.08 kB · 2 downloads

 

Next work... Calculate PI with it!

 

Guillaume.

 

XB RND sucks for being so slow that it is insane just to have that precision! 

I mean why do you need 9 decimal places of Random Numbers for plotting  on a 24x32 screen or even a 49152 pixel screen as 1 decimal place is more then enough.

What TI should have done was have a high precision version RAND and a non precision one RND.

This is why so many games on the TI99/4A sucked as being so slow.

 

By the way:

FOR I=1 TO 1000::A=RND::NEXT I

I did this in RXB and it took 22 seconds in GPL!

 

  • Like 1
Link to comment
Share on other sites

19 minutes ago, RXB said:

By the way:


FOR I=1 TO 1000::A=RND::NEXT I

I did this in RXB and it took 22 seconds in GPL!

 

Damned!! That's a good chrono !

 

So, you'll be able to test one of my programs with RXB! (the first in the list)

 

I have put three programs in this ZIP, all allow to save the current state. The second one even saves its current random seed.

 

DSK1.PIXBALONE

That uses only Extended Basic, full precision RND, 34 seconds for 200 points.

Speed = 5,88 point / sec

 

DSK1.PIXBLIKE

Thats uses MLC, but simulates the full precision RND, 34 seconds for 1000 points.

Speed = 29,41 point / sec

 

DSK1.PIRANDOM

Thats uses MLC with its fast RND routine from 0 to 4095, 15 seconds for 25000 points.

Speed = 1666,67 point / sec

DISK1.zip

 

Guillaume.

  • Like 2
Link to comment
Share on other sites

2 hours ago, moulinaie said:

Damned!! That's a good chrono !

 

So, you'll be able to test one of my programs with RXB! (the first in the list)

 

I have put three programs in this ZIP, all allow to save the current state. The second one even saves its current random seed.

 


DSK1.PIXBALONE

That uses only Extended Basic, full precision RND, 34 seconds for 200 points.

Speed = 5,88 point / sec

 


DSK1.PIXBLIKE

Thats uses MLC, but simulates the full precision RND, 34 seconds for 1000 points.

Speed = 29,41 point / sec

 


DSK1.PIRANDOM

Thats uses MLC with its fast RND routine from 0 to 4095, 15 seconds for 25000 points.

Speed = 1666,67 point / sec

DISK1.zip 11.71 kB · 1 download

 

Guillaume.

RXB uses the same RND as TI Basic now, so what ever speed TI Basic does it will be with a second  of that.

Just using my wrist watch on both RXB and TI Basic for 

FOR I=1 TO 1000::A=RND::NEXT I

both resulted in about 20 or up to 22 seconds, I would say both average about 21 seconds.

 

So if you test you program on TI Basic should be almost the same for RXB, unless you use a ton of slower TI Basic routines.

Link to comment
Share on other sites

  • 2 years later...

With all the discussion about benchmarking going on, I remembered that thread about the PI Calculator program I wrote in XB back in 2007. 

I wonder if I could engage the Forth gurus to take a shot at it for timing comparison.

 

Spoiler

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

 

 

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

 

Picalc.dsk

  • Like 1
Link to comment
Share on other sites

Some sort of fixed-point Forth implementation would be great. Build a fixed-point library that can do simple FP arithmetic, and then code an appropriate PI calculation algorithm using the FP library.

 

I'm not going to have the time to look at that though, unfortunately. Too much RL stuff in the way!

Edited by Willsy
typos
  • Like 1
Link to comment
Share on other sites

6 hours ago, Vorticon said:

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

 

Abasic 408 - MAME - GenMod Geneve (2:17)
 
The tool chain for doing this was super easy!
 
I cut the code from Atariage
Saved the code in notepad to a file
Used TI image tool to import an extended basic file to a hard drive image.
 
Ran the program:

image.png.1e02b4768fd726b2d59b5e952d135dc3.png

 

The run took 2:17

image.png.898dd26d1fa1efe018e2abc0255bdf26.png

 

image.png.7bb8a8f71312c57507ecf6b47d924ac4.png

 

 

 

 

 

 

  • Like 3
Link to comment
Share on other sites

22 hours ago, Willsy said:

Some sort of fixed-point Forth implementation would be great. Build a fixed-point library that can do simple FP arithmetic, and then code an appropriate PI calculation algorithm using the FP library.

I'm not going to have the time to look at that though, unfortunately. Too much RL stuff in the way!

I might take a crack at it. I'd like to use Lee's FbForth though as I have not had a chance to play with it yet. 

Just need some time to re-familiarize myself with Forth in general as it has been a while. Amazing how memory fails with disuse! 

  • Like 4
Link to comment
Share on other sites

2 minutes ago, Vorticon said:

I might take a crack at it. I'd like to use Lee's FbForth though as I have not had a chance to play with it yet. 

Just need some time to re-familiarize myself with Forth in general as it has been a while. Amazing how memory fails with disuse! 

That works for me. If you write a fixed-point library for fbForth I'll port it over to TF. I doubt it will need much alteration. :thumbsup:

  • Like 4
Link to comment
Share on other sites

As it happens, when you posted, I had this page open so I saw it immediately, and I also had Forth Dimensions magazine, Volume IV, issue I open, which has a ruck of fixed point articles! Whet your appetite (and reacquaint yourself with Forth here!). 

  • Like 2
  • Thanks 2
Link to comment
Share on other sites

4 minutes ago, Willsy said:

As it happens, when you posted, I had this page open so I saw it immediately, and I also had Forth Dimensions magazine, Volume IV, issue I open, which has a ruck of fixed point articles! Whet your appetite (and reacquaint yourself with Forth here!). 

Ah very cool bedtime reading! ? This actually looks very interesting. Thanks!

  • Like 3
Link to comment
Share on other sites

On 5/4/2022 at 8:56 AM, Vorticon said:

With all the discussion about benchmarking going on, I remembered that thread about the PI Calculator program I wrote in XB back in 2007. 

I wonder if I could engage the Forth gurus to take a shot at it for timing comparison.

 

  Reveal hidden contents


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

 

 

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

 

Picalc.dsk 180 kB · 2 downloads

 

I can't take credit for the code. It was written by Ed of DxForth fame in Australia. Ed makes DxForth for other legacy op. systems like CPM and DOS.

According to the credits some of it goes back 1991.

It uses Machin's algorithm. Most of the programs I have seen use the "spigot".

I don't think they are the same so this could be an invalid comparison.

It's pretty hairy code. It does some pretty neat dynamic memory allocation as well. 

 

The video is running on my indirect threaded system.  It would go a bit faster on direct threading. I will try it.

 

Spoiler

\ PI.FTH from DxForth
\
\ Revised 2015-02-09  es
\
\ Compute Pi to an arbitrary precision. Uses Machin's
\ formula:  pi/4 = 4 arctan(1/5) - arctan(1/239)
\
\ Compile with 16-bit DX-Forth: FORTH - INCLUDE PI.F BYE
\ Compile with CAMEL99 Forth: INCLUDE DSK*.PI  ( where * is your drive no.)
\
\ This 16-bit implementation allows up to 45,808 digits
\ to be computed before arithmetic overflow occurs.
\
\ The code can be used on 32-bit targets with appropriate
\ changes:
\
\   16-bit             32-bit
\
\   10000 Multiply     100000000 Multiply
\   <# # # # # #>      <# # # # # # # # # #>
\   4 +loop            8 +loop
\   525 um/mod         1050 um/mod
\                      remove 'digits > 45808' warning
\
\ Acknowledgements:
\
\   Roy Williams, Feb 1994
\   J. W. Stumpel, May 1991
\   E. Ford, Aug 2009
\   R. Bishop, Aug 1978
\
\ This code is PUBLIC DOMAIN. Use at your own risk.

\ Modified for Camel99 Forth  Mar 2021 Fox
NEEDS D.	  FROM DSK1.DOUBLE
NEEDS DUMP  FROM DSK1.TOOLS
NEEDS VALUE FROM DSK1.VALUES
NEEDS .R    FROM DSK1.UDOTR
NEEDS ELAPSE FROM DSK1.ELAPSE

DECIMAL
0 VALUE POWER  ( adr)
0 VALUE TERM   ( adr)
0 VALUE RESULT ( adr)
0 VALUE SIZE   ( n)

VARIABLE CARRY

: ADD ( -- )
  0 CARRY !
  RESULT
  0 SIZE 1- DO
    I CELLS OVER + ( res) DUP @ 0
    I CELLS TERM + @ 0  D+  CARRY @ M+
    ( hi) CARRY !  ( lo) SWAP ( res) !
  -1 +LOOP  DROP ;

: SUBTRACT ( -- )
  0 CARRY !
  RESULT
  0 SIZE 1- DO
    I CELLS OVER + ( RES) DUP @ 0
    I CELLS TERM + @ 0  D-  CARRY @ M+
    ( HI) CARRY !  ( LO) SWAP ( RES) !
  -1 +LOOP  DROP ;

0 VALUE FACTOR

\ scan forward for cell containing non-zero
: +INDEX ( ADR -- ADR INDEX )
    -1
    BEGIN 1+ DUP SIZE -
    WHILE
       2DUP CELLS + @
    UNTIL
    THEN ;

: (DIVIDE)
  ?DO
     I CELLS OVER + ( res)
     DUP @  CARRY @  FACTOR  UM/MOD
    ( quot) ROT ( res) !  ( rem) CARRY !
  LOOP ;

: DIVIDE ( ADR FACTOR -- )
  TO FACTOR   0 CARRY !  +INDEX
  ( adr index )  SIZE SWAP
  (DIVIDE)
  DROP ;

\ scan backward for cell containing non-zero
: -INDEX ( adr -- adr index )
    SIZE
    BEGIN 1- DUP
    WHILE
       2DUP CELLS + @
    UNTIL
    THEN ;

: MULTIPLY ( adr factor -- )
  TO FACTOR   0 CARRY !  -INDEX
  ( adr index )  0 SWAP
  DO
    I CELLS OVER + ( res)
    DUP @  FACTOR  UM*  CARRY @ M+
    ( hi) CARRY !  ( lo) SWAP ( res) !
  -1 +LOOP
  DROP ;

: COPY ( -- ) POWER TERM SIZE CELLS CMOVE ; \

: ZERO? ( result -- f )  +INDEX NIP SIZE = ;

0 VALUE PASS
VARIABLE EXP
VARIABLE SIGN

: DIVISOR ( -- N )
  PASS 1 = IF  5  ELSE  239  THEN ;

: ERASE  0 FILL ;

: INITIALIZE ( -- )
  POWER SIZE CELLS ERASE
  TERM  SIZE CELLS ERASE
  PASS 1 = IF  RESULT SIZE CELLS ERASE  THEN
  16  PASS DUP * / POWER !
  POWER  DIVISOR  DIVIDE
  1 EXP !  PASS 1- SIGN ! ;

0 VALUE NDIGIT

: CalcPi ( -- )
  NDIGIT 45800 U> IF
    ." Warning: digits > 45808 will be in error " CR
  THEN

  2 1+ 1
  DO
    I TO PASS
    INITIALIZE
    BEGIN
      COPY
      TERM  EXP @ DIVIDE
      SIGN @  DUP IF  SUBTRACT  ELSE  ADD  THEN
      0= SIGN !  2 EXP +!
      POWER  DIVISOR DUP *  DIVIDE
      POWER ZERO?
    UNTIL
  LOOP ;

\ Camel99 has OUT but I don't use in the Video driver
\ : CR  CR  OUT OFF ;
\ : #   #   OUT 1+! ;

DECIMAL
: (PRINT)
   ?DO
    0 OVER !
    DUP 10000 MULTIPLY
    DUP @  0 <# # # # # #> TYPE SPACE
    VCOL @ 3 + C/L @ > IF CR THEN
  4  +LOOP ;

: PRINT ( -- )
  CR
  RESULT  DUP @ 0 .R  [CHAR] . EMIT SPACE
  NDIGIT 0 (PRINT)
  DROP  CR ;

\ : GetNumber ( -- n )
\  CR ." How many digits do you want? "
\  PAD DUP 20 ACCEPT NUMBER? ABORT" Invalid" CR ;

: PI ( n -- )
( GetNumber ) DUP TO NDIGIT

  \ array size = ceil(ndigit / log10(2^16))
  109 UM* 525 UM/MOD SWAP ( rem) IF  1+  THEN
  2+  TO SIZE    ( extra for accurate last digits)

  50 ALLOT  ( expand the HOLD buffer space)

  HERE TO POWER   SIZE CELLS ALLOT
  HERE TO TERM    SIZE CELLS ALLOT
  HERE TO RESULT  SIZE CELLS ALLOT

  TICKER OFF
  CalcPi
  .ELAPSED
  PRINT
;
\ end

 

 

 

  • Like 2
Link to comment
Share on other sites

DTC is predictably a bit faster.

 

PI-100-DIGITS-DTC.png

 

Edit:  I am still perfecting this DTC system. I realized that I used a "colon" alias for CELLS and old code for VALUES.

Since I want this system to be optimized for speed not size, that is a faux pas.

So this the result with the changes.

image.png.f07af6c0c8e9afa3fb07fb0071209639.png

  • 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...