Threaded Mode | Linear Mode

Albert Chan · 06-23-2024, 06:00 PM

(06-23-2024 03:00 AM)Gerson W. Barbosa Wrote: \(\pi \approx \pi{'} \left[1+\frac{1}{6}\left(\frac{\pi{'}}{2^{n+1} }\right)^2\right]\left[1+\frac{3}{40}\left(\frac{\pi{'}}{2^{n+1} }\right)^4\right]\left[1+\frac{9}{280}\left(\frac{\pi{'}}{2^{n+1} }\right)^6\right] \)

Correction factor had error of O(z^8), where z=pi'/2^(n+1)
If we use Horner's scheme, we might as well use correct z^8 coefficient.
It made code shorter (less terms), also more accurate, error = O(z^10)

XCas> expand( (1+1/6*z^2) * (1+3/40*z^4) * (1+9/280*z^6) )

1 + 1/6*z^2 + 3/40*z^4 + 5/112*z^6 + 3/560*z^8 + 27/11200*z^10 + 9/22400*z^12

Replaced this with asin(z)/z, upto z^8 only

XCas> expand( 1+1/6*z^2 * (1+3²/(4*5)*z^2 * (1+5²/(6*7)*z^2 * (1+7²/(8*9)*z^2))) )

1 +1/6*z^2 + 3/40*z^4 + 5/112*z^6 + 35/1152*z^8

To make code stack-only, we can't store n, and had to recover 2^(n+1) factor.
Code was using 2^(n+1) = IP((355/113)/z), but this will fail if n is too big.

Instead, I use log to recover n from z (Yes, this is stupid, all we need is to store n)

z = pi' / 2^(n+1)
log2(z) = log2(pi') - (n+1)
-(n+1) ≈ log2(z) - 1.6515
-(n+1) = IP(log2(z)) - 2

n=1 --> IP(log2(z)) - 2 = IP(-0.50) - 2 = -2
n=9 --> IP(log2(z)) - 2 = IP(-8.35) - 2 = -10
n=99 --> IP(log2(z)) - 2 = IP(-98.35) - 2 = -100
n=999 --> IP(log2(z)) - 2 = IP(-998.35) - 2 = -1000
...

Code:

00 { 79-Byte Prgm }

01▸LBL "Vh"

02 2

03 0

04 SIGN

05▸LBL 00

06 2

07 STO× ST Z

08 RCL+ ST L

09 SQRT

10 ×

11 DSE ST Z

12 GTO 00

13 1/X          @ z = sin(pi/2^(n+1)) = pi'/2^(n+1)

14 X↑2

15 ENTER

16 ENTER

17 ENTER

18 49

19 ×

20 72

21 ÷

22 ×

23 +

24 25

25 ×

26 42

27 ÷

28 ×

29 +

30 .45

31 ×

32 ×

33 +

34 6

35 ÷            @ eps ≈ asin(z)/z - 1

36 X<>Y

37 SQRT

38 ENTER

39 LN           @ ln(z)   z   eps

40 2

41 LN

42 ÷            @ log2(z) z   eps   eps

43 IP

44 2

45 -

46 LASTX

47 X<>Y

48 Y↑X

49 ÷            @ pi'   eps   eps   eps

50 ×

51 LASTX

52 +

53 END

1 XEQ "Vh" --> 3.138316846974900201597524772067293
10 XEQ "Vh" --> 3.141592653589793238462643383274432
11 XEQ "Vh" --> 3.141592653589793238462643383279499
12 XEQ "Vh" --> 3.141592653589793238462643383279504

Gerson W. Barbosa · 06-24-2024, 09:42 PM

(06-23-2024 06:00 PM)Albert Chan Wrote:
(06-23-2024 03:00 AM)Gerson W. Barbosa Wrote: \(\pi \approx \pi{'} \left[1+\frac{1}{6}\left(\frac{\pi{'}}{2^{n+1} }\right)^2\right]\left[1+\frac{3}{40}\left(\frac{\pi{'}}{2^{n+1} }\right)^4\right]\left[1+\frac{9}{280}\left(\frac{\pi{'}}{2^{n+1} }\right)^6\right] \)

