Exponential operator (and the numerical error in using it)

Exponentiation to an integral power, or

r = x ** n

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.

Hi Damian —

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?)

-Brad

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.

Thanks again.