Π day

[Thomas Klemm]

(03-15-2022 04:56 PM)robve Wrote:  So I am not sure if this is the best choice of parameterization for this particular alternating series.

I only used what I wrote in the linked post, which was copied verbatim from Valentin's article.
But I agree with you. This series is not slowly converging like this one.
Therefore it's probably not that easy to decide if any form of acceleration is worth the effort.
That's why I gave an alternative solution for the HP-42S.

(03-15-2022 12:35 AM)robve Wrote:  In addition, I observed that summing the terms in reverse order may sometimes improve accuracy.

That's the point of Wirth's exercise.

(03-15-2022 12:35 AM)robve Wrote:  That's far more CPU power than the N=17 steps for the Sharp series to converge to 10 decimal places.

If you really care about CPU cycles you may want to avoid the expensive \(3^n\) calculation which usually uses \(\log\) and \(\exp\) under the hood.
For integer exponents this can be done more efficiently.
But I don't know if that is used by your Sharp calculator.

This can be achieved by rearranging the series like so:

\(
\begin{align}
1-{\frac {1}{3\cdot 3}}+{\frac {1}{5\cdot 3^{2}}}-{\frac {1}{7\cdot 3^{3}}}+\cdots =
1 - \frac{\frac{1}{3}-\frac{\frac{1}{5}-\frac{\frac{1}{7}-\frac{\frac{1}{9}-\cdots}{3}}{3}}{3}}{3}
\end{align}
\)

The calculation starts from inside out or rather top down like this:

Code:

def madhava(n):
    k = 2 * n + 1
    s = 0
    while k > 0:
        s = 1 / k - s / 3
        k -= 2
    return s

Cheers
Thomas

---

Quote:"I can count on my friends"

Is this count your friend?

[robve]

(03-15-2022 07:55 PM)Thomas Klemm Wrote:  

(03-15-2022 12:35 AM)robve Wrote:  In addition, I observed that summing the terms in reverse order may sometimes improve accuracy.

That's the point of Wirth's exercise.

Good point. I tried to convey this subtly as a minor contribution to that other thread by just restating a known fact that summing a series with monotonically decreasing magnitude should always be performed in reverse. It's easy to see why when you look at floating point addition that right shifts the mantissa of one of the operands until the exponents align. You want to start with the smallest terms, those who's exponents align most closely to the (initially) small sum.

(03-15-2022 07:55 PM)Thomas Klemm Wrote:  

(03-15-2022 12:35 AM)robve Wrote:  That's far more CPU power than the N=17 steps for the Sharp series to converge to 10 decimal places.

If you really care about CPU cycles you may want to avoid the expensive \(3^n\) calculation which usually uses \(\log\) and \(\exp\) under the hood.

To clarify, a better way is to perform full "strength reduction" as compilers do to optimize loops. Even better, for equidistant grids in one or more dimensions, functions can be more effectively optimized using the Chains of recurrences algebra, see https://dl.acm.org/doi/abs/10.1145/19034...NSgF9Z3eNA and our article https://dl.acm.org/doi/abs/10.1145/13755...aRdGHFiMLQ (our earlier work on the inverse CR algebra is used by the popular GCC compiler to optimize loops.)

After strength reduction we end up with one division, two multiplications, two additions and one negation per loop iteration:

S=0
P=1
X=1
Y=1
FOR I=0 TO N
  S=S+P/(X*Y)
  P=-P
  X=X*3
  Y=Y+2
NEXT I

Note that the first iteration can be stripped:

S=1
P=-1
X=3
Y=3
FOR I=1 TO N
  S=S+P/(X*Y)
  P=-P
  X=X*3
  Y=Y+2
NEXT I

and "operator strength reduction" then gives:

S=1
P=-1
X=3
Y=3
FOR I=1 TO N
  S=S+P/(X*Y)
  P=-P
  X=X+X+X
  Y=Y+2
NEXT I