Correction factor had error of O(z^8), where z=pi'/2^(n+1)
If we use Horner's scheme, we might as well use correct z^8 coefficient.
It made code shorter (less terms), also more accurate, error = O(z^10)

XCas> expand( (1+1/6*z^2) * (1+3/40*z^4) * (1+9/280*z^6) )

1 + 1/6*z^2 + 3/40*z^4 + 5/112*z^6 + 3/560*z^8 + 27/11200*z^10 + 9/22400*z^12

Replaced this with asin(z)/z, upto z^8 only

XCas> expand( 1+1/6*z^2 * (1+3²/(4*5)*z^2 * (1+5²/(6*7)*z^2 * (1+7²/(8*9)*z^2))) )

1 +1/6*z^2 + 3/40*z^4 + 5/112*z^6 + 35/1152*z^8

I could have used up to four correction terms in the format I’ve been using, but I have used only three as I wanted to compare the results with those of another program using COMB, as that one is limited to only three factors. But my program would be even longer, so it’s great to know there’s a way to do it more compactly.

(06-23-2024 06:00 PM)Albert Chan Wrote: To make code stack-only, we can't store n, and had to recover 2^(n+1) factor.
Code was using 2^(n+1) = IP((355/113)/z), but this will fail if n is too big.

Instead, I use log to recover n from z (Yes, this is stupid, all we need is to store n)

z = pi' / 2^(n+1)
log2(z) = log2(pi') - (n+1)
-(n+1) ≈ log2(z) - 1.6515
-(n+1) = IP(log2(z)) - 2

n=1 --> IP(log2(z)) - 2 = IP(-0.50) - 2 = -2
n=9 --> IP(log2(z)) - 2 = IP(-8.35) - 2 = -10
n=99 --> IP(log2(z)) - 2 = IP(-98.35) - 2 = -100
n=999 --> IP(log2(z)) - 2 = IP(-998.35) - 2 = -1000
...

IP((355/113)/z) works for the small n we would need for maximum accuracy, 14 or 15 in my example with three factors. This starts to fail at n=23. Also, if more factors are used even less iterations are required (11 or 12 in your example).

(06-23-2024 06:00 PM)Albert Chan Wrote:

Code:

00 { 79-Byte Prgm } 01▸LBL "Vh" 02 2 03 0 04 SIGN 05▸LBL 00 06 2 07 STO× ST Z 08 RCL+ ST L 09 SQRT 10 × 11 DSE ST Z 12 GTO 00 13 1/X @ z = sin(pi/2^(n+1)) = pi'/2^(n+1) 14 X↑2 15 ENTER 16 ENTER 17 ENTER 18 49 19 × 20 72 21 ÷ 22 × 23 + 24 25 25 × 26 42 27 ÷ 28 × 29 + 30 .45 31 × 32 × 33 + 34 6 35 ÷ @ eps ≈ asin(z)/z - 1 36 X<>Y 37 SQRT 38 ENTER 39 LN @ ln(z) z eps 40 2 41 LN 42 ÷ @ log2(z) z eps eps 43 IP 44 2 45 - 46 LASTX 47 X<>Y 48 Y↑X 49 ÷ @ pi' eps eps eps 50 × 51 LASTX 52 + 53 END

I don’t understand why one would want to run this for n=10000+ on Free 42. You code would be even more compact if is limited to 20 or so. Anyway, if that is desired you can save at least one more step. They say engineers round \(\pi\) to 4 just in case. Being an engineer, that’s exactly what I have done at step 42 :-)

Code:

00 { 79-Byte Prgm }

01▸LBL "VH"

02 2

03 0

04 SIGN

05▸LBL 00

06 2

07 STO× ST Z

08 RCL+ ST L

09 SQRT

10 ×

11 DSE ST Z

12 GTO 00

13 1/X

14 X↑2

15 ENTER

16 ENTER

17 ENTER

18 49

19 ×

20 72

21 ÷

22 ×

23 +

24 25

25 ×

26 42

27 ÷

28 ×

29 +

30 .45

31 ×

32 ×

33 +

34 6

35 ÷

36 LN1+X

37 E↑X

38 X<>Y

