Programming Exercise (HP-15C, 15C LE - and others)

[Gerson W. Barbosa]

(04-06-2014 06:24 AM)Thomas Klemm Wrote:  Now we have 3 different ways to calculate the same thing:

1. Borwein's formula using tangent numbers
2. Convergence acceleration using Euler's transformation
3. Gerson's method using continued fractions
   

Do you have a proof of your formula or is it still a conjecture? Congratulations to your discovery! Your method appears to be more efficient than the other two.

Isn't this the same method used by Valentin in his original problem, only in continued fraction format? Anyway, the formula appears to hold, even when we have only two terms of the series, albeit converging very slowly in this case:

\[\ln (2)= 1-\frac{1}{2}+\frac{1}{3}-\frac{1}{4}+\cdots +\frac{1}{n-1}-\frac{1}{n}+\frac{1}{2n+1+\frac{1^{2}}{2n+1+\frac{2^{2}}{2n+1+\frac{3^{2}}{2n+1+​\frac{4^{2}}{2n+1+... }}}}}\]

...ln(2) = 0.6931471805599453094172321214581766
0001e000 2 0.6923076923076923076923076923076923
0001e001 2 0.6931477851289172043889025021100493
0001e002 2 0.6931471805715443422837621214224286
0001e003 2 0.6931471805599454334834722951537730
0001e004 2 0.6931471805599453094184811843330782
0001e005 2 0.6931471805599453094172321339572391
0001e006 2 0.6931471805599453094172321214583016
0001e007 2 0.6931471805599453094172321214581766


Here is a similar formula for π/4:

\[\frac{\pi }{4}= 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots +\frac{1}{2n-3}-\frac{1}{2n-1}+\frac{1}{4n+\frac{1^{2}}{n+\frac{2^{2}}{4n+\frac{3^{2}}{n+\frac{4^{2}}{4n+...​ }}}}}\]

.......π = 3.141592653589793238462643383279503
00000009 2 3.141574838064079287182344939989334
00000099 2 3.141592651424974945646180479799372
00000999 2 3.141592653589572786959492672266442
00009999 2 3.141592653589793216377737847078192
00099999 2 3.141592653589793238460434495123538
00999999 2 3.141592653589793238462643162386710
09999999 2 3.141592653589793238462643383257414
99999999 2 3.141592653589793238462643383279500


About 2d/3 terms of the series and of the continued fraction, where d is the number of desired correct decimal places, appear to suffice:

0000010 10 0.6931471805599454035014297564573935
0000022 22 0.6931471805599453094172321214581766

0000011 10 3.141592653589793130085240933000310
0000021 22 3.141592653589793238462643383279502

This is the program to compute π:

Code:


001 LBL A
002 STO 00
003 STO 01
004 STO+ X
005 STO+ 01
006 STO+ 01
007 #000
008 DEC Y
009 RCL Y
010 DEC X
011 DEC X
012 RCL* Z
013 z<> L
014 1/x
015 +
016 RCL+ L
017 DSE Y
018 BACK 010
019 x<> Y
020 RCL+ 00
021 RCL Z
022 x²
023 x<> Y
024 /
025 RCL 01
026 RCL- 00
027 x<> 00
028 Rv
029 DSE Z
030 BACK 010
031 RCL+ 00
032 1/x
033 + 
034 STO+ X
035 STO+ X
036 END

Again, the second parameter has to be even. Also, due to a program limitation the first parameter has to be odd.

Cheers,

Gerson.

[Thomas Klemm]

(04-15-2014 04:41 PM)Gerson W. Barbosa Wrote:  Isn't this the same method used by Valentin in his original problem, only in continued fraction format?

It appears to be the same but I was just curious how Euler's transformation is related to your continued fraction.

Quote:Here is a similar formula for π/4:

\[\frac{\pi }{4}= 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots +\frac{1}{2n-3}-\frac{1}{2n-1}+\frac{1}{4n+\frac{1^{2}}{n+\frac{2^{2}}{4n+\frac{3^{2}}{n+\frac{4^{2}}{4n+...​ }}}}}\]