However, and this is an important point, there is a limit to how aggressive we can be with strength reduction (via Chains of Recurrences) to optimize floating point operations, because of potentially introducing small errors that accumulate over large grids (not a problem when N is small). See https://dl.acm.org/doi/abs/10.1145/38410...6bSQ-bA_eg Error accumulation is not a problem for whole numbers, so \( 3^i \) is always safe to optimize as shown above.

PS. Division is typically the most expensive operation, which should be avoided at all cost, when possible.

PPS. Summing terms in the order of increasing magnitude can be important, but methods such as Kahan summation work generally best. Also binary tree summation is generally better (the proof is in my lecture notes on parallel computing).

- Rob

[Thomas Klemm]

(03-15-2022 08:49 PM)robve Wrote:  functions can be more effectively optimized using the Chains of recurrences algebra

That sounds interesting.
Unfortunately, I do not have access to these papers.
Would it be possible to give this?

Many thanks in advance
Thomas

[Thomas Klemm]

(03-14-2022 03:35 AM)robve Wrote:  What's your favorite pi series?

Not a series but I like Vieta's Formula for Pi since only 2 is used:

\(
\begin{align}
\pi = 2
\times \frac{2}{\sqrt{2}}
\times \frac{2}{\sqrt{2 + \sqrt{2}}}
\times \frac{2}{\sqrt{2 + \sqrt{2 + \sqrt{2}}}}
\times \frac{2}{\sqrt{2 + \sqrt{2 + \sqrt{2 + \sqrt{2}}}}}
\times \cdots
\end{align}
\)

Code:

00 { 30-Byte Prgm }
01▸LBL "VIETA"
02 2
03 ENTER
04 ENTER
05 STO 00
06 0
07▸LBL 00
08 RCL+ ST T
09 SQRT
10 ÷
11 STO× 00
12 X<> ST L
13 VIEW 00
14 GTO 00
15 END

R00=2.82842712475
x: 1.41421356237

R00=3.06146745892
x: 1.84775906502

R00=3.12144515226
x: 1.96157056081

R00=3.13654849055
x: 1.99036945334

R00=3.14033115695
x: 1.99759091241

R00=3.14127725093
x: 1.99939763739

R00=3.14151380114
x: 1.99984940368

R00=3.14157294037
x: 1.99996235057

R00=3.14158772528
x: 1.99999058762

R00=3.14159142151
x: 1.99999764690

R00=3.14159234557
x: 1.99999941173

R00=3.14159257658
x: 1.99999985293

R00=3.14159263434
x: 1.99999996323

R00=3.14159264878
x: 1.99999999081

R00=3.14159265239
x: 1.99999999770

R00=3.14159265329
x: 1.99999999943

R00=3.14159265351
x: 1.99999999986

R00=3.14159265357
x: 1.99999999996

R00=3.14159265359
x: 1.99999999999

R00=3.14159265359
x: 2.00000000000

[Gerson W. Barbosa]

(03-17-2022 03:54 AM)Thomas Klemm Wrote:  

Code:

00 { 30-Byte Prgm }
01▸LBL "VIETA"
02 2
03 ENTER
04 ENTER
05 STO 00
06 0
07▸LBL 00
08 RCL+ ST T
09 SQRT
10 ÷
11 STO× 00
12 X<> ST L
13 VIEW 00
14 GTO 00
15 END

We can save at least one byte, one step and one register by using Katie Wasserman’s idea here:

Code:



00 { 29-Byte Prgm }
01▸LBL "VIETA"
02 2
03 RCL X
04 0
05▸LBL 02
06 +
07 SQRT
08 STO÷ ST Y
09 2
10 STO× ST Z
11 VIEW ST Z
12 X≠Y?
13 GTO 02
14 END

[Thomas Klemm]

(03-17-2022 11:39 AM)Gerson W. Barbosa Wrote:  We can save at least two bytes (…)

Another byte beats the dust:

Code:

