Small challenge

[John Keith]

(04-23-2023 02:13 PM)J-F Garnier Wrote:  The implementation of the exponentiation changed between the Capricorn platform (Series 80 / 75C) and the Saturn machines (starting with the 71B).

Does anyone know exactly how the code changed and, more importantly, is either implementation generally better than the other?

The exact value of 23^10 is 41426511213649, and we can see that a 1 ULP error in the least significant digit ( 4142651121365 instead of 41426511213649) can propagate to give the answer 4.14265112137 rather than the correct answer 4.14265112136. However, if 23^10 is computed with 15-digit precision internally, this error should not occur. Since it does occur, it would be nice to know why it occurs in the Saturns and not the others.

[J-F Garnier]

(04-23-2023 06:41 PM)John Keith Wrote:  

(04-23-2023 02:13 PM)J-F Garnier Wrote:  The implementation of the exponentiation changed between the Capricorn platform (Series 80 / 75C) and the Saturn machines (starting with the 71B).

Does anyone know exactly how the code changed and, more importantly, is either implementation generally better than the other?

Here is a comparison between the 75C, the 71B (or any Saturn), and the 35S for a few other known test cases:
- the 729^33.5 vs 3^201 test cited by Pr. Kahan in Mathematics Written In Sand, for 10-digit calculators,
- the 3^729 vs (3^9)^81 test, cited here , source unclear, more suited for 12-digit machines.
- the (1e44)^10.5 worst case for the Saturn (9 ulp error), personal find, see here.

                       75C                71B (Saturn)                35S

729^33.5   7.96841966628e95   7.96841966628e95   7.96841966628e95   correct
3^201       7.96841966627e95   7.96841966626e95   7.96841966626e95   75c slightly better

3^729        6.62818605424e347  6.62818605419e347  6.62818605419e347   75c perfect !
(3^9)^81    6.62818605424e347  6.62818605425e347  6.62818605424e347

(1e44)^10.5  9.99999999998e461 9.99999999991e461 9.99999999991e461   75c clearly better

For the last test, the 35S closely mimics the Saturn, but the 75C is better !

You are right to relate the problem of the exponentiation accuracy to the guard digits, this has been discussed several times (the "2^3=8 New Accuracy" story), see here for instance.

J-F

[Albert Chan]

(04-23-2023 06:41 PM)John Keith Wrote:  The exact value of 23^10 is 41426511213649, and we can see that a 1 ULP error in the least significant digit ( 4142651121365 instead of 41426511213649) can propagate to give the answer 4.14265112137 rather than the correct answer 4.14265112136.

With halfway round-to-even, it will take 2 ULP error to cause incorrect result.

Quote:However, if 23^10 is computed with 15-digit precision internally, this error should not occur.

x^y = 10^(y*log10(x)), 1 ULP exponent error will cause incorrectly rounded 12-digits result.

log10(23) ≈ 1.3617 27836 01759 28788 ...

23^10 ≈ 10 ^ 13.6172783601759 ≈ 4.14265112136463E13
23^10 ≈ 10 ^ 13.6172783601760 ≈ 4.14265112136558E13

Even if we pull out exponent (to reduce errors), we end up the same error.

23^10 = 2.3^10 * 10^10

2.3^10 ≈ 10 ^ 3.61727836017593 ≈ 4.14265112136491E+3
2.3^10 ≈ 10 ^ 3.61727836017594 ≈ 4.14265112136501E+3

[Albert Chan]

(04-23-2023 11:10 AM)robve Wrote:  It is common to find exponentiation by squaring for integer powers, which tends to be more accurate than repeated multiplication and log*exp closed forms.

Unfortunately, this is true only if we don't cause intermediate rounding errors.

(03-22-2022 07:11 PM)Albert Chan Wrote:  When we keep squaring, the errors also blows up exponentially.

x = 1e8 = 0b101 111 101 011 110 000 100 000 000
b = 1+1/x

Free42 (internally only do squarings)
b^(2^08) = 1.000002560003264002763521747927282
b^(2^13) = 1.000081923355125194292436470243333
...
b^(2^26) = 1.956365933428064586618947538663749

All terms to multiply has errors on the *same* side; losing 7 digits precision is normal.

We had this discussion in Free42 possible accuracy flaw thread
Free42/Plus42 were updated to avoid above explosive error issues.

This is result of Plus42 1.0.14 (RHS now all exact)

b^(2^08) = 1.000002560003264002763521747927281
b^(2^13) = 1.000081923355125194292436470243294
b^(2^26) = 1.956365933428064586618947538036231

(03-24-2022 05:01 PM)Thomas Okken Wrote:  

