Π day

[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