Very nice! Thanks for being so persistent.

Cheers
Thomas

[Gerson W. Barbosa]

(04-15-2014 06:25 PM)Thomas Klemm Wrote:  

(04-15-2014 04:41 PM)Gerson W. Barbosa Wrote:  Isn't this the same method used by Valentin in his original problem, only in continued fraction format?

It appears to be the same but I was just curious how Euler's transformation is related to your continued fraction.

Quote:Here is a similar formula for π/4:

\[\frac{\pi }{4}= 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots +\frac{1}{2n-3}-\frac{1}{2n-1}+\frac{1}{4n+\frac{1^{2}}{n+\frac{2^{2}}{4n+\frac{3^{2}}{n+\frac{4^{2}}{4n+...​ }}}}}\]

Very nice! Thanks for being so persistent.

If we take only one term of the series then the continued fraction becomes

\[\frac \pi 4 = 1-\cfrac{1}{4+\cfrac{1^2}{1+\cfrac{2^2}{4+\cfrac{3^2}{1+\cfrac{4^2}{4+\cfrac{5^2}{​1+\ddots}}}}}}\]

which resembles Brouncker's formula:

\[\frac 4 \pi = 1+\cfrac{1^2}{2+\cfrac{3^2}{2+\cfrac{5^2}{2+\cfrac{7^2}{2+\cfrac{9^2}{2+\ddots}}​}}}\]

The latter, equivalent to the Gregory series, doesn't include even squares and converges more slowly:

...........n .....result................n .....result

000000000001 3.00000000000000000000000001 2.66666666666667
000000000010 3.14513671429771000000000010 3.23231580940559
000000000100 3.14163153155977000000000100 3.15149340107099
000000001000 3.14159304589625000000001000 3.14259165433954
000000010000 3.14159265751639000000010000 3.14169264359054
000000100000 3.14159265362906000000100000 3.14160265348979
000001000000 3.14159265359019000001000000 3.14159365358879
000010000000 3.14159265358980000010000000 3.14159275358978
00000000000000000000000000000000100000000 3.14159266358979
00000000000000000000000000000001000000000 3.14159265458979


But what has intrigued me is this part of an article on William Brouncker's biography, which reads:



"This result, written up in around ten pages, was added by Wallis to his treatise Arithmetica Infinitorum and probably first discovered by Brouncker in 1654. Wallis told Huygens of this result and Huygens expressed strong doubts that it was true. However after Brouncker correctly computed the first 10 places in the decimal expansion of π using his continued fraction expansion, Huygens accepted the result."

From the table above, we know this would have required the computation of 10^9 terms of the continued fraction, an impossible task to do by hand. However, one hundred and twenty terms could well have been computed in a week or less. Then a simple weighted mean, the weights being two successive number of terms, might have given about 10 decimal places, provided, of course, he was able to explain why this apparently works - I am not :-)

Code:


{ For up to 18-digit figures, compile in Borland TurboBCD }
Program CF_PIWM;
var  p1, p2, pi, t: real;

Function Brouncker(t: real): real;
var i, r: real;
begin
  i := 2*t - 1;
  r := 0;
  repeat
    r := Sqr(i)/(r + 2);
    i := i - 2
  until i < 0;
  r := r + 1;
  Brouncker := 4/r
end;