00 { 26-Byte Prgm }
01▸LBL "VIETA"
02 2
03 0
04▸LBL 00
05 +
06 SQRT
07 STO÷ ST Y
08 2
09 STO× ST Z
10 VIEW ST Z
11 GTO 00
12 END

I also removed the check because the VIEW command stops when flag 21 is set.
Thus I simply hit the R/S button as long as I want.

---

Just noticed your edit: I don't think that the RCL X command is needed.

[Gerson W. Barbosa]

Thomas Klemm dateline=' Wrote:  I also removed the check because the VIEW command stops when flag 21 is set.
Thus I simply hit the R/S button as long as I want.

---

Just noticed your edit: I don't think that the RCL X command is needed.

I’d replaced RCL X with ENTER. It worked once only because there was a 2 on the stack I hadn’t noticed. Then I went back to the initial version. But you’re right, better use the R/S to stop, that’s what that key is for.

——

PS.: RCL X is necessary. Try running the program after clearing the stack.

[DavidM]

Back in March of 2004, Olivier De Smet posted this fairly compact RPL version of a pi digits program (see post #17):

Code:

%%HP: T(3)A(R)F(.);
\<< -105. CF \-> N
  \<< 280000 N ALOG * DUP 'R' STO 2 \-> y
    \<<
      DO y * y 1 + IQUOT 50 IQUOT DUP 'R' STO+ 2 'y' STO+
      UNTIL DUP 0 ==
      END DROP
    \>> 30336 N ALOG * DUP 'R' STO+ 2 \-> y
    \<<
      DO 9 * y * 6250 IQUOT y 1 + IQUOT DUP 'R' STO+ 2 'y' STO+
      UNTIL DUP 0 ==
      END DROP
    \>>
  \>> R
\>>

It uses exact integers for intermediate computations, thus providing digit counts in excess of the inherent approximate limitations. He attributed it as an RPL translation of a program originally posted by Katie Wasserman which was written for the HP 32SII. Katie, in turn, cited the algorithm as essentially being the Euler/Machin approach (Pi = 20 * arctan(1/7) + 8 * arctan(3/79)).

As soon as I saw Olivier's version some time last year, I wanted to try a SysRPL implementation based on his code that would avoid the inherent performance penalty of IQUOT. This seemed like a good time to try it.

This is essentially Olivier's program, but coded in SysRPL and providing a result rounded to the specified number of digits (Olivier's provides 6 additional digits). The code will compile on a 49g+/50g with the standard HP extable installed (use 256.06 MENU ASM). I've used some DCs and DEFINEs to simplify the code listing a bit.

Code:

!NO CODE
!ASM
DC Z0_      273AB
DC Z1_      273B6
DC Z2_      273C2
DC Z9_      27416
!RPL
DEFINE      z*                FPTR2 ^QMul
DEFINE      z+                FPTR2 ^QAdd
DEFINE      zInc1             Z1_ z+
DEFINE      zInc2             Z2_ z+
DEFINE      zIntDiv           FPTR2 ^ZDIVext DROP
DEFINE      zDupZero?         FPTR2 ^DupQIsZero?
DEFINE      zALOG             FPTR 001 00E1
::
   CK1NoBlame
   CK&DISPATCH1
   real ::
      ( convert digit count to BINT )
      COERCE
      
      ( "bad argument value" if digit count is 0 )
      DUP#0=case SETSIZEERR
      
      ( setup locals )
      FPTR2 ^#>Z zALOG                       ( mag - antilog of digit count )
      Z2_                                    ( num - initially 2 )
      OVER ZINT 280000 z* DUP4UNROLL         ( result - initially mag*280000 )
      {{ result num mag }}
      
      ( initial result already on stack from above )
      BEGIN
         num z*
         num zInc1 zIntDiv
         ZINT 50 zIntDiv
         DUP result z+ !result
         num zInc2 !num
         zDupZero?   
      UNTIL
      DROP
      
      ZINT 30336 mag z*
      DUP result z+ !result
      Z2_ !num
      
      BEGIN
         Z9_ num z* z*
         ZINT 6250 zIntDiv
         num zInc1 zIntDiv
         DUP result z+ !result
         num zInc2 !num
         zDupZero?   
      UNTIL
      DROP

      ( retrieve result, abandon LAMs )
      1GETABND
      
      ( round the result to the number of digits specified )
      ZINT 500000 z+
      ZINT 1000000 zIntDiv
   ;