(03-24-2022 04:53 PM)Albert Chan Wrote:  Can I assume new implementation is never worse than internal bid128_pow() ?

Yes: it only uses multiplications when the result of that calculation is exact, and bid128_pow() in all other cases.

[J-F Garnier]

(04-23-2023 06:41 PM)John Keith Wrote:  

(04-23-2023 02:13 PM)J-F Garnier Wrote:  The implementation of the exponentiation changed between the Capricorn platform (Series 80 / 75C) and the Saturn machines (starting with the 71B).

Does anyone know exactly how the code changed and, more importantly, is either implementation generally better than the other?

I traced the operations in my Emu75 and Emu71/DOS emulators.
First thing I immediately realized is that the HP-75C is using 16-digit extended accuracy for internal calculations,
when the Saturn is using only 15 digits, and this explains a lot !

Here are the intermediate results in extended precision on both machines for 3^729, computed as exp(729*ln(3)):

                        75c (16 digits)                71b (15 digits)
ln(3) =               1.098612288668109          1.09861228866810       exact value = 1.098612288668109(69)
729*ln(3) =        800.8883584390514          800.888358439044
exp(729*ln(3)) = 6.628186054237396e347    6.62818605418905e347
rounded to:        6.62818605424e347           6.62818605419e347

Each 71B step is correct, the loss of accuracy comes from the limited guard digits.
This confirms my opinion that 3 guard digits are just too short for the exponentiation on large numbers. It was fine for the 10-digit machines with numbers limited to 9.99..E99, but no more for the Saturn.

So the question is: why did the Saturn team (well, the initial 71B team) use only 15 digits for internal calculations instead of 16, since both the Capricorn and the Saturn CPUs have 64-bit, i.e. 16-nibble registers?
I don't have a clear answer, but managing only 15 digits simplifies the code (adding two 15-digit numbers naturally fits in a 16-digit register) and makes it faster, an important aspect for the quite slow HP-71B. Furthermore, this is how the 41C was doing: internal 13-digit numbers in 14-digit registers.
So it was probably a deliberate choice, a trade-off between efficiency and accuracy.

J-F

[brouhaha]

(04-23-2023 02:13 PM)J-F Garnier Wrote:  The implementation of the exponentiation changed between the Capricorn platform (Series 80 / 75C) and the Saturn machines (starting with the 71B).

I don't have direct knowledge of how this stuff evolved inside HP, but from looking at the code, the Saturn math routines appear to be more directly evolved from the Nut math routines. The Capricorn routines may also have done that, but it looks to me like they were a separate development line, rather than being an intermediate step between the Nut and Saturn routines.

Prior to Saturn, the documentation for the math routines was separate from the source code, and apparently included flow charts. I suspect that even for Saturn, there was some additional documentation on the math routines outside the source code. As far as I know, none of HP's internal math routine documentation outside of source code has ever been publicly seen.

[John Keith]

(04-24-2023 01:35 PM)J-F Garnier Wrote:  I traced the operations in my Emu75 and Emu71/DOS emulators.
First thing I immediately realized is that the HP-75C is using 16-digit extended accuracy for internal calculations,
when the Saturn is using only 15 digits, and this explains a lot !
...

Very informative, thanks! It's funny how thirty-odd years ago, 12 digits seemed like a luxury. Now, in a world with Free42 and NewRPL, a mere 12 digits seems rather dated.

[Christoph Giesselink]

(04-24-2023 01:35 PM)J-F Garnier Wrote:  I traced the operations in my Emu75 and Emu71/DOS emulators.
First thing I immediately realized is that the HP-75C is using 16-digit extended accuracy for internal calculations,
when the Saturn is using only 15 digits, and this explains a lot !

I had a look at the HP48G-R and HP48S-E sources.

The assembly entry for computing

** Name:(S) YX2-15 - Y to the X power