39 SQRT

40 ×

41 2

42 4             @ pi ~ 4

43 LASTX

44 ÷

45 LN

46 2

47 LN

48 ÷

49 IP

50 Y↑X

51 ×

52 END

Albert Chan · 06-25-2024, 12:41 AM

(06-24-2024 09:42 PM)Gerson W. Barbosa Wrote: They say engineers round \(\pi\) to 4 just in case.
Being an engineer, that’s exactly what I have done at step 42 :-)

pi = 4 ... and it actually work! Big Grin

Turns out I don't need to recover 2^(n+1). It was already there!

This version is going for simplicity, which turns out pretty good.

ε = asin(z)/z - 1 ≈ z^2/6 / (1 - 9/20*z^2) = 1 / (6/z^2 - 2.7)

Code:

00 { 38-Byte Prgm }

01▸LBL "Vh"

02 2

03 0

04 SIGN

05▸LBL 00

06 2

07 STO× ST Z

08 RCL+ ST L

09 SQRT

10 ×

11 DSE ST Z

12 GTO 00       @ 1/z     2^(n+1)

13 ÷     

14 2.7   

15 6     

16 LASTX        @ 1/z     6     2.7   pi'

17 X↑2   

18 ×     

19 -            @ -1/eps  pi'   pi'   pi'

20 ÷

21 -

22 END

1 --> 3.132559073643629893044600829969869
17 --> 3.141592653589793238462643383279402
18 --> 3.141592653589793238462643383279501
19 --> 3.141592653589793238462643383279502

Gerson W. Barbosa · 06-25-2024, 01:19 PM

(06-25-2024 12:41 AM)Albert Chan Wrote: [quote='Gerson W. Barbosa' pid='188722' dateline='1719265376']

Turns out I don't need to recover 2^(n+1). It was already there!

Code:

… 05▸LBL 00 06 2 07 STO× ST Z 08 RCL+ ST L 09 SQRT 10 × 11 DSE ST Z 12 GTO 00 @ 1/z 2^(n+1) 13 …

Yes, except that it is difficult to use it when we need all the stack for something else afterwards. Unless, of course, we use a numbered register, as in Vh and Vf below. The latter uses a formula to generate the coefficients in the format (1 + n/d), but this works only for the three first coefficients.

Code:

00 { 92-Byte Prgm }

01▸LBL "Vh"

02 2

03 0

04 SIGN

05▸LBL 00

06 2

07 STO× ST Z

08 RCL+ ST L

09 SQRT

10 ×

11 DSE ST Z

12 GTO 00

13 1/X

14 ×

15 STO 01

16 LASTX

17 X↑2

18 ENTER

19 ENTER

20 ENTER

21 9

22 ×

23 22400

24 ÷

25 27

26 11200

27 ÷

28 +

29 ×

30 3

31 560

32 ÷

33 +

34 ×

35 5

36 112

37 ÷

38 +

39 ×

40 3

41 40

42 ÷

43 +

44 ×

45 6

46 1/X

47 +

48 ×

49 1

50 +

51 RCL× 01

52 END

Code:

00 { 85-Byte Prgm }

01▸LBL "Vf"

02 2

03 0

04 SIGN

05▸LBL 00

06 2

07 STO× ST Z

08 RCL+ ST L

09 SQRT

10 ×

11 DSE ST Z

12 GTO 00

13 1/X

14 ×

15 STO 01

16 LASTX

17 X↑2

18 3

19▸LBL 01

20 X<>Y

21 3

22 RCL ST Z

23 +/-

24 NOT

25 Y↑X

26 RCL ST Y

27 R↑

28 Y↑X

29 STO× ST Y

30 CLX

31 2

32 LASTX

33 Y↑X

34 STO÷ ST Y

35 X<> ST L

36 ENTER

37 SIGN

38 STO+ ST Y

39 X<> ST L

40 STO+ ST Y

41 COMB

42 STO÷ ST Y

43 R↓

44 LASTX

45 SIGN

46 STO+ ST Y

47 R↓

48 STO× 01

49 X<> ST L

50 DSE ST X

51 GTO 01

52 X<> 01

53 END

