[VA] SRC #009  Pi Day 2021 Special

03162021, 10:29 PM
(This post was last modified: 03172021 03:37 AM by Gerson W. Barbosa.)
Post: #21




RE: [VA] SRC #009  Pi Day 2021 Special
(03152021 08:59 PM)JF Garnier Wrote: I don't know and don't think there is  any relation that can be used to get \(\pi\) from e. Et pourtant, il y en a – and yet there is at least one: http://oeis.org/wiki/A_remarkable_formula_of_Ramanujan Really remarkable, isn’t it? P.S.: Yet another one (I had forgotten about the Gaussian Integral) \[{\int_{\infty}^{ \infty}\rm{e}^{{{x}}^{2}}\rm{d}x}=\sqrt{\pi}\] 

03172021, 12:49 AM
Post: #22




RE: [VA] SRC #009  Pi Day 2021 Special
(03162021 08:28 PM)robve Wrote: This one. I quote: "Also, what is the fun of doing math and calc exercises if we don't implement numerical integration ourselves?" Regards. V. All My Articles & other Materials here: Valentin Albillo's HP Collection 

03172021, 02:31 AM
Post: #23




RE: [VA] SRC #009  Pi Day 2021 Special
(03162021 07:17 PM)Albert Chan Wrote:(03142021 07:00 PM)Valentin Albillo Wrote: [*]d. Conversely, the volume enclosed by the ndimensional sphere of radius R is given by: Yes, it is kind of weird. But this is connecting two seemingly unrelated formulae, which is nice. 1. Taylor series: $$ e^x = 1+\frac{x}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots $$ 2. The volume of an nball with radius R: $$ V_n = \frac{\pi^{\frac{n}{2}}}{\Gamma(\frac{n}{2}+1)}R^n $$ The latter simplifies to $$ \frac{\pi^k}{k!} $$ for k=2n and R=1 (the conditions stated in the challenge). Therefore, the answer is e as the π root of the sum: $$ \sum_{k=0}^\infty \frac{\pi^k}{k!} = \mathrm{e}^\pi $$ I recognized the Taylor series after simplifying the sum's terms. Whenever you see a factorial in a denominator in a term of a series sum, there may be a Taylor series lurking beneath. Nice piece of natural pie...  Rob "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03172021, 02:32 AM
(This post was last modified: 03172021 11:26 AM by Albert Chan.)
Post: #24




RE: [VA] SRC #009  Pi Day 2021 Special
(03162021 09:09 PM)robve Wrote: $$ \arctan u  \arctan v = \arctan\frac{uv}{1+uv} $$ This is not quite right, LHS has possible range of ±pi, RHS is limited to ±pi/2 Correct identity is more complicated: see Sum of ArcTangents We could use this instead: atan(u) ± atan(v) = atan2(u±v , 1∓uv) (*) y = atan(x)  atan((x1)/(x+1)) // y undefined when x = 1 = atan2(x  (x1)/(x+1), 1 + x*(x1)/(x+1)) = atan2((x²+1)/(x+1), (x²+1)/(x+1)) // numerator always positive. If x > 1, 4*y = 4*atan(1) = pi If x < 1, 4*y = 4*(atan(1)  pi) = 3*pi (*) Proof is trivial: (1+u*i) * (1±v*i) = (1∓u*v) + (u±v)*i Phase angle of two sides matched, and we have above atan2 identity 

03172021, 03:03 AM
(This post was last modified: 03172021 02:14 PM by robve.)
Post: #25




RE: [VA] SRC #009  Pi Day 2021 Special
(03172021 02:32 AM)Albert Chan Wrote:(03162021 09:09 PM)robve Wrote: $$ \arctan u  \arctan v = \arctan\frac{uv}{1+uv} $$ Right. I did not include the two necessary conditions since these hold, note the (mod π): $$ \arctan u  \arctan v = \arctan\frac{uv}{1+uv}\quad \pmod \pi,\quad uv\neq 1 $$ EDIT: the arctan identity comes from $$ \tan \alpha \pm \tan \beta = \frac{\tan \alpha \pm \tan \beta}{1\mp \tan \alpha\, \tan \beta} $$ Since 0<=e<π and 0<=(e1)/(e+1)<π we have (e.g. verify numerically) $$\tan\arctan\mathrm{e}=\mathrm{e}; \quad \tan\arctan\frac{\mathrm{e}1}{\mathrm{e}+1} = \frac{\mathrm{e}1}{\mathrm{e}+1}$$ Then $$ \alpha=\arctan\mathrm{e};\quad \beta=\arctan \frac{\mathrm{e}1}{\mathrm{e}+1};\quad \frac{\tan \alpha  \tan \beta}{1+ \tan \alpha\, \tan \beta} = \frac{\mathrm{e}  \frac{\mathrm{e}1}{\mathrm{e}+1}}{1+\mathrm{e}\frac{\mathrm{e}1}{\mathrm{e}+1}} = 1 $$ Generally, tan has period π $$ \tan(k\pi+\theta) = \theta $$ for any integer k.  Rob "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03172021, 03:42 AM
Post: #26




