Make EXPT use double-precision 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 single-float value here.

This behaviour actually occurs for all the following combinations
of argument types:

  (EXPT <(or rational single-float)> <(complex double-float)>)

  (EXPT <(or (complex rational) (complex single-float))>
        <(or (complex double-float) double-float)>)

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 DOUBLE-FLOAT or (COMPLEX DOUBLE-FLOAT) there,
as appropriate, before applying LOG. Add tests.

Fixes lp#741564.

Signed-off-by: Christophe Rhodes <csr21@cantab.net>

diff --git a/NEWS b/NEWS
index 16323bc..6a48974 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -57,6 +57,10 @@ changes relative to sbcl-1.0.53:
   * bug fix: on 64-bit systems setting the nursery size above 4Gb now works.
     (lp#870868)
   * bug fix: SB-POSIX tests failed on OpenBSD 5.0. (lp#892707)
+  * bug fix: With several combinations of argument types, for example (EXPT
+    <integer> <(complex double)>), EXPT now uses double-precision throughout
+    instead of partially calculating only to single-precision.  (lp#741564;
+    thanks to Lutz Euler)
 
 changes in sbcl-1.0.53 relative to sbcl-1.0.52:
   * enhancement: on 64-bit targets, in src/compiler/generic/early-vm.lisp,

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 93c957b..1c684ac 100644 (file)
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -323,8 +323,12 @@
                                    (coerce (* pow (%cos y*pi))
                                            rtype)
                                    (coerce (* pow (%sin y*pi))
-                                           rtype)))))))))))))
-      (declare (inline real-expt))
+                                           rtype))))))))))))
+             (complex-expt (base power)
+               (if (and (zerop base) (plusp (realpart power)))
+                   (* base power)
+                   (exp (* power (log base))))))
+      (declare (inline real-expt complex-expt))
       (number-dispatch ((base number) (power number))
         (((foreach fixnum (or bignum ratio) (complex rational)) integer)
          (intexp base power))
@@ -338,19 +342,33 @@
          (real-expt base power 'double-float))
         ((double-float single-float)
          (real-expt base power 'double-float))
-        (((foreach (complex rational) (complex float)) rational)
+        ;; Handle (expt <complex> <rational>), except the case dealt with
+        ;; in the first clause above, (expt <(complex rational)> <integer>).
+        (((foreach (complex rational) (complex single-float)
+                   (complex double-float))
+          rational)
          (* (expt (abs base) power)
             (cis (* power (phase base)))))
-        (((foreach fixnum (or bignum ratio) single-float double-float)
-          complex)
-         (if (and (zerop base) (plusp (realpart power)))
-             (* base power)
-             (exp (* power (log base)))))
-        (((foreach (complex float) (complex rational))
+        ;; The next three clauses handle (expt <real> <complex>).
+        (((foreach fixnum (or bignum ratio) single-float)
+          (foreach (complex single-float) (complex rational)))
+         (complex-expt base power))
+        (((foreach fixnum (or bignum ratio) single-float)
+          (complex double-float))
+         (complex-expt (coerce base 'double-float) power))
+        ((double-float complex)
+         (complex-expt base power))
+        ;; The next three clauses handle (expt <complex> <float>) and
+        ;; (expt <complex> <complex>).
+        (((foreach (complex single-float) (complex rational))
+          (foreach (complex single-float) (complex rational) single-float))
+         (complex-expt base power))
+        (((foreach (complex single-float) (complex rational))
+          (foreach (complex double-float) double-float))
+         (complex-expt (coerce base '(complex double-float)) power))
+        (((complex double-float)
           (foreach complex double-float single-float))
-         (if (and (zerop base) (plusp (realpart power)))
-             (* base power)
-             (exp (* power (log base)))))))))
+         (complex-expt base power))))))
 
 ;;; FIXME: Maybe rename this so that it's clearer that it only works
 ;;; on integers?

diff --git a/tests/arith.pure.lisp b/tests/arith.pure.lisp
index 9f04fcf..9e0e782 100644 (file)
--- a/tests/arith.pure.lisp
+++ b/tests/arith.pure.lisp
@@ -475,3 +475,68 @@
 (with-test (:name :logbitp-wide-fixnum)
   (assert (not (logbitp (1- (integer-length most-positive-fixnum))
                         most-negative-fixnum))))
+
+;; EXPT dispatches in a complicated way on the types of its arguments.
+;; Check that all possible combinations are covered.
+(with-test (:name (:expt :argument-type-combinations))
+  (let ((numbers '(2                 ; fixnum
+                   3/5               ; ratio
+                   1.2f0             ; single-float
+                   2.0d0             ; double-float
+                   #c(3/5 1/7)       ; complex rational
+                   #c(1.2f0 1.3f0)   ; complex single-float
+                   #c(2.0d0 3.0d0))) ; complex double-float
+        (bignum (expt 2 64))
+        results)
+    (dolist (base (cons bignum numbers))
+      (dolist (power numbers)
+        (format t "(expt ~s ~s) => " base power)
+        (let ((result (expt base power)))
+          (format t "~s~%" result)
+          (push result results))))
+    (assert (every #'numberp results))))
+
+(with-test (:name :bug-741564)
+  ;; The bug was that in (expt <fixnum> <(complex double-float)>) the
+  ;; calculation was partially done only to single-float precision,
+  ;; making the complex double-float result too unprecise. Some other
+  ;; combinations of argument types were affected, too; test that all
+  ;; of them are good to double-float precision.
+  (labels ((nearly-equal-p (x y)
+             "Are the arguments equal to nearly double-float precision?"
+             (declare (type double-float x y))
+             (< (/ (abs (- x y)) (abs y))
+                (* double-float-epsilon 4))) ; Differences in the two least
+                                             ; significant mantissa bits
+                                             ; are OK.
+           (test-complex (x y)
+             (and (nearly-equal-p (realpart x) (realpart y))
+                  (nearly-equal-p (imagpart x) (imagpart y))))
+           (print-result (msg base power got expected)
+             (format t "~a (expt ~s ~s)~%got      ~s~%expected ~s~%"
+                     msg base power got expected)))
+    (let ((n-broken 0))
+      (flet ((test (base power coerce-to-type)
+               (let* ((got (expt base power))
+                      (expected (expt (coerce base coerce-to-type) power))
+                      (result (test-complex got expected)))
+                 (print-result (if result "Good:" "Bad:")
+                               base power got expected)
+                 (unless result
+                   (incf n-broken)))))
+        (dolist (base (list 2                    ; fixnum
+                            (expt 2 64)          ; bignum
+                            3/5                  ; ratio
+                            2.0f0))              ; single-float
+          (let ((power #c(-2.5d0 -4.5d0)))       ; complex double-float
+            (test base power 'double-float)))
+        (dolist (base (list #c(2.0f0 3.0f0)      ; complex single-float
+                            #c(2 3)              ; complex fixnum
+                            (complex (expt 2 64) (expt 2 65))
+                                                 ; complex bignum
+                            #c(3/5 1/7)))        ; complex ratio
+          (dolist (power (list #c(-2.5d0 -4.5d0) ; complex double-float
+                               -2.5d0))          ; double-float
+            (test base power '(complex double-float)))))
+      (when (> n-broken 0)
+        (error "Number of broken combinations: ~a" n-broken)))))