Albert Chan · 06-25-2024, 02:46 PM

(06-25-2024 12:41 AM)Albert Chan Wrote: This version is going for simplicity, which turns out pretty good.

ε = asin(z)/z - 1 ≈ z^2/6 / (1 - 9/20*z^2) = 1 / (6/z^2 - 2.7)

Here is for next pade approximation

ε = asin(z)/z - 1 ≈ z^2/6 * (1 + 9/20*z^2 / (1 - 25/42*z^2))

ε ≈ (1 + 2.7/(x-25/7)) / x , where x = 6/z^2

Code:

00 { 50-Byte Prgm }

01▸LBL "Vh"

02 2

03 0

04 SIGN

05▸LBL 00

06 2

07 STO× ST Z

08 RCL+ ST L

09 SQRT

10 ×

11 DSE ST Z

12 GTO 00       @ 1/z   2^(n+1)

13 ÷

14 LASTX

15 X↑2

16 6           

17 ×            @ x = 6/z^2

18 DUP

19 .28

20 1/X          @ 25/7    x    x    pi'

21 -

22 2.7

23 X<>Y

24 ÷

25 LN1+X

26 E↑X

27 ÷            @ 1/eps   pi'  pi'  pi'

28 ÷

29 +

30 END

13 --> 3.141592653589793238462643383279481
14 --> 3.141592653589793238462643383279502

Gerson W. Barbosa · (This post was last modified: 07-01-2024 06:00 AM by Gerson W. Barbosa.)

(06-23-2024 06:00 PM)Albert Chan Wrote: ...

Code:

00 { 79-Byte Prgm } 01▸LBL "Vh" 02 2 03 0 04 SIGN 05▸LBL 00 06 2 07 STO× ST Z 08 RCL+ ST L 09 SQRT 10 × 11 DSE ST Z 12 GTO 00 13 1/X @ z = sin(pi/2^(n+1)) = pi'/2^(n+1) 14 X↑2 15 ENTER 16 ENTER 17 ENTER 18 49 19 × 20 72 21 ÷ 22 × 23 + 24 25 25 × 26 42 27 ÷ 28 × 29 + 30 .45 31 × 32 × 33 + 34 6 35 ÷ @ eps ≈ asin(z)/z - 1 36 X<>Y 37 SQRT 38 ENTER 39 LN @ ln(z) z eps 40 2 41 LN 42 ÷ @ log2(z) z eps eps 43 IP 44 2 45 - 46 LASTX 47 X<>Y 48 Y↑X 49 ÷ @ pi' eps eps eps 50 × 51 LASTX 52 + 53 END

1 XEQ "Vh" --> 3.138316846974900201597524772067293
10 XEQ "Vh" --> 3.141592653589793238462643383274432
11 XEQ "Vh" --> 3.141592653589793238462643383279499
12 XEQ "Vh" --> 3.141592653589793238462643383279504

Same results:

1 4 XEQ "VCW" --> 3.138316846974900201597524772067292
10 4 XEQ "VCW" --> 3.141592653589793238462643383274435
11 4 XEQ "VCW" --> 3.141592653589793238462643383279500
12 4 XEQ "VCW" --> 3.141592653589793238462643383279504

Now with eight coefficients:

1 8 XEQ "VCW" --> 3.141497389177767990537566794077536
2 8 XEQ "VCW" --> 3.141592652530259569168830657515667
6 8 XEQ "VCW" --> 3.141592653589793238462643383279186
7 8 XEQ "VCW" --> 3.141592653589793238462643383279508
8 8 XEQ "VCW" --> 3.141592653589793238462643383279503

No optimization attempt at this stage. No coin flipping for program label – it looks better this way:-)
Formulae later, if I find a suitable LaTeX editor (tired of editing these by hand).

Code:

00 { 108-Byte Prgm }

01▸LBL "VCW"

02 STO 01

03 X<>Y

04 2

05 STO 02

06 0

07 SIGN

08 STO 03

09▸LBL 00

10 2

11 STO× ST Z

12 RCL+ ST L

13 SQRT

14 ×

15 DSE ST Z

16 GTO 00

17 1/X

18 ×

19 X<> 01

