(EXPT <notadouble> <complexdouble>) uses intermediate singlefloat values
Affects  Status  Importance  Assigned to  Milestone  

 SBCL 
Medium

Unassigned 
Bug Description
Error detected by Maxima's testsuite, 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 "floatzeta", defined in
maxima
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 doublefloat. These "expt"s are calculated directly using
SBCL's "expt". Here is a trace:
Maxima 5.23.1 http://
using Lisp SBCL 1.0.46.42
(%i1) to_lisp();
Type (tomaxima) to restart, ($quit) to quit Maxima.
MAXIMA> (trace expt)
(EXPT)
MAXIMA> (floatzeta #$(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))]
...
#C(1.107214406
The test case expects the following result:
#C(1.107214408
As a simple reduced test case without maxima we have:
(expt 2 #C(3.0d0 1.0d0))
SBCL gives:
#C(6.153911236
Clisp (2.48) gives:
#C(6.153911210
and Clisp with more precision (setf (system:
(expt 2 #C(3.0l0 1.0l0))
#C(6.153911210
5.111690210
The definition of "expt" is in sbcl/src/
(defun expt (base power)
...
(numberdispatch ((base number) (power number))
...
(((foreach fixnum (or bignum ratio) singlefloat doublefloat)
complex)
(if (and (zerop base) (plusp (realpart power)))
(* base power)
(exp (* power (log base)))))
...))
The calculation of (log base), base being an integer, yields a
singlefloat 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
12.1.4.1 "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 'doublefloat))))
then compiled SBCL and maxima. This SBCL calculates:
* (expt 2 #C(3.0d0 1.0d0))
#C(6.153911210
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 x8664, 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
singlefloat calculations if all arguments are integer or rational or
singlefloat, be it real or complex, and to use doublefloat calculations
if at least one argument is doublefloat 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 singlefloat, 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 doublefloat would need to be defined.
This may be easy as "log" already contains the needed logic, only
unfortunately always coercing the result to singlefloat.
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
anyway, currently.
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))
#C(1.0 5.421011e19)
whereas with the above hack we have the less exact
#C(1.0d0 0.0d0)
(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
whereas
(log (/ (1+ (expt 2 64)) (expt 2 64))) > 5.421011e20
which result would be correct for the former ratio, too.)
Greetings
Lutz
tags:  added: floatingpoint library 
Lutz Euler (lutzeuler) wrote :  #1 
Christophe Rhodes (csr21cantab) wrote : [Christophe Rhodes] [Sbclcommits] master: Make EXPT use doubleprecision throughout in more cases  #2 
status fixcommited
done
Applied the patches, thank you.
commit b90e13dea92ee66
Author: Lutz Euler <email address hidden>
Date: Fri Jul 1 18:06:17 2011 +0200
Make EXPT use doubleprecision throughout in more cases
lp#741564 notes that a Maxima test case fails because the result of
(EXPT <fixnum> <(complex double)>) is much less precise than expected.
This is caused by EXPT using an intermediate singlefloat value here.
This behaviour actually occurs for all the following combinations
of argument types:
(EXPT <(or rational singlefloat)> <(complex doublefloat)>)
(EXPT <(or (complex rational) (complex singlefloat))>
<(or (complex doublefloat) doublefloat)>)
In all these cases the first step EXPT does is to calculate (LOG BASE)
in single precision.
Refine the type dispatch clauses in EXPT to separate these cases
and coerce BASE to DOUBLEFLOAT or (COMPLEX DOUBLEFLOAT) there,
as appropriate, before applying LOG. Add tests.
Fixes lp#741564.
Signedoffby: Christophe Rhodes <email address hidden>
Christophe Rhodes (csr21cantab) wrote : Re: [Bug 741564] [Christophe Rhodes] [Sbclcommits] master: Make EXPT use doubleprecision throughout in more cases  #3 
status fixcommitted
done
> status fixcommited
Oops.
Changed in sbcl:  
status:  Fix Committed → Fix Released 
Here is a patch.
After thinking a lot more about what the standard may require andprecision/rational and doubleprecision arguments, I have come to
looking around what other math functions do in SBCL when given mixed
single
the conclusion that the way I want to fix this bug is indeed to convert
the BASE to DOUBLEFLOAT or (COMPLEX DOUBLEFLOAT), respectively, in
the cases in question, before taking the LOG. This does not overdo the
amount of changes and leads to the most consistent behaviour.
I believe that I have found all combinations of argument types of EXPT
that wrongly use the singleprecision LOG for the BASE. The patch
should correct EXPT's behaviour for all of them.
For the details see the commit messages and the comments in the code.
The patch consists of two commits:
The first commit patches EXPT and adds tests.
I have tested that it compiles and passes SBCL's test suite under linux
on x86_64 and x86. Also, maxima's test suite passes with no unexpected
errors (only tested on x86_64).
The second commit adds a compiletime check to NUMBERDISPATCH (moreNUMBERDISPATCH).NUMBERDISPATCH
precisely GENERATE
I was prompted to do that when the NUMBERDISPATCH form in EXPT, being
already quite complex, grew 50 percent more complex when I changed it
and when I immediately made an error there that didn't hinder the build,
that SBCL's test suite didn't catch but that did lead to lots of new
errors in maxima's test suite. The new check in GENERATE
would have caught this type of error earlier.
Greetings
Lutz