is =cYX2-15 (#773D7 48SX).

It's a litte bit tricky to reach this breakpoint because it's in the covered area.

I found several references for 16 digit mantissa handling

Code:

*   16 digit result of x*ln(y)   *
 YX100  A=A-1  A              DEC EXP
        C=D    W
        BCEX   W              USE 16 DIGIT MANTISSA
        GOTO   EXP16M         Return E^(X*LN(Y))

but the used

=MULTF (#2B91E 48SX) multiplication subroutine

Code:

** Name:(S) MULTF   -  Multiply for finite args only
...
**      Reg. D has the 16 digit mant. of y*x if D(S)#0,
**        (mant of Inf & NaN is not put into D,but D(S)=0 here)
**      Results are truncated to 15 digits.

truncate the result to 15 digits.

[J-F Garnier]

(04-25-2023 07:13 PM)Christoph Giesselink Wrote:  I had a look at the HP48G-R and HP48S-E sources.

I found several references for 16 digit mantissa handling

You're right Christoph, in certain cases, the Saturn does use 16-digit BCD numbers.

I know these two cases:
-when subtracting two 15-digit BCD numbers, the larger mantissa is shifted left (moved to 16 digits) (1)
- the multiplication of two 15-digit mantissa can provide a mantissa on 15 or 16 digits, for instance 3.000.. x 3.000.. --> 9.000.. (on 15 digits), but 3.000.. x 4.000.. --> 12.000.. (on 16 digits). The result is always normalized on 15 digits but the possible 16-digit mantissa is kept (in D) and sometimes used to improved the accuracy of internal calculations, and this is the case for y^x computed as exp(x*ln(y)).

But these are exceptions, or tricks to improve the accuracy in a few special cases.
The Saturn generally uses 15-digit numbers and the Capricorn always uses 16-digit numbers for internal calculations.

A direct proof (71B, 75C: using Math ROMs):
define the vectors A= [1/3 , 1/21 ] and B = [1/7 , -1]
do DOT(A,B), that will compute 1/3*1/7 (w/ internal extended precision) - 1/21 (12-digit input value)
the Saturn returns -4.77e-14
and the Capricorn returns -4.762e-14 , i.e. 1 more significant digit.

J-F

(1) Excerpt from the subtraction (actual addition) code in math0.a:

Quote:* If signs differ, shift mantissa of larger magnitude left *
* to achieve additional digit of precision in subtract. *
* (Same signs could cause overflow out of 16 digit word) *

[Guenter Schink]

(04-22-2023 07:24 PM)Thomas Klemm Wrote:  The following produces the exact same sequence on an HP-48G:

'n^5^2'
'n'
5
24
1
SEQ

Also off topic is this code snipped for Python on Prime:
for i in (x**10 for x in range(5,25)):print(i)

Quote:9765625
60466176
282475249
1073741824
3486784401
10000000000
25937424601
61917364224
137858491849
289254654976
576650390625
1099511627776
2015993900449
3570467226624
6131066257801
10240000000000
16679880978201
26559922791424
41426511213649
63403380965376

Günter

[John Keith]

(04-25-2023 09:56 PM)Guenter Schink Wrote:  Also off topic is this code snipped for Python on Prime:
for i in (x**10 for x in range(5,25)):print(i)

Quote:9765625
60466176
282475249
1073741824
3486784401
10000000000
25937424601
61917364224
137858491849
289254654976
576650390625
1099511627776
2015993900449
3570467226624
6131066257801
10240000000000
16679880978201
26559922791424
41426511213649
63403380965376

Günter

Can also be done in the Prime CAS or in Exact mode on the 49/50, but no suspense regarding significant digits.

[J-F Garnier]

Some more clarifications:

(04-25-2023 08:49 PM)J-F Garnier Wrote:  -when subtracting two 15-digit BCD numbers, the larger mantissa is shifted left (moved to 16 digits)

I would better say: subtraction is done on 16 digits, but in any case the result is normalized to 15 digits.
An example of the benefit of the 16-digit subtraction:

Do 1 - 1/3000 in extended precision, that is:
  1.00000000000000
- 0.000333333333333333
Computing with 15 digits only, you would get, after alignement of the mantissa:
  1.00000000000000
- 0.00033333333333
= 0.99966666666667 i.e. only 14 significant digits

But doing the operation on 16 digits (by shifting both mantissa left 1 position):
  1.000000000000000
- 0.000333333333333
= 0.999666666666667 i.e. now 15 significant digits

Quote:-the multiplication of two 15-digit mantissa can provide a mantissa on 15 or 16 digits, for instance 3.000.. x 3.000.. --> 9.000.. (on 15 digits), but 3.000.. x 4.000.. --> 12.000.. (on 16 digits). The result is always normalized on 15 digits but the possible 16-digit mantissa is kept (in D) and sometimes used to improved the accuracy of internal calculations, and this is the case for y^x computed as exp(x*ln(y)).

Note that in the 3^729 case, the term x*ln(y) = 729 x 1.098... doesn't generate a 16-digit mantissa, no extra accuracy here.

Same for the 1e44^10.5 case, the term x*ln(y) is 10.5 x 101.31..
Here are the intermediate results in extended precision on both machines for 1e44^10.5:
                              75c (16 digits)                71b (15 digits)
ln(1e44) =                101.3137440917379        101.313744091738    exact value = 101.3137440917380(10)
10.5*ln(1e44) =        1063.794312963247        1063.79431296324
exp(10.5*ln(1e44)) = 9.999999999977472e461   9.99999999990893e461
rounded to:              9.99999999998e461          9.99999999991e461

The difference comes from the result of the step 10.5*ln(1e44) and the missing 16th digit on the 71b.

J-F

[Dave Britten]

Some interesting observations:

(04-22-2023 07:24 PM)Thomas Klemm Wrote:  The following produces the exact same sequence on an HP-48G:

'n^5^2'

...

A lot of calculators apparently get nested powers wrong and treat this as (n^5)^2 instead of the correct n^(5^2). I tried this on my Sharp PC-1475 in double-precision mode and couldn't figure out at first why I was getting such larger results! Turns out it's because the Sharp interprets it correctly.

Also, the TI-84 maintains the full 14-digit precision (and gets the correct result of ...5376 on the last term), but only displays 10-digit mantissas on the home screen. But if you use the seq() function to generate and store a list, then pull it up in the stat editor, you can scroll to the bottom elements and hit enter to see them displayed with the full precision.

[J-F Garnier]

To go back to the source of the challenge, the 9825/9835/9845 actually are using a square and multiply method to calculate exponentiation with integer powers.
See the comments in the exponentiation code in this US patent for the 9825 (page 518 of the pdf).
Integer power must be understood in the context of these machines as -32768 <= power < 32768

Let's test the example cited in the code comments:
                        HP-9845                    Capricorn / Saturn
.999^(-32768) = 1.73014224633e14      1.73014224721e14

The Capricorn/Saturn result is exact to 12 digits, but the accuracy of the 9845 doesn't exceed 10 digits (within +/1 ulp) for this test.

We can reproduce the 9845's result with this command line on the 75C/71B:
>P=.999 @ FOR I=1 TO 15 @ P=P*P @ NEXT I @ X=1/P @ DISP X
1.73014224633E14

This demonstrates that the 9835/9845 didn't use internal extended precision for the exponentiation.
Yet the square and multiply method was more accurate than the y^x=exp(x*ln(y)) formula,
as demonstrated with 1/(.999^32768) = 1.73014204403e14 computed on the 9845.

So no, the 9835 was not more accurate than the Saturn (most of the time) !

It seems the HP-85 (1979) was the first HP desktop to use extended precision for internal calculations, probably using the experience of the accuracy improvement for handheld calculators (HP-22, 1976).

J-F
(9845 tests conducted on the MAME 9845 emulator).

Edit: rephrase "multiplication method" by "square and multiply method."

[Gerson W. Barbosa]

(04-22-2023 05:29 PM)Gerson W. Barbosa Wrote:  [

I’m away from home and have only my HP 50g at hand otherwise I would try the HP-12C Platinum (a 12-digit calculator actually) and the HP-33s.

"Home, home again, I like to be here when I can"

HP 12C Platinum:

23 ENTER 10 y^x EEX 13 ÷ 10 × g FRAC f PREFIX -> 4265112136

hp 33s:

23 ENTER 10 y^x right-shift SHOW -> 414265112136

[J-F Garnier]

(04-28-2023 12:52 AM)Gerson W. Barbosa Wrote:  HP 12C Platinum:
23 ENTER 10 y^x EEX 13 ÷ 10 × g FRAC f PREFIX -> 4265112136

hp 33s:
23 ENTER 10 y^x right-shift SHOW -> 414265112136

Thanks!

So the 12C Platinum, the 33S and the 35S are in the same class.
Maybe are they using a square and multiply method for integer powers. Does anybody know?

J-F

Edit: rephrase "multiplication method" by "square and multiply method."

[EdS2]

(04-27-2023 07:31 PM)J-F Garnier Wrote:  To go back to the source of the challenge, the 9825/9835/9845 actually are using a multiplication method to calculate exponentiation with integer powers.
See the comments in the exponentiation code in this US patent for the 9825 (page 518 of the pdf).

Wow - thanks for that, what a document!

[robve]

(04-24-2023 01:07 PM)Albert Chan Wrote:  

(04-23-2023 11:10 AM)robve Wrote:  It is common to find exponentiation by squaring for integer powers, which tends to be more accurate than repeated multiplication and log*exp closed forms.

Unfortunately, this is true only if we don't cause intermediate rounding errors.

Yes and no. I should have provided more context.

You're right that simulating exponentiation by squaring in BASIC, C or Lua or other PL may produce rounding errors (and less accurate results compared to exp-log) if guard digits are rounded away.

However, calculators since the mid 70s use one or more guard digits. It is perfectly safe to use this method with small integer powers \( n \) to speed up \( x^n \) since we can perform multiplications with guard digits. I would not be surprised if some HP calculators use exponentiation by squaring if there is enough ROM space to add the necessary logic. But I'm not an expert on the history of HP math internals and so I can't tell.

As for exp-log in IEEE-754, library implementations are free to use extra guard digits. Or use higher precision floats, see e.g. Elementary Functions - Algorithms and Implementation. We can do that too, when computing \( x^n \) with squaring.

With IEEE-754 representations without guard digits, the method is indeed less effective. It depends on the mantissa value of \( x \) if we can avoid rounding completely, given a small integer power \( n \).

If all else fails, we can place an upper bound on the final error to bound integer powers within an acceptable range.

- Rob

[J-F Garnier]

(04-28-2023 07:13 AM)J-F Garnier Wrote:  

(04-28-2023 12:52 AM)Gerson W. Barbosa Wrote:  HP 12C Platinum:
23 ENTER 10 y^x EEX 13 ÷ 10 × g FRAC f PREFIX -> 4265112136

hp 33s:
23 ENTER 10 y^x right-shift SHOW -> 414265112136

So the 12C Platinum, the 33S and the 35S are in the same class.
Maybe are they using a square and multiply method for integer powers. Does anybody know?

I decided to have a closer look at this question. Below are my results.

Summary

The 35S, and presumably the 12CP and 33S too, are using a square and multiply method for exponentiation with positive integer powers up to about 127.

Method

I measured the duration of the 2^x calculation for different values of the power x on both the HP-32S (Saturn) and HP-35S.
Here is my test program for the 32S/35S:

Code:

A01 LBL A
    STO B
    2
    STO A
    100
    STO C
    RCL A
    ENTER
    ENTER
    ENTER
A11 RCL B
B01 LBL B
    y^x
    CLx
    LASTx
    y^x
    CLx
    LASTx
    y^x
    CLx
    LASTx
    y^x
    CLx
    LASTx
    y^x
    CLx
    LASTx
    DSE C
B18 GTO B

On the 32S, the execution time is almost independent from the power value, and is about 22s for the 500 y^x evaluations.

On the 35S, there is a big difference between a positive integer power and a fractional/negative power, and the execution time increases with larger integer powers, more precisely with the number of bits '1' in the binary representation.

Results for 2^x on the 35S

 x  time (500x)
  4  17s     positive power
 -4  38s     negative power
  4.5 31s     non-integer power
  7  19s     3 bits in binary representation of the power
  8  19s
 15  23s     4 bits ...
 16  21s
 31  27s     5 bits ...
 32  23s
 63  31s     6 bits ...
 64  25s
100  29s
120  32s
127  34s     7 bits ...
128  35s
>128 35s

The execution time is then quite constant for powers beyond 128.

These results strongly suggest an exponentiation by squaring for positive powers up to 127.

Benefit ?

The question is: what is the benefit of implementing the exponentiation in this way for these limited special cases, between 2 and 127?
As the initial challenge demonstrated, there are a few cases where the exponentiation by squaring provides better results, but this is globally rare.
Moreover, and more importantly, better results are not granted for the exponentiation by squaring.

Here is a counter example:
.96279338^63
32S:    9.17455514137e-2
35S:    9.17455514136e-2
exact:  9.174555141365119..e-2

For this case, the 32S provides the correct rounded value, but not the 35S.

So definitively the accuracy is not significantly improved, the only small benefit is a faster result for small positive integer powers.

More historic notes

I showed above that the HP 9825/35/45 series machines (1976-78) were using exponentiation by squaring for integer powers between -32768 and 32767.
The reason is quite clear, it is due to the modest accuracy of the log/exp function of these machines, without guard digits as were the early HP handhelds too.
Example on the 9845:
2^3 = 8 exact (exponentiation by squaring)
but 4^1.5 = 7.99999999988

However, the exp/log implementation became soon accurate enough on the later Series 80 (1980) and Series 200 BASIC (1981), so the exponentiation by squaring was abandoned and the exponentiation y^x solely relied on the formula exp(x*log(y)) for all power values even for x=2.
The accuracy was similar to the exponentiation by squaring method, in the range of +/1 ULP for not-too-large results, with a few worst cases for larger numbers.

Interesting enough, the HP BASIC manuals for the Series 200 (for instance here, dated 1982) included a section "Exponentiation vs. Multiply and SQR" in the chapter "Efficient Use of the Computer's Resources":

Quote:Exponentiation is very slow when compared to other mathematical operations. [...]
The most common examples are squaring a number or finding a square root. [...]
... dramatic savings can be gained when raising numbers to an integer power. [...] [by using] repeated multiplication.

J-F

[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