20 LASTX

21 X↑2

22 STO 04

23 1ᴇ-3

24 RCL× ST Z

25 1

26 +

27 X<>Y

28 0.5

29▸LBL 01

30 STO+ ST X

31 RCL ST Z

32 IP

33 STO÷ ST Y

34 STO+ ST X

35 +/-

36 NOT

37 ×

38 LASTX

39 SIGN

40 STO+ ST X

41 RCL+ ST L

42 RCL× 02

43 X<>Y

44 RCL÷ ST Y

45 X<>Y

46 X<> ST L

47 X<> ST Z

48 ×

49 STO+ 03

50 X<> ST L

51 RCL× 04

52 X<>Y

53 4

54 STO× 02

55 R↓

56 ISG ST Z

57 GTO 01

58 RCL 01

59 RCL× 03

60 END

Edit to change Subject – this looks indeed accelerated although further tests are still needed.

Gerson W. Barbosa · (This post was last modified: 07-03-2024 12:37 PM by Gerson W. Barbosa.)

DECIMAL BASIC PROGRAM

Code:

! Accelerated Viète by Chan-Wasicki Method

OPTION ARITHMETIC DECIMAL_HIGH

INPUT  PROMPT "number of decimal digits: ":nd

LET tm = TIME 

LET n = CEIL(4/3*SQR(nd))         ! n = number of factors in Viète's Product 

LET m = n                         ! m = number of coefficients in the correction polynomial

PRINT "n, m: ";n                

LET p = 2

LET s = 0

LET v = 1

FOR k = 1 TO n

   LET s = SQR(2 + s) 

   LET v = s*v

   LET p = p + p

NEXT k

LET pv = p/v

LET n = 1                          

LET w = 1/(v*v)

LET x = w

LET c = 1

LET p = 4

FOR k = 1 TO m

   LET n = 2*n*(2*k - 1)/k         ! n = coefficient numerator

   LET d = p*(2*k + 1)             ! d = coefficient denominator

   LET c = c + n*x/d

   LET x = x*w

   LET p = 4*p

NEXT k

LET r = TIME - tm

LET r$ = STR$(pv*c - INT(pv*c))

PRINT

PRINT STR$(INT(pv*c));

PRINT r$(0:1);

FOR i = 2 TO nd + 1

   PRINT r$(i:i);

   IF MOD((i - 1),10) = 0 THEN PRINT " ";

   IF MOD((i - 1),50) = 0 THEN 

      PRINT USING " (%%%%)": i - 1

      PRINT REPEAT$(" ",1 + LEN(STR$(INT(pv*c))));

   END IF

NEXT i

IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT

PRINT 

PRINT "Runtime: ";

PRINT  USING "0.##": r;

PRINT " seconds"

END

Code:

number of decimal digits: 1000

n, m:  43 

3.1415926535 8979323846 2643383279 5028841971 6939937510  (0050)

  5820974944 5923078164 0628620899 8628034825 3421170679  (0100)

  8214808651 3282306647 0938446095 5058223172 5359408128  (0150)

  4811174502 8410270193 8521105559 6446229489 5493038196  (0200)

  4428810975 6659334461 2847564823 3786783165 2712019091  (0250)

  4564856692 3460348610 4543266482 1339360726 0249141273  (0300)

  7245870066 0631558817 4881520920 9628292540 9171536436  (0350)

  7892590360 0113305305 4882046652 1384146951 9415116094  (0400)

  3305727036 5759591953 0921861173 8193261179 3105118548  (0450)

  0744623799 6274956735 1885752724 8912279381 8301194912  (0500)

  9833673362 4406566430 8602139494 6395224737 1907021798  (0550)

  6094370277 0539217176 2931767523 8467481846 7669405132  (0600)

  0005681271 4526356082 7785771342 7577896091 7363717872  (0650)

  1468440901 2249534301 4654958537 1050792279 6892589235  (0700)

  4201995611 2129021960 8640344181 5981362977 4771309960  (0750)

  5187072113 4999999837 2978049951 0597317328 1609631859  (0800)

  5024459455 3469083026 4252230825 3344685035 2619311881  (0850)

  7101000313 7838752886 5875332083 8142061717 7669147303  (0900)

  5982534904 2875546873 1159562863 8823537875 9375195778  (0950)

  1857780532 1712268066 1300192787 6611195909 2164201818  (1000)

