Π day
03-20-2022, 07:28 PM
Post: #41
 Frido Bohn Member Posts: 59 Joined: Jan 2015
RE: Π day
(03-20-2022 07:39 AM)Ángel Martin Wrote:  Sorry about that, here's the missing code in the [INTZL] routine:

Hi Ángel,

Thanks for providing - works like a charm now!

KR
Frido
03-21-2022, 07:24 AM
Post: #42
 Thomas Klemm Senior Member Posts: 2,105 Joined: Dec 2013
RE: π day
Just stumbled upon this visual proof of Gregory’s Theorem:
Quote:Let $$I_k$$ and $$C_k$$ denote the areas of regular $$k$$-gons inscribed in and circumscribed around a given circle.
Then for all $$n$$, $$I_{2n}$$ is the geometric mean of $$I_n$$ and $$C_n$$, and $$C_2n$$ is the harmonic mean of $$I_2n$$ and $$C_n$$; that is,

\begin{align} I_{2n}=\sqrt{I_{n}C_{n}} \end{align}

and

\begin{align} C_{2n}=\frac{2}{\frac{1}{I_{2n}}+\frac{1}{C_{n}}}. \end{align}

So I thought I give it a try using the || operator of the WP-34S:
Code:
001 RCL× Y 002 √ 003 || 004 STO+ x 005 RCL L 006 END

Initialize the stack:

4
ENTER
2
RTN

Calculate the sequence:

R/S
R/S
R/S
R/S

We end up with the same sequence as when using Vieta's formula.
Hint: Select YDON in MODE to display both $$I_k$$ and $$C_k$$.

3.3137085
2.82842712475

3.18259788
3.06146745892

3.15172491
3.12144515226

3.14411839
3.13654849055

3.14222363
3.14033115695

3.14175037
3.14127725093

3.14163208
3.14151380114

3.14160251
3.14157294037

3.14159512
3.14158772528

3.14159327
3.14159142151

3.14159281
3.14159234557

3.14159269
3.14159257658

3.14159266
3.14159263434

3.14159266
3.14159264878

3.14159265
3.14159265239

3.14159265
3.14159265329

3.14159265
3.14159265351

3.14159265
3.14159265357

3.14159265
3.14159265359

3.14159265
3.14159265359
03-21-2022, 04:03 PM
Post: #43
 Frido Bohn Member Posts: 59 Joined: Jan 2015
RE: Π day
(03-21-2022 07:24 AM)Thomas Klemm Wrote:  So I thought I give it a try using the || operator of the WP-34S:

Hi Thomas,

How would the || operator translate into e.g., the HP41?

KR
Frido
03-21-2022, 05:27 PM
Post: #44
 Thomas Klemm Senior Member Posts: 2,105 Joined: Dec 2013
RE: Π day
(03-21-2022 04:03 PM)Frido Bohn Wrote:  How would the || operator translate into e.g., the HP41?

It calculates the resistor if X and Y are in parallel:

\begin{align} \frac{1}{R} = \frac{1}{X} + \frac{1}{Y} \end{align}

Something like:
Code:
1/X X<>Y 1/X + 1/X

But if one of X or Y is 0 you might want to rather use:

\begin{align} R = \frac{X \cdot Y}{X + Y} \end{align}

Thus maybe:
Code:
RCL× Y X<>Y RCL+ L /

However, you can't simply use this code snippet instead of || as register L won't contain the correct value.
But I have no idea how the WP-34S calculates the result.
03-21-2022, 05:54 PM
Post: #45
 Thomas Klemm Senior Member Posts: 2,105 Joined: Dec 2013
RE: π day
I came up with this for the HP-42S:
Code:
00 { 16-Byte Prgm } 01 RCL× ST Y 02 SQRT 03 X<>Y 04 RCL× ST Y 05 LASTX 06 RCL+ ST Z 07 ÷ 08 STO+ ST X 09 X<>Y 10 END
03-21-2022, 06:33 PM
Post: #46
 Thomas Klemm Senior Member Posts: 2,105 Joined: Dec 2013
