Topic: X to the Y power (x@3 != x*x*x, for math lovers only)


Author: ark@alice.att.com (Andrew Koenig)
Date: 29 Feb 92 16:55:33 GMT
Raw View
In article <22305@alice.att.com>, I asked:

>  Is x@3 guaranteed equal to x*x*x?  What if x*x*x is not
>  the best possible approximation to the infinite-precision
>  cube of x?  This problem becomes worse for higher powers.

Most people will not find it immediately obvious that x@3 should be
anything but x*x*x, so here's an example.  In this example I will follow
Joe Buck's suggestion and use @ for exponentiation, mostly to see how
it looks in print.

Assume we're using IEEE single-precision floating-point.  It is possible
to construct an example for double-precision as well, but the numbers
would be bigger, hence harder to follow.

The relevant characteristics of IEEE single-precision are that the fraction
is 24 bits long and longer results are rounded to the nearest available
value, with rounding to even in case of an exact tie.  Thus, if you
have to round this way:

 ...001101|011000101...

(where the | represents the boundary between the low-order bits of the
actual result and the remaining bits of the infinite-precision result)
then you leave it alone and give ...001101 as the final result; if
instead you have to round this way:

 ...001101|101100...

then you must round up to ...001110, and finally if you have

 ...001101|10000000....

(with all the trailing bits being 0), you must still round up to ...001110
because it's a tie and you must therefore round to even.  Similarly,

 ...001100|10000000....

would be rounded to ...001100 because of the round-to-even rule.

Given that background, let's look at 4097.0 @ 3:

First, 4097 is 2@12 + 1.  Recall that (x+y)@3 is x@3 + 3*x@2*y +
3*x*y@2 + y@3, so 4097@3 is 2@36 + 3*2@24 + 3*2@12 + 1.  In binary,
that is:

 1000000000011000000000011000000000001

We must round that to 24 bits, however, which means we get this:

 100000000001100000000001|1000000000001

so the final result according to IEEE rounding is

 100000000001100000000010 * 2@13

or 2@36 + 3*2@24 + 2@14.

Now let's calculate x*x*x.  That would presumably be evaluated as

 x1 = x * x;
 result = x1 * x;

where x1 is also a single-precision floating-point number.  In other
words, the result of x*x must be rounded before multiplying again by x.

OK, here we go.  Recall that (a+b)@2 is a@2 + 2*a*b + b@2, so
x@2 is 2@24 + 2*2@12 + 1, or

 100000000001000000000000|1

which, by the round-to-even rule becomes 2@24 + 2*2@12, or 2@24 + 2@13.
Now we must multiply that by 2@12 + 1, and we get 2@36 + 2@25 + 2@24 + 2@13,
or 2@36 + 3*2@24 + 2@13.

You will notice that that is *not* the same as our correctly rounded
result for x@3, which was 2@36 + 3*2@24 + 2@14.

So that is the point of the question: as a user, should I be entitled to
assume that x@3 is the same as x*x*x?  Even if x*x*x is not the best
possible approximation to `x cubed' ?
--
    --Andrew Koenig
      ark@europa.att.com




Author: cimshop!davidm@uunet.UU.NET (David S. Masterson)
Date: 1 Mar 92 22:54:41 GMT
Raw View
>>>>> On 29 Feb 92 16:55:33 GMT, ark@alice.att.com (Andrew Koenig) said:

[...much on binary-level calculations deleted...]

> So that is the point of the question: as a user, should I be entitled to
> assume that x@3 is the same as x*x*x?  Even if x*x*x is not the best
> possible approximation to `x cubed' ?

I'm sorry, but the value of this question eludes me.  What is being discussed
is the value of adding an operator to the language, not the precision of the
calculations of that operator.  The above question can be shown to be moot by
simply asking:

 Can I assume that pow(x,3) is the same as x*x*x in ANSI-C?

Other standards have avoided this question for a long time and gotten
reasonable results.  You, as a user of a standard conforming compiler, are
only allowed to assume that pow(x,3) (or x@3) will give a reasonable answer to
the question 'what is x raised to the power of 3?'  The precision of the
calculations will get continuously better as the hardware and compiler develop
and, so, even if we attempted to provide a definition for 'reasonable', it
would most likely change over time.
--
====================================================================
David Masterson     Consilium, Inc.
(415) 691-6311     640 Clyde Ct.
uunet!cimshop!davidm    Mtn. View, CA  94043
====================================================================




Author: pgh@stl.stc.co.uk (Peter Hamer)
Date: 3 Mar 92 16:25:39 GMT
Raw View
In the referenced article ark@alice.att.com (Andrew Koenig) writes: ... a great
deal of painful truth about numerical computation.

But to emphasise something I only implied in another posting,

1) The C++ language committee should try to provide the syntactic hooks to
   enable numerical people to write C++ expressions using infix notation for
   operators which are traditionally infix.

2) The numerical community should be left to develop libraries which support
   solutions to numerical problems: it's their problem; they have the expertise;
   extension by library is part of the C++ culture.

The point that should be at issue is: how can hook(s) be provided to enable
a significant potential user community to develop something(s) that they
would really like.

A general solution may not be possible, but a special one for the power operator
doesn't seem technically unrealistic. [Rejection because the compiler producers
decided that they didn't think it worth the effort would be another matter.]

Regards, Peter