Efficient floating point exponentiation is allowed AI-00868/00 1
90-02-01 BI WI
!standard 04.05.06 (06) 90-02-01 AI-00868/00
!class binding interpretation 90-02-01 (provisional classification)
!status work-item 90-02-01
!status received 90-01-29
!references AI-00137
!topic Efficient floating point exponentiation is allowed
!summary 90-02-01 (DRAFT)
Efficient methods for performing floating point exponentiation can be used
even if the result does not lie within the correct model interval.
!question 90-02-01 (DRAFT)
AI-00137 emphasizes that it is not proper to use factored approaches to
computing the result of exponentiations. For example, A**4 being computed as
(A**2)*(A**2) is one simple case. A more general case is the common approach
of computing successive powers of 2, which allows A**B to be computed using
O(log(B)) multiplications instead of O(B) multiplications.
A major part of the reasoning in coming to this decision was the pragmatic
observation that it was very unusual to have large powers in real
exponentiations. The discussion section in the commentary said:
In practice, performing N-1 multiplications is not inefficient,
because 90% of the time the exponent is 2, 7% of the time, it is
3, and only in 3% of the cases is the exponent larger than three
[Wichmann, B. A., "Algol 60 Compilation and Assessment," Academic
Press, 1973]. So the possibility of reducing the number of
multiplications does not arise very often.
This observation is significantly false in two situations (and they may well
be others).
1) In writing conversion routines in Ada, large powers of 10 or
2 are often needed, leading to expressions such as 2.0**N
and 10.0**N, where N can be large (say of the order of 100).
2) The computation of compound interest rates where interest
compounds daily requires exponentation by large powers. For
example, to compute the annual rate of return for an account
paying 0.03% interest daily, the expression 1.0003**365 may
be computed.
Computing these cases by successive multiplication is not only appallingly
slow, but is also much less accurate than using the logarithmic method, and
indeed for expressions of the form 1.0003**365, it is probably more
appropriate to treat the situation as a call on a math library containing a
specialized exponentiation routine. As with trigonometric functions,
requiring model number accuracy of such routines is certainly not always
reasonable. Note that with the current AI, 1.0003**365 is only model number
accurate in the sense of being within the (very wide) model interval of the
result of 365 successive multiplications. If a math library does contain such
DRAFT DRAFT DRAFT
Efficient floating point exponentiation is allowed AI-00868/00 2
90-02-01 BI WI
an exponentiation routine, then it will in general be more accurate than
doing the successive multiplications. However, it is not clear whether it can
easily meet the requirement that its result ALWAYS be in the model interval
implied by successive multiplication. The analysis required to assure this is
complex, and we simply don't know enough to understand the consequences of
such a requirement.
As an example of our lack of understanding, consider the case of A**4 being
computed as (A**2)**2. Brian Wichman has done an analysis showing that this
substitution is valid for rounded 32-bit IEEE floating point arithmetic,
i.e., there are no cases in which the result of the streamlined computation
is outside the model number interval implied by the sequence of three
multiplications.
However, Brian has not succeeded in carrying out the analysis for the case of
64-bit IEEE floating point. An implementor is in the somewhat strange
position of generating inefficient code JUST IN CASE there is an example
where 64-bit IEEE floating-point gives the "wrong" result. Shouldn't the
accuracy requirements for exponentiation be determined based on what is
feasible using efficient methods for computing the result?
!recommendation 90-02-01 (DRAFT)
[State an accuracy requirement after consultation with the NRG.]
!discussion 90-02-01 (DRAFT)
[To be provided after discussion by the ARG.]
!appendix 90-01-29
*****************************************************************************
!section 04.05.06 (06) Software Leverage, Inc. 83-10-29 83-00192
!version 1983
!topic Exponentiation with floating point operand
The RM demands an inefficient implementation of exponentiation for
floating point operands.
The problem is that the RM defines the model interval of the result in
terms of the model intervals of a series of multiplications from left
to right. Unfortunately this sometimes excludes model numbers
obtained for calculations used in usual algorithms, e.g. the one that
uses the identity
x**(2*n+e) = (x*x)**n * x**e
where e = 0 or 1, so x**e is either 1.0 or x.
The simplest case of this is calculating x**4 as sqr(sqr(x)) where
sqr(x) is shorthand for x*x. Even assuming we have a machine that
rounds, there are cases in which this method will fail. Adding a
DRAFT DRAFT DRAFT
Efficient floating point exponentiation is allowed AI-00868/00 3
90-02-01 BI WI
guard bit doesn't help; the only result we have been able to obtain
is that with a number of guard bits equal to the number of bits in the
exponent, one can construct an algorithm using the above which works.
There is no mathematical result we can find to improve this.
It was presumably not the intent to require that, for example, x**128
take 127 multiplications to calculate, but this follows from 4.5.6(6).
Deviation from the manual is difficult to test for here, but imple-
menters should be able to use power trees a la Knuth without fear. In
this case, we believe the language in the RM should be changed; there
remains the question of what to do until the next Standard. Since
conformance is very hard to enforce anyway, it presumably would suf-
fice to publish some form of statement to the effect that deviation in
this area is allowed if the implementation instead conforms to a
changed version of 4.5.6(6), perhaps along the following lines:
"Exponentiation with a positive exponent delivers the left operand for
an exponent of one, and otherwise (for exponent n) is equivalent to
multiplication of the results of exponentation of the same left
operand to powers k and n-k, for some positive k less than n. For an
operand of a floating point type..." (only the first sentence
changes).
One might wish to also change 4.5.7(9), for the sake of clarity, to be
something like this:
For the result of exponentiation, the model interval associated with
the expression x**n is determined as follows:
1. If n = 1, the model interval is that of the expression x.
2. If n > 1, the model interval is the union of the model inter-
vals for expressions (x**k)*(x**(n-k)) for all k in the range
1 <= k < n.
3. If n = 0, the model interval consists of the single number
1.0.
4. If n < 0, the model interval is that for the expression
1.0/(x**(-n)).
This might need changing for the sake of rigor (evaluation of x might
cause side effects); but the idea is there.
Counterexamples where sqr(sqr(x)) doesn't give the answer required by
the RM for x**4 are given below; the list is complete for mantissas
up to 17 long. They are expressed in octal notation with "b" as the
mantissa length (some values of b don't come up in Ada, but the exam-
ples suggest that counterexamples can be found for arbitrarily long
mantissas).
The situation would presumably be worse if machine rounding wasn't
assumed.
DRAFT DRAFT DRAFT
Efficient floating point exponentiation is allowed AI-00868/00 4
90-02-01 BI WI
b = 14:
0.43642
0.44650
0.64452
0.64032
b = 15:
0.45360
0.62234
0.63025
0.63116
b = 16:
0.442350
0.445264
0.453420
0.454100
0.611304
0.641400
0.642074
b = 17:
0.427562
0.443354
0.447526
0.447626
0.451304
0.455276
0.455472
0.455626
0.457266
0.624100
0.624720
0.631352
0.643616
0.647376
0.652222
These examples also show that a small number of guard bits doesn't
suffice. The number 0.62234 fails for b = 14 when 1 guard bit is used
(so the machine mantissa length is 15).
*****************************************************************************
!section 04.05.06 (06) Jean D. Ichbiah 84-03-01 83-00312
!version 1983
!topic Exponentiation by real exponent
There is nothing wrong with the definition given by the RM.
Nothing prevents one from having a function EXP that uses repeated
squaring.
(Most exponentiations are by small integers - so why bother with
pathologies that would make the language more complex. The problem was
DRAFT DRAFT DRAFT
Efficient floating point exponentiation is allowed AI-00868/00 5
90-02-01 BI WI
known; the design choice deliberate).
*****************************************************************************
!section 04.05.06 (06) R. Dewar 90-01-29
!version 1983
!topic Allow efficient exponentiation methods
A previous AI emphasizes the fact that it is not proper to use factored
approaches to computing the result of exponentiations. The example of
A**4 being computed as (A**2)*(A**2) is simple one case. Another more
general case is the general logarithmic approach of computing successive
powers of 2, which allows A**B to be computed using O(log(B)) multiplications
instead of O(B) multiplications.
A major part of the reasoning in the ARG coming to this decision was (as I
remember) the pragmatic observation that it was very unusual to have large
powers in real exponentiations.
This observation is significantly false in two situations (and they may
well be others).
1. In writing conversion routines in Ada, large powers of 10 or 2 are
often needed, leading to expressions such as 2.0**N and 10.0**N,
where N can be large (say of the order of 100).
2. The computation of compound interest rates where interest compounds
daily requires exponentation by large powers. For example, to compute
the annual rate of return for an account paying 0.03% interest daily,
the expression 1.0003**365 may be computed.
COmputing these cases by successive multiplication is not only appallingly
slow, but is also much less accurate than using the logarithmic method, and
indeed for expressions of the form 1.0003**365, it is probably more appropriate
to treat the situation as a call on a math library containing a specialized
exponentiation routine. As with trigonometric functions, requiring model
number accuracy of such routines is certainly not always reasonable. Note
that with the current AI, 1.0003**365 is only model number accurate in the
sense of being within the (very wide) model interval of the result of
365 successive multiplications. If a math library does contain such an
exponentiation routine, then it will in general be more accurate than
doing the successive multiplications. However, it is not clear whether it
can easily meet the requirement that its result must ALWAYS be in the model
interval implied by successive multiplication. The analysis required to assure
this is complex, and we simply don't know enough to understand the
consequences of such a requirement.
As an example of our lack of understanding, consider the case of A**4 being
computed as (A**2)**2. Brian Wichman has done an analsys that shows that
this substitution IS valid for rounded 32-bit IEEE floating point arithmetic,
i.e. there are no cases which arise in which the result of the streamlined
computation is outside the model number interval implied by the sequenceof
DRAFT DRAFT DRAFT
Efficient floating point exponentiation is allowed AI-00868/00 6
90-02-01 BI WI
three multiplications.
However, Brian has not succeeded in carrying out the analysis for the case
of 64-bit IEEE floating point. An implementor is in the somewhat strange
position (if he is being VERY honest) of generating inefficient code JUST
IN CASE there is an example where 64-bit IEEE floating-point gives the
"wrong" result. This seems plain silly.
Recommendation
The AI should be withdrawn. The NRG should be asked to recommend appropriate
error bounds for the result of computing the result of an exponentiation.
Relation to Numeric Standards Development
The NRG is currently developing a package of basic numeric functions. This
package does NOT currently include an exponentiation operation, on the
grounds that this operation is included in the language. However, when I
it was pointed out that the semantics of exponentiation required repeated
multiplication (which is in the general case unacceptably inaccurate), there
was discussion of whether to put such an operation into the package. The
general response was that it would be better to fix the Ada operation.
DRAFT DRAFT DRAFT