RE: π day
(03-21-2022 05:27 PM)Thomas Klemm Wrote:  But I have no idea how the WP-34S calculates the result.

From xrom/parallel.wp34s:
Code:
/**************************************************************************/ /* The real parallel operator.  *  *    par(x, y) = (x . y) / (x + y) = 1 / ( 1/x + 1/y )  */     XLBL"PARL"         xIN DYADIC         x=0?             xOUT xOUT_NORMAL         ENTER[^]         RCL+ Z         x=0?             JMP ret_neginf         /         [times]         xOUT xOUT_NORMAL

Thus a bit shorter:
Code:
00 { 6-Byte Prgm } 01 ENTER 02 RCL+ ST Z 03 ÷ 04 × 05 END
03-21-2022, 10:45 PM
Post: #47
 Albert Chan Senior Member Posts: 2,680 Joined: Jul 2018
RE: Π day
(03-21-2022 07:24 AM)Thomas Klemm Wrote:  We end up with the same sequence as when using Vieta's formula.

Yes. For 2^n sided polygon, we have 2^(n+1) right triangles, with acute angle x = pi/2^n

area(inscribled triangle) ≤ area(sector) ≤ area(circumscribed triangle)
sin(x)*cos(x)/2 ≤ 1^2*x/2 ≤ 1*tan(x)/2
sin(2x)/(2x) ≤ 1 ≤ tan(x)/x