RE: [VA] SRC #009  Pi Day 2021 Special
(03172021 12:49 AM)Valentin Albillo Wrote:(03162021 08:28 PM)robve Wrote: "Me/ourselves" and "implement integration" therefore "all programmers (should) write integration"? As in "Socrates is a man, Socrates is mortal, therefore all men are mortal"? Quotations matter. Not a programmer, a D in Ph'y after passing some BS, humbly not wanting to be a P in the A, just taking some late snacks on vintage stuff to honor those that came before us.  Rob "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03172021, 08:27 AM
Post: #27




RE: [VA] SRC #009  Pi Day 2021 Special
(03162021 10:29 PM)Gerson W. Barbosa Wrote:(03152021 08:59 PM)JF Garnier Wrote: I don't know and don't think there is  any relation that can be used to get \(\pi\) from e. OK, I see what you (and Valentin probably too) mean and I agree of course. By relation, I was (wrongly) limiting myself to finite expressions, like the arctan expression in Valentin's post. There are obviously many infinite sums and integrals involving e and producing pi in a way or another. JF 

03172021, 10:32 AM
Post: #28




RE: [VA] SRC #009  Pi Day 2021 Special
(03142021 07:00 PM)Valentin Albillo Wrote: [*]e. \(\pi\) also features in a song by Kate Bush (included in her 2005's album "Aerial") about a man who's utterly obsessed with the calculation of \(\pi\) (that could describe some of us here at the MoHPC). She sings more than a hundred digits of \(\pi\) and Love that song ;) https://www.youtube.com/watch?v=W8RE2NyAiJg 

03172021, 01:31 PM
Post: #29




RE: [VA] SRC #009  Pi Day 2021 Special
I love Kate Bush.
Greetings, Massimo +×÷ ↔ left is right and right is wrong 

03182021, 02:56 AM
Post: #30




RE: [VA] SRC #009  Pi Day 2021 Special
VA's posts are fascinating and responses are brilliant. His challenges encourages sleuthing using our advanced and vintage HP calculators and perhaps by writing some code to figure this all out.
To return the favor I hereby post two small and related challenges. These two "counter" challenges "invert" VA's common objective (if I may so): instead of writing code and (CAS) expressions, let's figure out what the given code does, find its formula and finally investigate who discovered it (online searching is permitted!). The first question of each of these two should be easy to answer. If you do not have a HP71B (I recently acquired mine ), then any Basiccapable machine can be used instead. This code is simple enough to easily convert to HPPL, RPN, Forth. a. Consider the following HP71B program: 10 P=SQR(2) 20 Q=P/2 30 DISP 2/Q 40 P=SQR(2+P) 50 Q=Q*P/2 60 IF P<2 GOTO 30 1. What constant does it compute? 2. What is the algebraic formula computed by this code for this constant? 3. Who discovered the formula? 4. Anything else that is interesting about this formula? b. Consider the following HP71B program: 10 P=1 20 FOR I=2 TO 1000 STEP 2 30 P=P*I/(I1) 40 DISP 2*P*P/(I+1) 50 NEXT I 1. What constant does it compute? 2. What is the algebraic formula computed by this code for this constant? 3. Who discovered the formula? 4. Anything else that is interesting about this formula? I will reply with the answers after VA posted his conclusions of the pi day challenge.  Rob "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03182021, 02:44 PM
Post: #31




RE: [VA] SRC #009  Pi Day 2021 Special
Hi, robve
Thanks for spending the time to add code challenges. Both codes are calculating pi by 2/(sin(x)/x) at x=pi/2, in 2 different ways. (a) pi is via Viète's formula: sin(x)/x = cos(x/2) cos(x/4) cos(x/8) ... At x = pi/2: cos(x/2) = cos(pi/4) = √(2) / 2 cos(x/4) = √((1 + cos(x/2))/2) = √(2 + √(2)) / 2 cos(x/8) = √((1 + cos(x/4))/2) = √(2 + √(2 + √(2))) / 2 ... Generated outputs, 2/Q = 2^k * sin(pi/2^k), k = 2, 3, 4, ... limit(2^k * sin(pi/2^k), k=∞) = 2^k * (pi/2^k) = pi // sin(ε) ≈ ε (b) pi is via Wallis products With roots of sin(x) = 0, ±pi, ±2pi, ±3pi, ..., and limit(sin(x)/x, x=0) = 1: sin(x)/x = (1(x/pi)²) * (1(x/(2pi))²) * (1(x/(3pi))²) ... At x=pi/2: LHS = 2/pi ≈ 0.63662 RHS = (11/4) * (11/16) * (11/36) ... = (1*3)/(2*2) * (3*5)/(4*4) * (5*7)/(6*6) ... BTW, it is more efficient (and accurate !) to calculate P=sin(x)/x, at x=pi/2: Bonus: code is shorter, and easier to understand. 10 P=1 20 FOR I=2 TO 1000 STEP 2 30 P=PP/(I*I) 40 DISP 2/P 50 NEXT I 

03202021, 09:36 PM
(This post was last modified: 03212021 01:33 AM by robve.)
Post: #32




RE: [VA] SRC #009  Pi Day 2021 Special
(03172021 02:32 AM)Albert Chan Wrote: We could use this instead: atan(u) ± atan(v) = atan2(u±v , 1∓uv) (*) A bit late to reply. Initially I also thought about matching angles between complex numbers in polar coordinates with atan2. Subtracting the angles gives the resulting angle of the complex quotient expressed in polar coordinates: \( \arctan \mathrm{e}  \arctan\frac{\mathrm{e}  1}{\mathrm{e} + 1} = \mathrm{atan2}(\mathrm{e},1)  \mathrm{atan2}(\mathrm{e}1,\mathrm{e}+1) \) then simplifying the quotient of the corresponding complex coordinates \( \frac{\sqrt{\mathrm{e}^2+1}}{\sqrt{(\mathrm{e}1)^2+(\mathrm{e}+1)^2}} \cdot \frac{1+\mathrm{e}\mathrm{i}}{\mathrm{e}+1+(\mathrm{e}1)\mathrm{i}} = \frac{1}{\sqrt{2}}\cdot \frac{1}{1\mathrm{i}} \) where the modulus of the denominator is \( \sqrt{2} \) and angle \( \mathrm{atan2}(1,1) \) i.e. representing \( \frac{1\angle 0}{\sqrt{2}\angle\arctan{}1} \). This gives \( \arctan \mathrm{e}  \arctan\frac{\mathrm{e}  1}{\mathrm{e} + 1} = 0  \arctan {}1 = \arctan 1 = \frac{\pi}{4} \). Just another way to prove this equation, which is an identity that holds for other values than e by the way (with constraints). (03172021 02:32 AM)Albert Chan Wrote: We could use this instead: atan(u) ± atan(v) = atan2(u±v , 1∓uv) Sure, but this is the same formula I had used in my previous post, albeit yours is in disguise using atan2 instead of atan, i.e. atan(u) ± atan(v) = atan((u±v)/(1∓uv)) = atan2(u±v , 1∓uv) the latter by definition and only if 1∓uv is nonzero, i.e. the necessary constraint I mentioned before.  Rob Minor edit: fix typo and \( \LaTeX \). Second edit: comment on atan2 versus arctan. "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03212021, 06:48 PM
(This post was last modified: 03212021 09:38 PM by robve.)
Post: #33




