Small challenge

[robve]

(05-16-2023 06:57 PM)J-F Garnier Wrote:  I decided to have a closer look at this question. Below are my results.

Very nice work!

There is not really much to add to your in-depth assessment of the inner-workings of these HP machines. Quite interesting stuff. Nice to see these trade-offs in speed and accuracy.

Rounding error

With exponentiation by squaring, we need to perform up to \( k \) multiplications, where \( k \le 2\lceil\log_2 n\rceil \). The rounding error in the result is about \( k\varepsilon \), since \( (1+\varepsilon)^k \approx 1+k\varepsilon \) for small MachEps \( \varepsilon \ll 1 \). For floating point, \( \varepsilon \) is 0.5 ULP.

When limiting \( n\le 127 \) we have \( k=14 \) so the error is at most \( 14\varepsilon \) or only 7 ULP. However, this all depends on the accuracy of BCD multiplication, which may or may not be accurate to the last guard digit or +/- 0.5 ULP.

Another consideration is when \( x \) has only a few leading digits in the mantissa, i.e. squaring is exact. However, each multiplication doubles the number of leading digits in the mantissa, so subsequent steps in exponentiation by squaring are likely affected by round-off issues eventually. Interesting special cases arise when \( x \) is a power of 10 for BCD machines and a power of 2 for binary floating point.

Benefit?

Computing \( x^n \) with exp-log based \( x^n = \exp(n\log x) \) is always approximate. But we don't need to be concerned about the approximation error when the exp and log approximation algorithms (e.g. power series, CORDIC, polynomials, and so on) use a higher internal accuracy with guard digits to ensure that e.g. \( 3^3 = 27 \) is exact after rounding off. With exponentiation by squaring, \( 3^3 = 27 \) is guaranteed to be exact. Remember the warnings that came with early BASIC pocket computer manuals that stated that expressions like X^2 <= 9 may fail when X=3, because 3^2 is 9.000000001 (or something along those lines.)

Perhaps it should be said that if guard digits aren't important to keep as exact as possible in chain calculations, or if the exponentiation-by-squaring algorithm uses a higher accuracy internally (e.g. one or more extra guard digits that can be tossed away to produce an accurate result), then exponentiation by squaring makes sense, because exponentiation by squaring is very fast (faster than exp-log approximations in software) and at least as accurate as exp-log approximations or more accurate, e.g. when \( x \) is also a small integer. The latter is what I alluded to before, that exponentiation by squaring can be more accurate (but only when these conditions are met.)

I've seen the range -32768 .. 32767 is in the assembly listings. That range of \( n \) in exponentiation by squaring of \( x^n \) is rather large. I did some quick tests and the results aren't great. Since that range is also used in Forth850, I'm making some changes to the Forth850 math routines to maintain accuracy with exponentiation by squaring limited to \( |n|\le 16 \) for IEEE single precision floating point. Since exp and log are implemented as power series in Forth, calculating \( x^n \) with exp-log is considerably slower than exponentiation by squaring when \( |n|\le 16 \).

- Rob