begin
  ClrScr;
  WriteLN('   n          p1 = p(n-1)            p2 = p(n)      pi~((n-1)*p1+n*p2)/(2*n-1)',#10);
  t:=10;
  repeat
    p1 := Brouncker(t - 1);
    p2 := Brouncker(t);
    pi := (p1*t + p2*(t + 1))/(2*t + 1);
    Writeln(t:5:0,'     ',p1:18:16,'     ',p2:18:16,'     ',pi:18:16);
    t := t + 10
  until t > 150
end.

...n .........p1 = p(n-1)............p2 = p(n)......pi~((n-1)*p1+n*p2)/(2*n-1)

00010.....3.0418396189294022.....3.2323158094055927.....3.1416128615597877
00020.....3.0916238066678386.....3.1891847822775947.....3.1415940624679576
00030.....3.1082685666989461.....3.1738423371907494.....3.1415929418669117
00040.....3.1165965567938323.....3.1659792728432150.....3.1415927463990754
00050.....3.1215946525910105.....3.1611986129870501.....3.1415926919989117
00060.....3.1249271439289964.....3.1579849951686658.....3.1415926722399041
00070.....3.1273076679812339.....3.1556764623074750.....3.1415926637057950
00080.....3.1290931417757212.....3.1539378622726156.....3.1415926595412395
00090.....3.1304818853613080.....3.1525813328751202.....3.1415926573157661
00100.....3.1315929035585528.....3.1514934010709906.....3.1415926560399270
00110.....3.1325019323081855.....3.1506014798194977.....3.1415926552663559
00120.....3.1332594649198298.....3.1498569752932738.....3.1415926547753764
00130.....3.1339004596806044.....3.1492261301786887.....3.1415926544516736
00140.....3.1344498875489983.....3.1486847629938381.....3.1415926542312845
00150.....3.1349260609930860.....3.1482150975379365.....3.1415926540770475

Madhava's correction terms (14th century), might as well have been used. This appears to be quite an ancient game!

Regarding the generic continued fraction for Ln(2), if no term of the original series is taken, that is, if n = 0, then the resulting continued fraction will be equivalent to the original series itself. A WP 34S program takes about 20 seconds to evaluate all 10000 terms though.

Cheers,

Gerson.

Edited to fix a couple of typos.

P.S.: As shown above, the very slowly convergent Brouncker's formula (equivalent to the Gregory series) can be sped up by a simple weighted mean, although not nearly as effectively as the other methods. Anyway, 3000 terms would give 15 correct digits. The plain series or continued fraction would require 10^14 terms to yield the same accuracy.

Code:


10  DEFDBL A-Z
15  CLS: KEY OFF
20  T=1000
25  WHILE T<3001
30    N=T-1: GOSUB 100: P1=B
35    N=T: GOSUB 100: P2=B
40    PI=(P1*T+P2*(T+1))/(2*T+1)
45    PRINT USING"####";N;:PRINT"   ";
50    PRINT USING"#.###############";P1;: PRINT"   ";
55    PRINT USING"#.###############";P2;: PRINT"   ";
60    PRINT USING"#.###############";PI
65    T=T+125
70  WEND
90  END
100 I=2*N-1
110 R=0
120 WHILE I>0
130   R=I*I/(R+2)
140   I=I-2
150 WEND
160 R=R+1
170 B=4/R
180 RETURN

This is a GW BASIC version of the Turbo Pascal code above.

RUN
1000...3.140592653839793...3.142591654339543...3.141592653590043
1125...3.142481542303099...3.140704554297768...3.141592653589638
1250...3.140792653717793...3.142392013973691...3.141592653589896
1375...3.142319926220898...3.140865909499706...3.141592653589724
1500...3.140925986997201...3.142258876034188...3.141592653589843
1625...3.142208038146917...3.140977647497886...3.141592653589757
1750...3.141021225065012...3.142163755770525...3.141592653589820
1875...3.142125986885201...3.141059604587147...3.141592653589773
2000...3.141092653621043...3.142092403683528...3.141592653589809
2125...3.142063241799034...3.141122286729639...3.141592653589781
2250...3.141148209167297...3.142036900569207...3.141592653589803
2375...3.142013706202711...3.141171778187556...3.141592653589785
2500...3.141192653605793...3.141992493637787...3.141592653589800
2625...3.141973605956924...3.141211846292098...3.141592653589788
2750...3.141229017238178...3.141956157758083...3.141592653589798
2875...3.141940479666230...3.141244948454266...3.141592653589790
3000...3.141259320265719...3.141925875839790...3.141592653589796
Ok 
PRINT TAB(47) 4*ATN(1#)
...............................................3.141592653589793
Ok