Skip to content
Snippets Groups Projects
Commit 2600065b authored by ram's avatar ram
Browse files

Fixed EXPT to correctly return complex results for negative numbers

raised to fractional powers.
parent bbadce44
No related branches found
No related tags found
No related merge requests found
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
;;; Scott Fahlman (FAHLMAN@CMUC). ;;; Scott Fahlman (FAHLMAN@CMUC).
;;; ********************************************************************** ;;; **********************************************************************
;;; ;;;
;;; $Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.5 1990/10/24 16:42:48 ram Exp $ ;;; $Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.6 1991/01/03 13:16:27 ram Exp $
;;; ;;;
;;; This file contains all the irrational functions. Actually, most of the ;;; This file contains all the irrational functions. Actually, most of the
;;; work is done by calling out to C... ;;; work is done by calling out to C...
...@@ -96,6 +96,11 @@ ...@@ -96,6 +96,11 @@
(defparameter *intexp-maximum-exponent* 10000) (defparameter *intexp-maximum-exponent* 10000)
;;; This function precisely calculates base raised to an integral power. It
;;; separates the cases by the sign of power, for efficiency reasons, as powers
;;; can be calculated more efficiently if power is a positive integer. Values
;;; of power are calculated as positive integers, and inverted if negative.
;;;
(defun intexp (base power) (defun intexp (base power)
(when (> (abs power) *intexp-maximum-exponent*) (when (> (abs power) *intexp-maximum-exponent*)
(cerror "Continue with calculation." (cerror "Continue with calculation."
...@@ -113,46 +118,57 @@ ...@@ -113,46 +118,57 @@
(setq base (* base base)) (setq base (* base base))
(setq power nextn))))) (setq power nextn)))))
;;; This function calculates x raised to the nth power. It separates
;;; the cases by the type of n, for efficiency reasons, as powers can
;;; be calculated more efficiently if n is a positive integer, Therefore,
;;; All integer values of n are calculated as positive integers, and
;;; inverted if negative.
;;; EXPT -- Public
;;;
;;; If an integer power of a rational, use INTEXP above. Otherwise, do
;;; floating point stuff. If both args are real, we try %POW right off,
;;; assuming it will return 0 if the result may be complex. If so, we call
;;; COMPLEX-POW which directly computes the complex result. We also separate
;;; the complex-real and real-complex cases from the general complex case.
;;;
(defun expt (base power) (defun expt (base power)
"Returns BASE raised to the POWER." "Returns BASE raised to the POWER."
(if (zerop power) (if (zerop power)
;; This is wrong if power isn't an integer. (1+ (* base power))
(typecase (realpart base) (labels ((real-expt (base power rtype)
(single-float (coerce 1 'single-float)) (let* ((fbase (coerce base 'double-float))
(double-float (coerce 1 'double-float)) (fpower (coerce power 'double-float))
(t 1)) (res (coerce (%pow fbase fpower) rtype)))
(number-dispatch ((base number) (power number)) (if (and (zerop res) (minusp fbase))
(((foreach fixnum bignum ratio (complex rational)) integer) (multiple-value-bind (re im)
(intexp base power)) (complex-pow fbase fpower)
(((foreach single-float double-float) integer) (%make-complex (coerce re rtype) (coerce im rtype)))
(coerce (%pow (coerce base 'double-float) res)))
(coerce power 'double-float)) (complex-pow (fbase fpower)
'(dispatch-type base))) (let ((pow (%pow (- fbase) fpower))
(((foreach fixnum bignum ratio single-float) (fpower*pi (* fpower pi)))
(foreach ratio single-float)) (values (* pow (%cos fpower*pi))
(coerce (%pow (coerce base 'double-float) (* pow (%sin fpower*pi))))))
(coerce power 'double-float)) (declare (inline real-expt))
'single-float)) (number-dispatch ((base number) (power number))
(((foreach fixnum bignum ratio single-float double-float) double-float) (((foreach fixnum (or bignum ratio) (complex rational)) integer)
(%pow (coerce base 'double-float) (coerce power 'double-float))) (intexp base power))
(((complex rational) ratio) (((foreach single-float double-float) integer)
(* (expt (abs base) power) (real-expt base power '(dispatch-type base)))
(cis (* power (phase base))))) (((foreach fixnum (or bignum ratio) single-float)
(((complex float) (foreach integer ratio)) (foreach ratio single-float))
(* (expt (abs base) power) (real-expt base power 'single-float))
(cis (* power (phase base))))) (((foreach fixnum (or bignum ratio) single-float double-float)
(((foreach fixnum bignum ratio single-float double-float) complex) double-float)
(if (minusp base) (real-expt base power 'double-float))
(/ (exp (* power (log (- base))))) ((double-float single-float)
(exp (* power (log base))))) (real-expt base power 'double-float))
(((foreach (complex float) (complex rational)) complex) (((foreach (complex rational) (complex float)) rational)
(exp (* power (log base))))))) (* (expt (abs base) power)
(cis (* power (phase base)))))
(((foreach fixnum (or bignum ratio) single-float double-float)
complex)
(if (minusp base)
(/ (exp (* power (truly-the float (log (- base))))))
(exp (* power (truly-the float (log base))))))
(((foreach (complex float) (complex rational)) complex)
(exp (* power (log base))))))))
(defun log (number &optional (base nil base-p)) (defun log (number &optional (base nil base-p))
"Return the logarithm of NUMBER in the base BASE, which defaults to e." "Return the logarithm of NUMBER in the base BASE, which defaults to e."
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment