Will have errors. Assume we define ULP as a unit in the last place, a bit like the machine epsilon.
The computation could be implemented as
r = x * x; // if n == 2
t = x * x, r = x * t; // if n == 3
t = x * x, r = t * t; // if n == 4
t = x * x; r = t * (t * x); // if n == 5
t = x * x * x, r = t * t; // if n == 6
The error the each of the above will never exceed (n / 2) * ULP.
Beyond that, the errors just grow depending on how you do the calculations.
Instead of this, .....
If x is one of e (Napiers constant) or 10, I could use exp or exp10. In this case, I know the error I get would never exceed 1*ULP.
If x is 2, I would just add to n the IEEE 754 exponent bias of the type of x, a function of the number of bits in that type, and inject it into a floating point representation. The error here would be 0*ULP, i.e. zilch
I could also use the new version of the pow routine which accepts an integral exponent and know that the error never exceeded 4*ULP, independent of x
Could you publish how you interpret the exponentiation operator so that users know what error to expect please? Or is it documentated in a place I have not been able to find. I am not really flash at searching HTML pages.
The implementation of the exponentiation operator is not well-documented today, but ultimately ought to be. The lack of documentation also reflects, in part, a lack of confidence that our implementation is as it should be—i.e., we didn't want to specify and lock in how exponentiation works without more confidence that what we are doing is reasonable. And even with that confidence, there's a little bit of a question about how to specify it without closing off opportunities to further optimize it later.
In a nutshell, I believe the implementation effectively is as follows:
complex**complex combinations call cpow() / cpowf(), based on width
real**real combinations call pow()
real**param integral combinations are:
strength-reduced for exponents in the range of 0–8
0: returns 1
1: returns the base
2: returns the base*itself
3: cubes the base via multiplication
...
8: squares the base via multiplication 3 times
call pow() otherwise
int ** [param] int:
uses various heuristics to compute the result... (where I'm guessing you're more focused on the floating point cases?)
As always, I'm interested in your feedback on how to make it better (e.g., what is this new pow routine you refer to?)
Thanks heaps. That's extremely useful. A bit of cut and paste on that and it is ready to hit the prime time Chapel documentation.
One day I will
Do an error review of x7 and x8 against pow and you can decide whether to trim to 0-7;
Ditto cases of -6..0 (of you still like 0-8 or -4..0 if you prefer 0-6) and you can maybe expand to -6..+8
look at ... proc pow(x : real(?w), n : integral) return /* IEEE 754 */ pown(x, n) ... in this context.
The last case is the interface to an (long) recommended IEEE 754 recommended of which I know of any published implementations. We might see one out of a French project but by 2023. But it would certainly be the fallback routine to call in the case of real**param integral and would be the go-to routine for any case of real real and int arguments. With an FMA used inside such a routine, I am guessing that you would probably get close to doubling the effective range, but I am not prepared to guarantee that without testing.
That investigation would even be a good student project or, with some enhancement of scope, a good undergraduate thesis. Probably not at the level of an Master's.
Some libraries use a new table-driven version of the pow() routine. Julia's own openlibm project does not. But versions of MUSL and GLIBC do and have for a while.
I think this information should be in the documentation. It also gets tied up with type promotion rules (which are not explained in context).
That said, given this information, the use of integral powers for floating point data is questionable. Its errors bounds for x8 and x9 will be very different.
And I am not sure I like the error at x8 either, or x7. It is quite scary. Then again, so is any form of exponentiation. It is almost like an "Are you sure you want to do this" situation with this operator. I would be tempted to assume and ensure the existence of
pow(x : complex(?w), y : complex(w))
pow(x : real(?w), y: real(w))
pow(x : real(?w), n : int(8))
pow(x : int(?w), n : int(8))
with param n versions for the last two, each with their own known and published errors, and rely on them.
Was there an internal discussion paper on the algorithm to do exponentiation of a floating point number to an integral power. Did people consult any of the references like
Approximating x^n efficiently - Chandrasekhar Narayanaswami et al
Information Processing Letters, Volume 50, Issue 4, 25 May 1994, pp205-210
Computing correctly rounded integer powers in floating-point arithmetic
P Kornerup, C Lauter, V Lefevre, N Louvet, JM Muller
ACM Transactions on Mathematical Software (TOMS) 37 (1), pp1-23
Or some on addition-chain exponentiation. I think Wikipedia's summary is quite adequate. That said, I think both of these approaches are arguably overkill but I could be wrong.
Unlike when the exponential operator was implemented in Chapel, the language now has an FMA (thanks again Michael). We could exploit that instruction and probably come up with a simple and fast algorithm for x**n which has an error like O(log(n)) * ULP((1.0).
Once you have those 4 algorithms, it makes documentation a lot easier. Also, I think the case of complex and imag exponentiation is missing from the current doco.
I went and searched our old Sourceforge mailing list
(Chapel Mailing Lists) but didn't find anything
there. Looking through old internal emails, I see a discussion about
exponentiation from 2009 due to a bug in the implementation because of
converting between types (the bug looked to be due to using the C pow
function but converting the result back to an integer without a cast
when one was necessary for int(32)). It looks like the result of that
conversation was to handle things more manually for integers, but my
hope is that we still just rely on C's implementation for real(32) and
the like. There's also talk about looking for algorithms to support
integer exponentiation in algorithm books, though I suspect the notes on
the actual work done is long since lost to the sands of time.
I'm not sure how relevant any of that is to what you're wondering about,
though.
In answer to your question about where the implementation of
exponentiation is, it looks to be handled by the compiler today.
In general, operators for basic types have their Chapel code
implementation in the internal module ChapelBase
(chapel/modules/internal/ChapelBase.chpl at main · chapel-lang/chapel · GitHub
for the exponentiation operator), and it looks like reals there are
handled by the internal __primitive( ) call, which basically just
means "the compiler handles it". I followed the trail for a little
while but got stymied by its support (and that of many operators)
seeming to be handled via macros - hoping someone more familiar with the
dyno side of that implementation can speak more on that.
Thanks for investigating this Lydia. My fishing results got nowhere.
Anyway, Brad's earlier comments said that simpler algorithms are used for cases of a param integer exponent in the range 0..8. I cannot verify that for the smaller numbers but my testing says otherwise for the higher values in 0..8 unless it is a really smart algorithm. I suspect that that these days it does as the immediately preceding post suggests. But I cannot conclusively prove it which is why I wanted to look at the source.
I think the exponential operator is one of the best (and most expressive and readable operators in existence) and very useful for engineers which often use formula which raises a variable to a small real (non-integral) power. I also think it is one of the worst. Should we even be suggesting to people that they raise a real number to a power as small as say 20. Is that practical (for anything other than where the real number is a power of two - and Chapel already has smarter ways to achieve that). And I grew up using (and teaching) this operator with Fortran. I also grew up using and teaching the Pascal routine pow() which accepts two integral arguments.
I am personally inclined to dismiss the exponential operator, especially in the context of some class notes and a book I am (ever so slowly) putting together. But, my bias is misplaced unless I teach (especially novice programmers) things like Horner's rule and Estrin's formula and fmas and the (as yet missing in Chapel) hypot() routine and so on. See the ASIDE below.
Certainly x**n for n in [-1, 0, +1, +2] with a trivial algorithm satisfies the definition of a basic mathematical operation being accurate to better than ULP(1.0). In fact, it is correctly rounded, i.e. better than ULP(0.5). Outside that range, you either need a simple algorithm, probably based on fma to get the sort of accuracy demanded by an IEEE 754 mathematical operation or the use of the low level pow routine (the latter of which slows things down beyond what most people expect). And that simple algorithm may have problems if the exponent is a negative number. And what do you really do for integers? What are these heuristics that Brad mentioned. I always worry when I see the word heuristic.
At a very low priority level, i.e. after summer breaks, recovery from ChapelCon and the like, it would be nice to find out what the real story is (and enhance the documentation to note the accuracy level or (my suspected) performance hit that one takes if one uses it). I was surprised by the accuracy of exponentiation but I suspect that this comes with a performance hit for using pow().
ASIDE: I am not pushing for the inclusion of hypot() - I just mentioned it as this opens a big can of worms about whether a language like Chapel needs a 3D version as well, faster and less accurate versions, and for us Linear Algebra users, its enhancement to a Givens rotation.
I am going to take stab in the dark and say that if Chapel is passing the calculation for integral exponents off to LLVM, then it is using the powi intrinsic in LLVM these days. I have no idea of its accuracy. Even the LLVM notes say nothing about it. And it has different characteristics when compiling with fast math flags. Somebody needs to verify this as it is outside my skill set.
I might run into the Discourse multiple reply limit shortly.
I will summarize what seems to be behind the Chapel exponential operator which Lydia's exploration seemed to reveal (which is slightly different to what Brad mentioned quite a while ago). But either way, the exponentiation operator is not a conventional arithmetic operator like multiplication, division, or even residual which all have a rounding error which does not exceed 0.5 * ULP(1), or what some people call "correctly rounded". I think the documentation needs to be updated to reflect that fact because it appears in
intermixed with other arithmetic operators all of which do have that single rounding error accuracy. Also, for integral cases, it does not necessarily call a mathematical library routine that would hopefully guarantee the slightly higher rounding error which does not exceed 1.0 * ULP(1). What I think is happening after it has been passed to LLVM's powi intrinsic could have an error of about 6 binary digits for a real(64) in the range [1,2] raised to an integral power (I hope I am able to leverage/quote Moshier's work from the 1980s). I did ask some questions of the LLVM team of this in the last 24 hours and got some useful/helpful replies. If one considers that 6 out of 52 binary digits means several, their reply of several ULPs of error agrees with Moshier's experimentally derived precise data. My own, far less thorough, experiments with non-extreme numbers never gets anywhere near that limit.
I also think the overhead of the exponentiation operator needs to be noted in that same point in the documentation as that aspect is probably just as important for many people as the accuracy issues.
Some words of wisdom like the following need to be added to the
definition of the exponentiation operator (as well as the complex
and imag cases being documented):
This operator is unlike other mathematical operators for floating
point data which are correctly rounded, i.e. have a rounding error
which is not greater than 0.5 * ULP (1). At best, exponentiation of
a floating point number achieves a rounding error which is for a
floating point exponent is guaranteed to not exceed 1 * ULP(1)
with an underlying robust mathematical library. For cases of an
integral exponent of a floating point datum, this will be optimized
for speed (at a sacrifice in accuracy). In this case, the rounding
error can be significantly higher than 1 * ULP(1). For an integral
exponent of n, the error can exceed log(n) * ULP(1). although in
most cases it is only several ULPs.
A more succinct version of the above would be better but I cannot
come up with tighter words for now.
That said, if a Chapel user thinks the above is too technical and
needs a lot of simplification, I would humbly suggest that this user
needs to learn a lot more about rounding errors and their potential
bad impact on the quality of their results.