Lest you think I am an expert in continued fractions, which definitely I AM NOT,
I will disclose my reply to someone who wanted to know how I had found the above formula.
Please do not be too disappointed at my "method".
-------
If you compute the Wallis product to 10000 digits,
3.140807746030394560485868201553507
and take its difference to pi,
0.000078534907871191484085332565629
you’ll notice the first significant digits are the same of those of pi/4,
0.7853981633
If you divide pi/4 by that difference you’ll get
10000.62500468726562573348383657523
Multiply this by 8 (because the fractional part is .625, whose multiplication by 8 will lead to a near-integer) and get
80005.00003749812500586787069260184
That’s the ‘8*n + 5’ part in the formula.
Take the fractional part of that,
0.00003749812500586787069260184
Invert it:
26668.00006249687481715340569058372
Multiply by 3:
80004.00018749062445146021707175116
That’s the first term of the continued fraction:
3/(8*n + 4)
Again, take the fractional part of that and invert it:
5333.600028938469868026787967169474
If you multiply this by 3 you’ll notice you’ll have to multiply it again by 5 so that the result is a near integer. So, let’s multiply it by 15:
80004.0004340770480204018195075421
Hence the next term of the continued fraction:
15/(8*n + 4)
So far we have
(pi/4)/(pi - Wn) = 8n + 5 + 3/(8n + 4 + 15/(8n + 4 + ...))
If you keep doing this for the next terms you’ll soon run out of precision to get meaningful results. So, now that you know the first
two terms, start over with W(1000) instead of W(10000). Thus, you can go up to 6 terms, which suffices for unveiling the pattern:
The first part before the CF starts is 8*n + 5, then the numerators are 3, 15, 35, 63, 99... and the denominators are 8*n + 4 + ...
Isolate pi and there you are.
-------
HP 50g program
(263 bytes, CHKSUM = # 9159h)
Code:
« PUSH RAD -105 CF -3 CF DUP 10 * 21 / CEIL 8 OVER * 4 + SWAP 0 1 ROT OVER
FOR i i DUP + SQ * LASTARG NIP 1 - DUP UNROT / 4 ROLLD PICK3 ROT + / ROT -1
STEP UNROT + 1 - 2 4 ROT / + * EXPAND FXND DUP SIZE R→I ALOG OVER - PICK3 *
SWAP IQUOT + →STR DUP HEAD 0 I→R →STR TAIL + SWAP TAIL + 1 ROT 2 + SUB POP
»
Time to compute 314 decimal digits of \(\pi\) on the real 50g: 2979.4707 seconds (14.9503 seconds on the emulator).
The formula is linearly convergent (about 2.1 digits per iteration).
Blame it on the slowness of EXPAND to process one algebraic expression containing 150 product and 151 continued fraction terms.
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631