Runtime: 0.26 seconds

The last three digits should be 989, but we are here at the limit of DECIMAL BASIC precision.
Like number of iterations for the Viète's Product and the correction polynomial, not sure whether this gives optimal efficiency. Anyway, sublinear convergence rate (the number of iterations depends on the square root of the number of digits).

Edited to fix a typo

----------

For one thousand digits, only 11 iterations of the first loop and 165 of the second loop would suffice. I have added six iterations to the latter as a safety margin for this range. Almost 2.5 times faster. Perhaps not too far from the ideal m/n ratio. The last two digits should be 89.

Code:

 ! Accelerated Viète by Chan-Wasicki Method

OPTION ARITHMETIC DECIMAL_HIGH

INPUT  PROMPT "number of decimal digits: ":nd

LET tm = TIME 

LET n = CEIL(9/26*SQR(nd))         ! n = number of factors in Viète's Product

LET m = 15*n + 6                   ! m = number of coefficients in the correction polynomial

PRINT "n, m: ";n;m                

LET p = 2

LET s = 0

LET v = 1

FOR k = 1 TO n

   LET s = SQR(2 + s) 

   LET v = s*v

   LET p = p + p

NEXT k

LET pv = p/v

LET n = 1                          

LET w = 1/(v*v)

LET x = w

LET c = 1

LET p = 4

FOR k = 1 TO m

   LET n = 2*n*(2*k - 1)/k         ! n = coefficient numerator

   LET d = p*(2*k + 1)             ! d = coefficient denominator

   LET c = c + n*x/d

   LET x = x*w

   LET p = 4*p

NEXT k

LET r = TIME - tm

LET r$ = STR$(pv*c - INT(pv*c))

PRINT

PRINT STR$(INT(pv*c));

PRINT r$(0:1);

FOR i = 2 TO nd + 1

   PRINT r$(i:i);

   IF MOD((i - 1),10) = 0 THEN PRINT " ";

   IF MOD((i - 1),50) = 0 THEN 

      PRINT USING " (%%%%)": i - 1

      PRINT REPEAT$(" ",1 + LEN(STR$(INT(pv*c))));

   END IF

NEXT i

IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT

PRINT 

PRINT "Runtime: ";

PRINT  USING "0.##": r;

PRINT " seconds"

END

Code:

number of decimal digits: 1000

n, m:  11  171 

3.1415926535 8979323846 2643383279 5028841971 6939937510  (0050)

  5820974944 5923078164 0628620899 8628034825 3421170679  (0100)

  8214808651 3282306647 0938446095 5058223172 5359408128  (0150)

  4811174502 8410270193 8521105559 6446229489 5493038196  (0200)

  4428810975 6659334461 2847564823 3786783165 2712019091  (0250)

  4564856692 3460348610 4543266482 1339360726 0249141273  (0300)

  7245870066 0631558817 4881520920 9628292540 9171536436  (0350)

  7892590360 0113305305 4882046652 1384146951 9415116094  (0400)

  3305727036 5759591953 0921861173 8193261179 3105118548  (0450)

  0744623799 6274956735 1885752724 8912279381 8301194912  (0500)

  9833673362 4406566430 8602139494 6395224737 1907021798  (0550)

  6094370277 0539217176 2931767523 8467481846 7669405132  (0600)

  0005681271 4526356082 7785771342 7577896091 7363717872  (0650)

  1468440901 2249534301 4654958537 1050792279 6892589235  (0700)

  4201995611 2129021960 8640344181 5981362977 4771309960  (0750)

  5187072113 4999999837 2978049951 0597317328 1609631859  (0800)

  5024459455 3469083026 4252230825 3344685035 2619311881  (0850)

  7101000313 7838752886 5875332083 8142061717 7669147303  (0900)

  5982534904 2875546873 1159562863 8823537875 9375195778  (0950)

  1857780532 1712268066 1300192787 6611195909 2164201968  (1000)