I(2^n)
= sin(2x)/(2x) * pi
= sin(pi/2^(n-1)) / (pi/2^(n-1)) * pi
= 2 / (cos(pi/4) * cos(pi/8) * ... * cos(pi/2^(n-1)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // Vieta's formula for pi

Of course, the formula does not required 2^n sided polygons
Example, we can start with hexagon: I(6) = 3/2*√3 ≈ 2.598, C(6) = 2*√3 ≈ 3.464

sin(2x)/(2x) = 1 - 2/3*x^2 + 2/15*x^4 + ...
tan(x)﻿ /﻿ x﻿ ﻿ ﻿ ﻿ ﻿ = 1 + 1/3*x^2 + 2/15*x^4 + ...

Circumscribed polygon area about twice as accurate as inscribed. (opposite sign)
We can correct for this.

Code:
10 N=6 @ A=SIN(PI/N*2)*N/2 @ B=TAN(PI/N)*N 20 A=SQRT(A*B) @ B=B*A/((B+A)*.5) @ N=N+N 30 DISP N,A,B,(A+2*B)/3 40 IF N<100 THEN 20
Code:
>RUN  12              3                    3.21539030918        3.14359353945  24              3.10582854123        3.15965994211        3.14171614182  48              3.13262861329        3.14608621514        3.14160034786  96              3.13935020305        3.14271459965        3.14159313412  192             3.14103195089        3.14187304998        3.14159268362

We can apply Richardson extrapolation, to remove O(x^4) error

>RES
﻿ 3.14159268362
>RES + (RES-3.14159313412)/(2^4-1)
﻿ 3.14159265359

For unit circle, circumscribed polygon perimeter = its area.
Inscribed n-sided polygon perimeter = Inscribed 2n-sided polygon area.
This make complicated polygon perimeter code un-necessary.
03-24-2022, 01:36 AM
Post: #48
 Gerson W. Barbosa Senior Member Posts: 1,603 Joined: Dec 2013
RE: Π day
(03-21-2022 10:45 PM)Albert Chan Wrote:
Code:
>RUN  12              3                    3.21539030918        3.14359353945  24              3.10582854123        3.15965994211        3.14171614182  48              3.13262861329        3.14608621514        3.14160034786  96              3.13935020305        3.14271459965        3.14159313412

Code:
                 a                     b             48   3.13935020304686722   3.14608621513143498                 c                     d   96   3.14103195089050965   3.14271459964536831    pi ~ (640*c*(28*c-139*a)/(c-4*a)+4*b*(3737+14256*b/(d-4*b))+1056*d-2568*a-6453*(a^2*b)^(1/3))/11655  pi ~ 3.141592653589793(44)
03-26-2022, 03:59 PM
Post: #49
 Albert Chan Senior Member Posts: 2,680 Joined: Jul 2018
RE: Π day
(03-24-2022 01:36 AM)Gerson W. Barbosa Wrote:
Code:
                 a                     b             48   3.13935020304686722   3.14608621513143498                 c                     d   96   3.14103195089050965   3.14271459964536831    pi ~ (640*c*(28*c-139*a)/(c-4*a)+4*b*(3737+14256*b/(d-4*b))+1056*d-2568*a-6453*(a^2*b)^(1/3))/11655  pi ~ 3.141592653589793(44)

How is formula derived ? With angle x = pi/n, it has error of O(x^10)
All dimensionless variables relative to size of pi (closer to 1 the better)

XCAS> a,b := sin(x)/x, tan(x)/x;
XCAS> c,d := a(x=x/2), b(x=x/2);
XCAS> f := (640*c*(28*c-139*a)/(c-4*a) + 4*b*(3737+14256*b/(d-4*b))
+1056*d - 2568*a - 6453*(a^2*b)^(1/3)) / 11655
:;
XCAS> float(series(f,x,0,13))

1 + 4.42801658079e-05*x^10 + 3.58291832057e-05*x^12 + x^14*order_size(x)

---

Another way to get pi, via asin().

r = (b-a)/(b+b) = (1-a/b)/2 = (1-cos(x))/2 = sin(x/2)^2
x/2 = asin(sqrt(r))
pi = 2n * asin(sqrt(r))

The problem with above is tiny errors of asin() get magnified by 2n.
A better way is to start with good pi estimate, then correct to pi.

A = sin(x)/x * pi
pi = A * (2*x/2) / (2*sin(x/2)*cos(x/2))
= A/sqrt(1-r) * asin(sqrt(r)) / sqrt(r)
= A/sqrt(1-r) * (1 + r/6 * (1 + r*3^2/(4*5) * (1 + r*5^2/(6*7) * (1 + r*7^2/(8*9) + ...

XCAS> r := (b-a)/(b+b) :;
XCAS> float(series(a/sqrt(1-r)*(1+r/6*(1+r*9/20*(1+r*25/42*(1+r*49/72)))),x,0,13))

1.0 - 2.18478116122e-05*x^10 + 5.77706557054e-06*x^12 + x^14*order_size(x)

With the terms shown, we already get relative error of pi to O(x^10)
And, it is trivial to extend for more precisions.

Tested with 48-sided polygon perimeter, inscribed and circumscribed.

lua> n, A, B = 48, 3.1393502030468672, 3.1460862151314350
lua> r = (B-A)/(B+B)
lua> A/sqrt(1-r) * (1+r/6*(1+r*9/20*(1+r*25/42*(1+r*49/72))))
3.141592653589793
03-26-2022, 05:37 PM
Post: #50
 Gerson W. Barbosa Senior Member Posts: 1,603 Joined: Dec 2013
RE: Π day
(03-26-2022 03:59 PM)Albert Chan Wrote:
(03-24-2022 01:36 AM)Gerson W. Barbosa Wrote:
Code:
                 a                     b             48   3.13935020304686722   3.14608621513143498                 c                     d   96   3.14103195089050965   3.14271459964536831    pi ~ (640*c*(28*c-139*a)/(c-4*a)+4*b*(3737+14256*b/(d-4*b))+1056*d-2568*a-6453*(a^2*b)^(1/3))/11655  pi ~ 3.141592653589793(44)

How is formula derived ?

The hard way:

Glad to know there are easier and better ways. Thank you very much!
03-26-2022, 11:24 PM
Post: #51
 Albert Chan Senior Member Posts: 2,680 Joined: Jul 2018
RE: Π day
(03-26-2022 03:59 PM)Albert Chan Wrote:  r = (b-a)/(b+b) = (1-a/b)/2 = (1-cos(x))/2 = sin(x/2)^2
x/2 = asin(sqrt(r))
...
A = sin(x)/x * pi
pi = A * (2*x/2) / (2*sin(x/2)*cos(x/2))
= A/sqrt(1-r) * asin(sqrt(r)) / sqrt(r)
= A/sqrt(1-r) * (1 + r/6 * (1 + r*3^2/(4*5) * (1 + r*5^2/(6*7) * (1 + r*7^2/(8*9) + ...

Above is for n-sided polygon perimeter
For n-sided polygon area, except for defintion of r, formula is exactly the same.

A = sin(2x)/(2x) * pi = sin(x)*cos(x)/x * pi
B = tan(x)/x * pi       = sin(x)/cos(x)/x * pi

r = (B-A)/B = 1-A/B = 1-cos(x)^2 = sin(x)^2

• n-sided inscribed polygon perimeter = 2n-side inscribed polygon area
• (r for n-sided polygon perimeter)     = (r for 2n-sided polygon area)

Proof for n-sided polygon perimeters implied proof of 2n-sided polygon area. QED

Note: we do not have to worry about odd-sided polygon area.
Geometrically, n required to be positive integer. But, algebraically, it does not.

From above defined A, B, this showed x = pi/n can be anything.
In other words, formula for doubling of polygon sides does not require integer sides.

A2 = sqrt(A*B) = sin(x)/x * pi

B2 = 2/(1/A2 + 1/B)
= 2*A2 / (1 + A2/B)
= 2*sin(x)/x * pi / (1+cos(x))
= sin(x)/(1+cos(x)) / (x/2) * pi
= tan(x/2) / (x/2) * pi
03-27-2022, 01:44 PM (This post was last modified: 03-30-2022 01:20 AM by Albert Chan.)
Post: #52
 Albert Chan Senior Member Posts: 2,680 Joined: Jul 2018
RE: Π day
If arugment for asin is big, taylor series take a long time to converge.

Ratio of coefficient = (2k+1)^2 / ((2k+2)*(2k+3)), approach 1 when k is huge.
For asin(1), this implied dropped term is almost same size as the last kept term.

We can use the "half-angle" formula for asin()
(05-31-2021 07:30 PM)Albert Chan Wrote:  $$\displaystyle\arcsin(x) = 2\arcsin\left( {x \over \sqrt{2\sqrt{1-x^2}+2}} \right)$$

Code:
function asinq(x) -- = asin(sqrt(x))     if x > 4e-4 then return 2 * asinq(0.5*x/(sqrt(1-x)+1)) end     return sqrt(x) * (1+x/6*(1+x*9/20/(1-x*25/42))) end

lua> asinq(1/4), pi/6
0.5235987755982989      0.5235987755982988
lua> asinq(2/4), pi/4
0.7853981633974483      0.7853981633974483
lua> asinq(3/4), pi/3
1.0471975511965979      1.0471975511965976
lua> asinq(4/4), pi/2
1.5707963267948966      1.5707963267948966

Another way is to use Carlson Elliptic Integrals: RC(1-x,1) = asin(√x)/(√x)

Code:
function RC(x,y, verbose) -- RC(1-x,1) = asin(sqrt(x))/sqrt(x)     local k = y + 2*sqrt(x*y)     if k==x+2*y then return sqrt(3/k) end     if verbose then print(sqrt(3/k)) end     return RC((x+k)*0.25, (y+k)*0.25, verbose) end

This is a simplified version of RF, RC(x,y) = RF(x,y,y)

It has no problem getting asin(1)/1 = pi/2, 1 sqrt per iteration

lua> RC(1-1, 1, true)
1.7320508075688772
1.5764775210064272
1.5711174700143078
1.5708159330462599
1.5707975451380392
1.570796402832061
1.5707963315455151
1.5707963270917837
1.5707963268134517
1.5707963267960565
1.5707963267949692
1.5707963267949012
1.570796326794897
1.5707963267948968
1.5707963267948968

Another example, from area of inscribed and circumsribed square, for pi

lua> A, B = 2, 4
lua> r = (B-A)/B -- = 1/2
lua> c = RC(1-r, 1, true) -- = asin(sqrt(1/2)) / sqrt(1/2) = pi/4 * sqrt(2)
1.1147379454918027
1.1109478170877694
1.1107345982528842
1.1107215960382895
1.1107207883059862
1.110720737898786
1.1107207347495225
1.110720734552712
1.1107207345404118
1.1107207345396428
1.1107207345395949
1.1107207345395917
1.1107207345395915

lua> A * c/sqrt(1-r)
3.141592653589793

lua> A * asinq(r)/sqrt(r-r*r)
3.141592653589793
03-27-2022, 04:00 PM
Post: #53
 Albert Chan Senior Member Posts: 2,680 Joined: Jul 2018
RE: Π day
Here is getting C(x) = asin(√x)/(√x), using "half-angle" formula for asin

Code:
function C(x)   -- = asin(sqrt(x)) / sqrt(x)     if x < 4e-4 then return 1+x/6*(1+x*9/20/(1-x*25/42)) end     local y = 0.5/(sqrt(1-x)+1)     y = 2*sqrt(y) * C(x*y)     print('C( '..x..' ) = ',y)     return y end

Redo both previous examples in 1 go, with less calculations (2 sqrt per recursion)

Code:
lua> C(1) C( 0.0006022718974138037 ) =    1.0001004058641836 C( 0.002407636663901557 ) =     1.0004017081549652 C( 0.009607359798384776 ) =     1.0016081890839748 C( 0.038060233744356624 ) =     1.0064545427995637 C( 0.14644660940672624 ) =      1.0261721529770307 C( 0.5 ) =      1.1107207345395913 C( 1 ) =        1.5707963267948963
03-31-2022, 02:04 AM
Post: #54
 ttw Member Posts: 284 Joined: Jun 2014
RE: Π day
This (new) paper has some nice error estimates for various summation orders.

As a bonus, there's also a nice skeleton code of each method.
04-02-2022, 03:14 AM
Post: #55
 Ren Member Posts: 184 Joined: Mar 2016
RE: Π day
(03-19-2022 11:18 AM)Ángel Martin Wrote:  Heresy or new truth?

Here's an interesting reading for y'all:
https://www.theverge.com/tldr/2018/3/14/...iday-truth

I have eaten pie,
I'm not even sure what I would bring to share with cow-orkers on June 28th to celebrate Tau Day.

10B, 10BII, 12C, 14B, 15C, 16C, 17B, 18C, 19BII, 20b, 22, 29C, 35, 38G, 39G, 41CV, 48G, 97
04-02-2022, 11:12 AM (This post was last modified: 04-02-2022 11:14 AM by floppy.)
Post: #56
 floppy Senior Member Posts: 369 Joined: Apr 2021
RE: Π day
(03-14-2022 11:06 PM)ttw Wrote:  One can often apply various sequence transformations to speed up these series. A couple of new methods are given here.

https://arxiv.org/pdf/1702.07199.pdf

https://people.mpim-bonn.mpg.de/zagier/f...lltext.pdf

In most record-setting attempts, the desired result is a fraction with really big integers, even bigger than 355/113. An (to me, anyway) amusing feature is that most of the time is taken by the final long division. Third order, first order, fourth order, etc., methods all take about the same amount of time as the time for arithmetic grows pretty fast.
Interesting. When I see PI, I am remembering of page 51 of the "exact" PI formula