RE: [VA] SRC #009  Pi Day 2021 Special
(03182021 02:44 PM)Albert Chan Wrote: Thanks for spending the time to add code challenges. You're very welcome! Because this is going to be a long post, I will post the two parts a and b separately. a. 10 P=SQR(2) 20 Q=P/2 30 DISP 2/Q 40 P=SQR(2+P) 50 Q=Q*P/2 60 IF P<2 GOTO 30 As Albert replied correctly, the code computes Viète's formula, published in 1593: $$ \frac{2}{\pi} = \frac{\sqrt{2}}{2} \cdot \frac{\sqrt{2+\sqrt{2}}}{2} \cdot \frac{\sqrt{2+\sqrt{2+\sqrt{2}}}}{2} \cdots $$ In recurrence form: $$ \lim_{n\rightarrow\infty} \prod_{i=1}^n \frac{a_i}{2} = \frac{2}{\pi};\qquad a_1 = \sqrt{2},\quad a_n = \sqrt{2+a_{n1}} $$ In functional form: $$ \pi \approx \mathrm{viete}(n) = 2/\prod_{i=1}^n v(n)/2;\qquad v(n) = \begin{cases} \sqrt{2} & n=1 \\ \sqrt{2+v(n1)} & n>1 \end{cases} $$ Viète obtained his formula by comparing the areas of regular polygons with \( 2^n \) and \( 2^{n+1} \) sides inscribed in a circle. The first term in the product, \( \frac{\sqrt{2}}{2} \), is the ratio of areas of a square and an octagon, the second term is the ratio of areas of an octagon and a hexadecagon, etc.  Wikipedia See also "Mathematics: From the Birth of Numbers" by Jan Gullberg. This directly relates to Archimedes's famous work (ca.225BC) on approximating the area of a circle by polygons inside and outside the circle squeezing the circle: "At each stage, he needed to approximate sophisticated square roots, yet he never faltered. When he reached the 96gon, his estimate was \( \frac{6336}{2017\frac{1}{4}} > 3\frac{10}{71} \) ."  "Journey Through Genius" by William Dunham. Note that Euclid's Elements does not prove the ratio of the radius of the circle to its circumference \( 2\pi \). Sometimes confused with Euclid VI.33 that proves an important property of angles and arcs but does not compare them to circles with different radii "in equal circles [emphasis mine], angles, whether at the center or circumference, are in the same ratio to one another as the arcs on which they stand."  Oliver Byrne's Euclid 1847 VI.33 [image credit: Oliver Byrne's Euclid 1847] Viète's formula can also be written as $$ \frac{\pi}{2} = \frac{1}{\sqrt{\frac{1}{2}}} \cdot \frac{1}{\sqrt{\frac{1}{2}+\frac{1}{2}\sqrt{\frac{1}{2}}}} \cdot \frac{1}{\sqrt{\frac{1}{2}+\frac{1}{2}\sqrt{\frac{1}{2}+\frac{1}{2}\sqrt{\frac{1}{2}}}}} \cdots $$ In recurrence form: $$ \lim_{n\rightarrow\infty} \prod_{i=1}^n \frac{1}{a_i} = \frac{\pi}{2};\qquad a_1 = \sqrt{\frac{1}{2}},\quad a_n = \sqrt{\frac{1+a_{n1}}{2}} $$ In functional form: $$ \pi \approx \mathrm{viete}(n) = 2/\prod_{i=1}^n v(n);\qquad v(n) = \begin{cases} \sqrt{1/2} & n=1 \\ \sqrt{(1+v(n1))/2} & n>1 \end{cases} $$ While not based on Viète's formula, with our calculators we can quickly integrate the unit quarter circle defined by \( 1=x^2+y^2 \), i.e. integrate \( y = f(x) = \sqrt{1x^2} \) to get \( \pi \), using a square root: $$ \pi = 4\int_0^1 \sqrt{1x^2} \,dx $$ Proof (I've simplified this somewhat): $$ 4 \int_0^1 \sqrt{1x^2} \,dx = 4 \int_0^{\frac{\pi}{2}} \sqrt{1\sin^2 \theta} \cos \theta \,d\theta = 4 \int_0^{\frac{\pi}{2}} \sqrt{\cos^2 \theta} \cos \theta \,d\theta = 4 \int_0^{\frac{\pi}{2}} \cos^2 \theta \,d\theta = \pi $$ Programs I wrote today to illustrate Viete's formula My own functional programming language Husky: v(1) := (r,r) where r := sqrt(2)/2; v(n) := (t,q) where t := p*q where q := sqrt(2+2*r)/2 where (p,r) := v(n1). viete(n) := 2/p where (p,r) := v(n). > viete(20). and a listbased version: prod := foldr(*, 1). v(1) := [sqrt(2)/2]; v(n) := [sqrt(2+2*x)/2, x  xs] where x.xs := v(n1). viete(n) := 2/prod(v(n)). > viete(20). Haskell: v 1 = (r,r) where r = (sqrt 2)/2 v n = (t,q) where t = p*q q = (sqrt (2+2*r))/2 (p,r) = v (n1) viete n = 2/r where (r,_) = v n main = putStrLn (show (viete 20)) and a listbased version: prod = foldr (*) 1 viete n = 2/(prod (v n)) v 1 = [(sqrt 2)/2] v n = (sqrt (2+2*x))/2 : x : xs where x:xs = v (n1) main = putStrLn (show (viete 20)) Prolog: viete(N, P) : v(N, R, _), P is 2/R. v(1, R, R) : R is sqrt(2)/2, !. v(N, T, Q) : M is N1, v(M, P, R), Q is sqrt(2+2*R)/2, T is P*Q. ? viete(20, P). C (version with convergence check) #include <stdio.h> #include <math.h> int main() { double p = sqrt(2), q = p/2; while (p < 2) q *= (p = sqrt(2+p))/2; printf("%.17g\n", 2/q); } My own MiniC Clike language: int main() { float p, q; p = sqrt(2.0); q = p/2; while (p < 2.0) q *= (p = sqrt(2+p))/2; print 2/q; } $ ./minic viete.c $ java viete Java (version with convergence check): import java.lang.*; public class Viete { public static void main(String[] arg) { double p = Math.sqrt(2), q = p/2; while (p < 2) q *= (p = Math.sqrt(2+p))/2; System.out.println(2/q); } } Python (version with convergence check): from math import sqrt def viete(): p = sqrt(2) q = p/2 while p < 2: p = sqrt(2+p) q *= p/2 print(2/q) if __name__ == "__main__": viete() HPPL (version with convergence check): EXPORT viete() BEGIN LOCAL p = √2, q = p/2; REPEAT p := √(2+p); q := q*p/2; UNTIL p >= 2; RETURN 2/q; END; HP71B FORTH with MultiMod (forthcoming: where can I find the pdf manual?) Edit: Thanks to rprosperi's help to locate the missing FORTH manual, here is my HP71B FORTH program: FVARIABLE P FVARIABLE Q : VIETE 2. SQRT P STO 2. F/ Q STO BEGIN 2. P RCL F+ SQRT P STO 2. F/ Q RCL F* Q STO 2. P RCL X>=Y? UNTIL 2. Q RCL F/ F. ; VIETE HP71B FORTH significantly extends ANS FORTH. Couldn't be happier with MultiMod installed on my HP71B!  Rob Minor edit: fixed typo and added listbased functional versions of Husky and Haskell programs and HP71B FORTH. "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03212021, 07:00 PM
Post: #34




RE: [VA] SRC #009  Pi Day 2021 Special
(03212021 06:48 PM)robve Wrote: HP71B FORTH with MultiMod (forthcoming: where can I find the pdf manual?) Not sure which manual you're asking about. MultiMod Manual is here. The 71B Forth/Assembler Manual is part of the MoHPC Document Set (which I presume you have by now) but is also available here if you don't have that. Bob Prosperi 

03212021, 09:31 PM
Post: #35




RE: [VA] SRC #009  Pi Day 2021 Special
(03212021 07:00 PM)rprosperi Wrote:(03212021 06:48 PM)robve Wrote: HP71B FORTH with MultiMod (forthcoming: where can I find the pdf manual?)The 71B Forth/Assembler Manual is part of the MoHPC Document Set (which I presume you have by now) but is also available here if you don't have that. Great, many thanks! I had done some searches online but didn't find it before. I immediately ordered the MoHPC Document Set and used the link to the pdf to figure out how to use floating point in HP71B FORTH and its editor. It's very easy. Within minutes after skimming the documentation I had my program entered as SCREEN, running and spitting out \( pi \). The program: FVARIABLE P FVARIABLE Q : VIETE 2. SQRT P STO 2. F/ Q STO BEGIN 2. P RCL F+ SQRT P STO 2. F/ Q RCL F* Q STO 2. P RCL X>=Y? UNTIL 2. Q RCL F/ F. ; VIETE  Rob "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03222021, 04:06 PM
Post: #36




RE: [VA] SRC #009  Pi Day 2021 Special
Thanks, robve, nice to see the algorithm expressed in so many languages. For completeness here is an RPL version:
\<< 2 \v/ DUP 2 / DO SWAP 2 + \v/ SWAP OVER * 2 / UNTIL OVER 2 \>= END SWAP DROP 2 SWAP / \>> It gets the result 3.1415926536 after 19 iterations. Takes about 1.2 seconds on an HP28S. 

03222021, 05:33 PM
Post: #37




RE: [VA] SRC #009  Pi Day 2021 Special
(03212021 06:48 PM)robve Wrote: Proof (I've simplified this somewhat): A comment about the last (missing) step. Instead of using halfangle formula, cos(x/2)^2 = (1+cos(x))/2, then integrate, fold the integral. \(\displaystyle \int _a^b f(x)\;dx = \int _a^b {f(x) + f(a+bx) \over 2}\;dx \) \(\displaystyle 4 \int_0^{\pi \over 2} \cos^2 θ\;dθ = 2 \int_0^{\pi \over 2} (\cos^2 θ \;+\; \sin^2 θ)\;dθ = 2 \int_0^{\pi \over 2} 1\;dθ = \pi \) This is the same trick used in [VA] Short & Sweet Math Challenge #25. (02282021 02:18 AM)Valentin Albillo Wrote: My original solution for "Concoction the Third: Weird Integral" 

03232021, 10:20 PM
Post: #38




RE: [VA] SRC #009  Pi Day 2021 Special
Hi, all: Thanks to all of you who posted messages and/or comments to this SRC #009  \(\pi\) Day 2021 Special (namely JF Garnier, Gerson W. Barbosa, Albert Chan, PeterP, robve, Maximilian Hohmann, Ángel Martín and Massimo Gnerucci), for your interest and valuable contributions, much appreciated. Now, these are my original results and comments for the various parts of this SRC:
a. The root of this equation in [3,4] is indeed X = \(\pi\) ~ 3.14159265359. The base integral is: so that we have: I(\(\pi\)) = \(\pi\) \(\pi\)^{\(\pi\)}/\(\pi\)! = 15.9359953238 , I(e) = \(\pi\) e^{e}/e! = 11.1735566407 and this simple HP71B commandline expression will evaluate I(X) and \(\pi\) X^{X}/X! for any given X ≥ 0: >INPUT X @ INTEGRAL(0,PI,0,(SIN(IVAR)*EXP(IVAR/TAN(IVAR))/IVAR)^X);PI*X^X/GAMMA(X+1) ?1 3.14159265359 3.14159265359 { \(\pi\) } ?2 6.28318530717 6.2831853072 { 2 \(\pi\) } ?EXP(1) 11.1735566407 11.1735566407 { \(\pi\) e^{e} / e! } ?PI 15.9359953238 15.9359953238 { \(\pi\) \(\pi\)^{\(\pi\) }/ \(\pi\)! } Additional Notes:
b. The root of this equation in [3,4] is also X = \(\pi\) ~ 3.14159265359. This time the base integral is: which is a representation of the ubiquitous Lambert W function as a parametric definite integral, which though not new (there are other various integral representations) seems to me rather awesome nevertheless: a definite integral solves a transcendental equation, W_{0}(X) e^{W0(X) } = X Particularizing its value for X = 1 and exchanging W_{0}(1) and \(\pi\) we get: from where the equation is then obtained. This HP71B program lets us try different values of X and returns the difference between the LHS and the RHS, to see how it approaches 0 when X approaches \(\pi\): 1 DEF FNF(T)=LN(1+SIN(T)/T*EXP(T/TAN(T))) 2 DEF FNI(X)=INTEGRAL(0,X,0,FNF(IVAR)) 3 DESTROY ALL @ W=FNROOT(0,1,FVAR*EXP(FVAR)1) 4 INPUT X @ DISP FNI(X)/WX @ GOTO 4 [RUN] ? .1 .131342478087 ? .5 .63098466783 ? 1 1.1030252775 ? 1.5 1.27459847697 ? 2 1.08227897742 ? 2.5 .6404117652 ? 3 .14159265358 ? 3.14 .00159265358 ? 3.1415 .00009265358 ? PI .00000000001 For X > \(\pi\), it gives an error: ERR L1:LOG(neg). As expected, for X = \(\pi\) it returns ~ 0 but notice that for X ≥ 3 it returns ~ \(\pi\)  X. Using a graphing calculator to plot the above values into a continuous graph will show why. c. Though this seems to be a striking finite evaluation which gives \(\pi\) in terms of e in a much simpler and direct way than Euler's formula and without involving complex values, the magic is tarnished somewhat by the fact that e isn't really needed here, as the basic identity is: \(\pi\) = 4 * ( Arctan X  Arctan \(\frac{X  1}{X + 1}\) ) which is valid for all X > 1, as some of you explained and proved. That said, there are many interesting particular cases which you can use to trick your colleagues into believing you've discovered a brandnew, amazing evaluation for \(\pi\). For instance:  Using \(\pi\) itself (!) to get \(\pi\): \(\pi\) = 4 * (Arctan \(\pi\)  Arctan \(\frac{\pi  1}{\pi + 1}\)) >4*(ATN(PI)ATN((PI1)/(PI+1))) > 3.14159265359  Using the Golden Ratio \(\phi\) to get \(\pi\): \(\pi\) = 4 * (Arctan \(\phi\)  Arctan \(\frac{\phi  1}{\phi + 1}\)) >P=(1+SQR(5))/2 @ 4*(ATN(P)ATN((P1)/(P+1))) > 3.14159265360  Using the current year (2021) to get \(\pi\): \(\pi\) = 4 * (Arctan 2021  Arctan \(\frac{2020}{2022}\)) >4*(ATN(2021)ATN(2020/2022)) > 3.14159265358 You can also trick them by using your phone number, your birthday or any number personally related to you or the person being tricked ! Additional Notes:
d. The summation for even dimensions from 0 to infinity of the volumes enclosed by the respective ndimensional unit spheres (R = 1) is e^{\(\pi\)} ~ 23.1406926328 (aka Gelfond's constant) and thus its \(\pi\)th root is e, or conversely we could say that \(\pi\) is its natural logarithm. The \(\pi\)th root of the summation is readily obtained with the HP71B by executing this from the command line: >V=0 @ FOR D=0 TO 50 STEP 2 @ V=V+PI^(D/2)/FACT(D/2) @ NEXT D @ V^(1/PI) 2.71828182846 which agrees with e ~ 2.71828182846. For the 12digit HP71B we stop at dimension 50 because \(\pi\)^{25}/25! ~ 1.73.10^{13}, so adding further terms won't affect the result. Also, as I read for the first time in some Martin Gardner's book many decades ago, the volume enclosed by the ndimensional unit sphere tends to 0 with growing n, and reaches a maximum for a (fractional) dimension between 5 (Vol = 8 \(\pi\)^{2}/15 ~ 5.264) and 6 (Vol = \(\pi\)^{3}/6 ~ 5.168). See if you can find this unique dimension and the corresponding maximum volume. Additional Notes:
e. The song "\(\pi\)" by Kate Bush is indeed awesome, as is most of her music ("Cloudbusting", mentioned in this thread, certainly is, as is the video for it, almost a whole tragic movie told in a few minutes), and part of it appears in The Simpsons' 26thseason finale, "Mathlete's Feat", featuring about one minute of the song or so. \(\pi\) also appears in at least two other episodes which I watched at the time: in one of them, Prof. Frink unexpectedly says aloud that \(\pi\) is exactly equal to 3 as a way to get some much needed attention (he sure gets it !), and in another Homer and Marge are visiting some school for "Snotty Girls and Mama's Boys" (i.e., gifted children) and two of them are singing a handclapping song with lyrics they've concocted to help them remember a few digits of \(\pi\). There may be many more ... (episodes featuring \(\pi\), that is). Additional Notes:
f. First, those 31,415,926,535,897 decimal places were reported by Emma Haruka Iwao on 2019's \(\pi\) Day after 121 days of computation. Then, on January 2020 Tim Mullican computed 50 trillion digits in some 300 days, which if I'm not wrong it's the current world record as of 2021. The resulting string of digits passes all normality tests, where we consider a real number R to be normal in some base B if any string of N baseB digits appears in the baseB representation of R with frequency B^{N}, e.g: in base 10 every digit 09 appears 1/10th of the time, every 2digit sequence 0099 appears 1/100th of the time, etc. It can be proved that almost all real numbers are normal in any and all bases B, but rigurously proving that a "naturallyoccurring" real R (i.e., not artificially defined), say \(\pi\), is normal for even just one base B (say base 2 or 10) is excruciatingly difficult and not a single result is known so far, though the expectations are that all irrational numbers actually are, e.g. if you compute a trillion bits of \(\sqrt{2}\), you'll find about the same number of 0's and 1's and about the same number of 00's, 01's, 10's and 11's, etc. Not all is lost, however. If we can't (yet) prove than \(\pi\) is normal in base 2 or 10, say, we can try to estimate the credibility of the decision "\(\pi\) is not normal", which somewhat resembles probabilistic primality testing, where we can't rigurously prove that a certain number is prime in reasonable time but we may certify it as a probable prime if we can quickly prove that the probability of it being composite is extremely small. In the case of \(\pi\)'s normality, this has been done by analyzing the first 16 trillion bits of \(\pi\), and the result is that the decision "\(\pi\) is not normal" has credibility 4.3497.10^{3064}, which in my dictionary is akin to "impossible". Second, the bilingual joke about "Fear of number \(\pi\)" being called "Trescatorcephobia" is a pun on 3.14 being usually said as "tres catorce" in Spanish. Finally, re the "peer reviewed" papers demonstrating that \(\pi\)'s value is some algebraic number, I'm astonished that anyone would give them credit (let alone publish them) except for nonacademic reasons. Next category: papers succinctly proving Riemann's Hypothesis or Fermat's Last Theorem in a few pages, almost casually. Additional Notes:
That's all for now. Again, thanks for your interest and participation in this SRC, glad you enjoyed it ! Best regards. V. All My Articles & other Materials here: Valentin Albillo's HP Collection 

03242021, 05:53 PM
(This post was last modified: 03242021 08:57 PM by robve.)
Post: #39




RE: [VA] SRC #009  Pi Day 2021 Special
(03232021 10:20 PM)Valentin Albillo Wrote: Again, thanks for your interest and participation in this SRC, glad you enjoyed it ! Thank you for posting the pi day special and for the conclusion! I moved the remainder of the long post as a separate topic on Adaptive Simpson and Romberg methods, as suggested: https://www.hpmuseum.org/forum/thread16523.html  Rob "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

03262021, 03:48 AM
(This post was last modified: 03262021 01:25 PM by robve.)
Post: #40




RE: [VA] SRC #009  Pi Day 2021 Special
b. answer to what this program computes:
10 P=1 20 FOR I=2 TO 1000 STEP 2 30 P=P*I/(I1) 40 DISP 2*P*P/(I+1) 50 NEXT I As Albert correctly points out and explains in his reply, the program computes Wallis Product. $$ \frac{\pi}{2} = \prod_{n=1}^\infty \frac{4n^2}{4n^21} = \frac{2\cdot 2\cdot 4\cdot 4\cdot 6\cdot 6\cdot 8\cdot 8\cdots}{1\cdot 3\cdot 3\cdot 5\cdot 5\cdot 7\cdot 7\cdot 9\cdots} $$ The proof is "typically featured in modern calculus textbooks" according to this Wikipedia article. Another way to express the Wallis Product in compact form with factorials: $$ \frac{\pi}{2} \approx \begin{cases} \displaystyle\frac{n\,(n1)!!^2}{n!!^2} & \text{if $n\ge 3$ is odd} \\ \displaystyle\frac{n!!^2}{(n+1)\,(n1)!!^2} & \text{if $n\ge 2$ is even} \end{cases} $$ where \( n!! \) is the double factorial. Of course, we would never want our programs to use such factorials in the numerator and denominator, because of overflow and loss of floating point precision. There are many ways this product can be implemented in code. The code shown above instantiates this product's numerator directly as 2*P*P to make this relation of this code to the product form more readily apparent. Wallis product very slowly converges to \( \pi \), much slower than Viète's formula and other hiistorical infinite series formulas for \( \pi \). It is not a practical algorithm to compute \( \pi \) to many decimals. Efficient algorithms exist to compute \( \pi \) in many decimals, such as the following short C program (160 characters when compacted) written by Dik T. Winter at CWI (Centrum voor Wiskunde en Informatica  the national research institute for mathematics and computer science) Amsterdam to compute \( \pi \) up to 800 decimal digits: int a = 10000, b, c = 2800, d, e, f[2801], g; main() { for (; bc; ) f[ b++] = a/5; for (; d=0, g = c*2; c=14, printf("%.4d", e+d/a), e = d % a) for (b = c; d += f[ b]*a, f[ b] = d % g, d /= g, b; d *= b) ; } 31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185 Other programs to compute the Wallis Product Haskell (functional PL): main = putStrLn (show (wallis 1000)) where wallis n = (2*p^2/(n+1)) where (_,p) = until done next (n,1.0) next (k,p) = (k2,p*k/(k1)) done (k,p) = (k==0) C (imperative PL): int main() { int i, n = 1000; // must be even double p = 1; for (i = 2; i <= n; i += 2) p *= (double)i/(i1); printf("%g\n", 2*p*p/(n1)); } Prolog (logic PL): wallis(W) : prod(1000, P), W is 2*P*P/1001. prod(0, 1) : !. prod(K, Q) : M is K2, prod(M, P), Q is P*K/(K1).  Rob "I count on old friends"  HP 71B,PrimeTi VOY200,Nspire CXII CASCasio fxCG50...Sharp PCG850,E500,2500,1500,14xx,13xx,12xx... 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)