Runtime: 0.11 seconds

----------

Using a recurrent formula also for the denominator might be better...

Code:

! Accelerated Viète by Chan-Wasicki Method

OPTION ARITHMETIC DECIMAL_HIGH

INPUT  PROMPT "number of decimal digits: ":nd

LET n = CEIL(9/26*SQR(nd))         ! n = number of factors in Viète's Product

LET m = 15*n + 6                   ! m = number of coefficients in the correction polynomial

PRINT "n, m: ";n;m

LET m = 2*m - 1 

LET tm = TIME               

LET p = 2

LET s = 0

LET v = 1

FOR k = 1 TO n

   LET s = SQR(2 + s) 

   LET v = s*v

   LET p = p + p

NEXT k

LET pv = p/v

LET n = 1                          ! n = coefficient numerator

LET d = 1                          ! d = coefficient denominator

LET w = 1/(v*v)

LET x = w

LET c = 1

FOR k = 1 TO m STEP 2

   LET n = 4*n*k/(k + 1)               

   LET d = d*(4*k+8)/k                               

   LET c = c + n/d*x

   LET x = x*w

NEXT k

LET r = TIME - tm

LET r$ = STR$(pv*c - INT(pv*c))

PRINT

PRINT STR$(INT(pv*c));

PRINT r$(0:1);

FOR i = 2 TO nd + 1

   PRINT r$(i:i);

   IF MOD((i - 1),10) = 0 THEN PRINT " ";

   IF MOD((i - 1),50) = 0 THEN 

      PRINT USING " (%%%%)": i - 1

      PRINT REPEAT$(" ",1 + LEN(STR$(INT(pv*c))));

   END IF

NEXT i

IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT

PRINT 

PRINT "Runtime: ";

PRINT  USING "0.##": r;

PRINT " seconds"

END

...and faster!

Code:

  ...

  1857780532 1712268066 1300192787 6611195909 2164201968  (1000)

Runtime: 0.09 seconds

Albert Chan · 07-01-2024, 08:29 PM

I just realized I had done this before, in Arc SOHCAHTOA method

(03-31-2022 11:07 PM)Albert Chan Wrote: asinq(x) = 2 * asinq(x/2/(1+sqrt(1-x)))

lua> k, x = 4, .5 -- 4 * asinq(.5) = 4 * pi/4 = pi
lua> FMT = '%g * asinq(%.17g)\n'
lua> for i=1,5 do printf(FMT,k,x); k=k+k; x=0.5*x/(1+sqrt(1-x)) end
4 * asinq(0.5)
8 * asinq(0.14644660940672624)
16 * asinq(0.038060233744356624)
32 * asinq(0.009607359798384776)
64 * asinq(0.0024076366639015569)

Argument inside asinq, x = z^2 in Free42 code. Undo asinq, we have pi

lua> k*sqrt(x) -- = pi', with n=6
3.141277250932773
lua> _ * (1 + x/6*(1 + 9/20*x*(1 + 25/42*x*(1 + 49/72*x))))
3.141592653589793

Gerson W. Barbosa · (This post was last modified: 07-01-2024 09:06 PM by Gerson W. Barbosa.)

(07-01-2024 08:29 PM)Albert Chan Wrote: I just realized I had done this before, in Arc SOHCAHTOA method

(03-31-2022 11:07 PM)Albert Chan Wrote: asinq(x) = 2 * asinq(x/2/(1+sqrt(1-x)))

Pi day celebrations appear to give these kind of interesting byproducts. In this thread from 2009 started by Egan Ford there’s a link about Vieta Acceleration. It looks quite complicated and his program for 5000 digits requires a lot of memory. I would like to check how the above compares to it.

A good think about the correction polynomial it that its coefficients can be easily determined, what was something I was counting on in order to make it a valid method for computing \(\pi\). Wolfram Alpha returned a formula including COMB(2n, n) for the coefficients, which of course would not be practical for that purpose. I had already quit, when I evaluated the first six or seven such COMBs on the 42S and noticed this pattern:

1 2 3 4 5 6 7 8 9 10
3 6 10 15 21
10 20 35 56
35 70 126
126 252 462

Problem solved!