Π day

[Frido Bohn]

(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

[Thomas Klemm]

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

[Frido Bohn]

(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

[Thomas Klemm]

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

[Thomas Klemm]

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

[Thomas Klemm]

(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

[Albert Chan]

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

[Gerson W. Barbosa]

(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)

[Albert Chan]

(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

[Gerson W. Barbosa]

(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:

https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=188443

Glad to know there are easier and better ways. Thank you very much!

[Albert Chan]

(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

[Albert Chan]

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)
see https://www.hpmuseum.org/forum/thread-17...#pid148498

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

[Albert Chan]

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

[ttw]

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.

[Ren]

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

[floppy]

(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
http://semjonadlaj.com/SP/Computer+Algeb..._37-56.pdf
we have there an highspeed-convergency