HP-42S (Free42, DM-42) mini-challenge – Wallis Product

[Gerson W. Barbosa]

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.141592653589793238462643383279502884197169399375105820974944592307816406286208​99862803482534211706798214808651328230664709384460955058223172535940812848111745​02841027019385211055596446229489549303819644288109756659334461284756482337867831​6527120190914564856692346034861045432664821339360726024914127372458700660631