(EXPT <not-a-double> <complex-double>) uses intermediate single-float values
Error detected by Maxima's test-suite, analysis by Lutz Euler below.
Maxima's calculation of the zeta function with the argument of the
test case is done using the lisp function "float-zeta", defined in
terms each containing a quotient of something by "expt" of a small
integer and the argument to the zeta function, the latter converted
to complex double-float. These "expt"s are calculated directly using
SBCL's "expt". Here is a trace:
Maxima 5.23.1 http://
using Lisp SBCL 126.96.36.199
Type (to-maxima) to restart, ($quit) to quit Maxima.
MAXIMA> (trace expt)
MAXIMA> (float-zeta #$(3+%i)$)
0: (EXPT 4 3.0)
0: EXPT returned 64.0
0: (EXPT 2 #C(-2.0 -1.0))
0: EXPT returned #C(0.1923097243
0: (EXPT #C(7.2421875 1.0) #C(1.25 0.5))
0: EXPT returned #C(4.4185283038
0: (EXPT 1 #C(3.0 1.0))
0: EXPT returned #C(1.0 0.0)
0: (EXPT -1 0)
0: EXPT returned 1
0: (EXPT -1 0)
0: EXPT returned 1
0: (EXPT 21 #C(3.0 1.0))
0: EXPT returned #C(-9217.
0: (EXPT 2 #C(3.0 1.0))
0: EXPT returned #C(6.1539112363
... [many more calls to (expt <small int> #C(3.0 1.0))]
The test case expects the following result:
As a simple reduced test case without maxima we have:
(expt 2 #C(3.0d0 1.0d0))
Clisp (2.48) gives:
and Clisp with more precision (setf (system:
(expt 2 #C(3.0l0 1.0l0))
The definition of "expt" is in sbcl/src/
(defun expt (base power)
(number-dispatch ((base number) (power number))
(((foreach fixnum (or bignum ratio) single-float double-float)
(if (and (zerop base) (plusp (realpart power)))
(* base power)
(exp (* power (log base)))))
The calculation of (log base), base being an integer, yields a
single-float as required by the standard. For the purpose of "expt"
it may be a bug to use (log base) as the standard says in section
188.8.131.52 "Rule of Float and Rational Contagion":
"When rationals and floats are combined by a numerical function,
the rational is first converted to a float of the same format."
For a quick test I replaced the calculation in the SBCL source as follows,
(exp (* power (log (coerce base 'double-float))))
then compiled SBCL and maxima. This SBCL calculates:
* (expt 2 #C(3.0d0 1.0d0))
This seems to differ from Clisp's result only in the last bit of the
real and imaginary part each.
The zeta value calculated by maxima is then:
(%i1) zeta (3 + %i), numer=true;
(%o1) 1.107214408431409 - .1482908671781753 %i
The resulting maxima passes its test suite with no unexpected errors.
In case it matters, the above is under x86-64, except when using
SBCL 0.8.12, which is x86.
If "expt" indeed needs to be modified a correct fix seems a bit more
complicated as currently there are two clauses dealing with complex
exponents, both not distinguishing between single and double floats.
The standard seems to require to look at all argument types and to use
single-float calculations if all arguments are integer or rational or
single-float, be it real or complex, and to use double-float calculations
if at least one argument is double-float real or complex?
Also, the current implementation of "expt" allows (with a complex
exponent) bignums or ratios as the base that are too big or too small
to be represented as a single-float, provided their logarithm can
(and the exponent is such that the result fits in its type).
If one doesn't want to lose this property, a function for log of
ratios and integers returning a double-float would need to be defined.
This may be easy as "log" already contains the needed logic, only
unfortunately always coercing the result to single-float.
On the other hand I believe the standard does not require this, so
hopefully nobody depends on this accidental property. And, when the
exponent is real, bignum/ratio bases too large and too small don't work
On the third hand, as "log" (intentionally, seemingly) deals with bignums
and ratios that don't fit in a float, maybe "expt" should, too, even in
the case of a real exponent?
Moreover, "log" is implemented the way it is to be more exact in the case
that its argument is a ratio very near to one. This carries over to "expt":
The original SBCL calculates:
* (expt (/ (expt 2 64) (1+ (expt 2 64))) #c(10 10))
whereas with the above hack we have the less exact
(Sorry, can't resist: Did you notice that "log" returns zero for rational
arguments very near to one if the integer length of numerator and
denominator differs by one?
(log (/ (expt 2 64) (1- (expt 2 64)))) -> 0.0
(log (/ (1+ (expt 2 64)) (expt 2 64))) -> 5.421011e-20
which result would be correct for the former ratio, too.)
|tags:||added: floating-point library|
|Christophe Rhodes (csr21-cantab) wrote : [Christophe Rhodes] [Sbcl-commits] master: Make EXPT use double-precision throughout in more cases||#2|
|Christophe Rhodes (csr21-cantab) wrote : Re: [Bug 741564] [Christophe Rhodes] [Sbcl-commits] master: Make EXPT use double-precision throughout in more cases||#3|
|Changed in sbcl:|
|status:||Fix Committed → Fix Released|