;
@

Performance improved, as expected. Generating 1000 digits with Olivier's original RPL program takes about 231 seconds (you actually get 1006 digits, but the last 4 are not accurate). Generating 1006 digits using the SysRPL version takes about 111 seconds, and all resulting digits appear to be correct (at least they match Wolfram Alpha's output for 1006 digits).

The usual advantages of SysRPL apply here, but a significant part of the improvement is in bypassing the overhead of IQUOT and going straight to the internal exact integer division code.

I'm sure NewRPL on the 50g and Python on the Prime would perform much better, but I thought this was reasonable performance for the task at hand on the original 50g platform.

[robve]

(03-17-2022 03:40 AM)Thomas Klemm Wrote:  

(03-15-2022 08:49 PM)robve Wrote:  functions can be more effectively optimized using the Chains of recurrences algebra

That sounds interesting.
Unfortunately, I do not have access to these papers.
Would it be possible to give this?

Many thanks in advance
Thomas

This link shows related articles on this subject, including PDFs:
https://scholar.google.com/scholar?q=rel...as_sdt=0,6

- Rob

[Xorand]

Back before we had a calculator with a Pi key (let alone programmability), we used the old standby 355 / 113. It gives a fairly decent approximation (to at least 5 digits) and can be handily entered on a four-banger calculator.

On the other hand, it was pretty straightforward to memorize pi out to about 8 digits of precision.

[Steve Simpkin]

(03-18-2022 03:06 AM)Xorand Wrote:  Back before we had a calculator with a Pi key (let alone programmability), we used the old standby 355 / 113. It gives a fairly decent approximation (to at least 5 digits) and can be handily entered on a four-banger calculator.

On the other hand, it was pretty straightforward to memorize pi out to about 8 digits of precision.

Those with Sinclair Scientific calculators didn’t need a Pi key (there was not enough ROM available to store it anyway) or need to memorize the value. It was printed on the front of the calculator. Of course it was only shown to 5 decimal places which was generally more accurate than the calculator. If you haven’t seen it be sure to check out Ken Shirriff’s blog post on reverse Engineering this model.
http://files.righto.com/calculator/sincl...lator.html

[Ángel Martin]

(03-17-2022 02:10 PM)Gerson W. Barbosa Wrote:  

Thomas Klemm dateline=' Wrote:  I also removed the check because the VIEW command stops when flag 21 is set.
Thus I simply hit the R/S button as long as I want.

---

Just noticed your edit: I don't think that the RCL X command is needed.

I’d replaced RCL X with ENTER. It worked once only because there was a 2 on the stack I hadn’t noticed. Then I went back to the initial version. But you’re right, better use the R/S to stop, that’s what that key is for.

——

PS.: RCL X is necessary. Try running the program after clearing the stack.

You guys have come up with the perfect VIETA User code, so there's no point for my trying to improve on it... so here comes an MCODE version instead.

In manual mode (not running a program) it will display the values afforded by the successive iterations, so you'll see the progress being made...

Will be included in the PIE_ROM as well.

Cheers,
ÁM

[MeindertKuipers]

This discussion reminds me about an article from NASA's JPL. It explains that only 15 decimals of pi are really needed to calculate the circumference of the orbit of a distance where Voyager (in 2016) is down to an error of about 1.5 inch.

https://www.jpl.nasa.gov/edu/news/2016/3...ally-need/

[Ángel Martin]

