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

[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