(03-18-2022 10:48 AM)MeindertKuipers Wrote:  This discussion reminds me about an article from NASA's JPL. It explains that only 15 decimals of pi are really needed to calculate the circumference of the orbit of a distance where Voyager (in 2016) is down to an error of about 1.5 inch.

https://www.jpl.nasa.gov/edu/news/2016/3...ally-need/

Very nice, a real down-to-earth (pun intended) perspective on the issue!

Of course that doesn't preclude us from attempting math explorations, never mind their actual usefulness in the engineering world, though.

[Frido Bohn]

(03-18-2022 09:07 AM)Ángel Martin Wrote:  You guys have come up with the perfect VIETA User code, so there's no point for my trying to improve on it... so here comes an MCODE version instead.

In manual mode (not running a program) it will display the values afforded by the successive iterations, so you'll see the progress being made...

Will be included in the PIE_ROM as well.

Hi Ángel,

Your MCODE needs Library 4 to work properly, right?
After doing some iterations it stops. After pressing R/S it asks 'RadIUS'. Enter '1' and you get a 'Pi' that isn't different from the key-bound Pi.

Great achievement!
Thanks!

[Ángel Martin]

(03-19-2022 09:45 AM)Frido Bohn Wrote:  Hi Ángel,

Your MCODE needs Library 4 to work properly, right?
After doing some iterations it stops. After pressing R/S it asks 'RadIUS'. Enter '1' and you get a 'Pi' that isn't different from the key-bound Pi.

Hi Frido, yes it does to present the "RUNNING..." message - sorry about that but it's kind of my standard these days/ - and as part of the PIE ROM (forthcoming) this is a worth-having dependency for may other reasons as well.

I have no idea about those strange prompts/questions you're getting, the VIETA function does not prompt for any parameter at all. It'ss likely that when you press R/S the calculator defaults to your current FOCAL program in memory, which is generating the prompts...

Best,
ÁM

[Ángel Martin]

Heresy or new truth?

Here's an interesting reading for y'all:
https://www.theverge.com/tldr/2018/3/14/...iday-truth

[Frido Bohn]

(03-19-2022 11:17 AM)Ángel Martin Wrote:  I have no idea about those strange prompts/questions you're getting, the VIETA function does not prompt for any parameter at all. It'ss likely that when you press R/S the calculator defaults to your current FOCAL program in memory, which is generating the prompts...

Hi Ángel,

You are right, the R/S instruction indeed jumps into the current FOCAL program step. Just by coincidence it had to do something with Pi. Sorry for the confusion.

KR
Frido

[Frido Bohn]

(03-18-2022 09:07 AM)Ángel Martin Wrote:  You guys have come up with the perfect VIETA User code, so there's no point for my trying to improve on it... so here comes an MCODE version instead.

Hi Ángel,

The first 3 instructions of your code (391, 08C, 079) point to Address 879 of the same port where a routine 'INITLZ' is apparently located. Of course, by simply copying the 'Vieta' MCODE into my ROM the program runs into nirvana.

I will wait for the PIE_ROM then.

KR
Frido

[Ángel Martin]

Sorry about that, here's the missing code in the [INTZL] routine:

Code:

INITLZ    A879    18C    ?FSET 11    
    A87A    3B5    ?C XQ    Stack lift
    A87B    051    ->14ED    [R_SUB]
    A87C    2CC    ?FSET 13    
    A87D    2A9    ?NC XQ    Show "RUNNING" - leaves F8 as-is
    A87E    13C    ->4FAA    [RUNMSG] 
    A87F    1A0    A=B=C=0    zero trinity
    A880    070    N=C ALL    k=0
    A881    2A0    SETDEC    
    A882    001    ?NC GO    iniital sum = 1
    A883    062    ->1800    [ADDONE]

The call to 0x4FAA you can simply remove, it's a cosmetic effect to show "RUNNING..." in the LCD but can be left out without any problems.