computing expontential as repeated multiplication for small integer exponent

An alternative is to use the Russian peasant algorithm. I removed the check
for Euler's constant. Maybe we should use expt for all cases and allow some
Maxima tests to fail for some versions of Lisp. I don't know what is best.
Using GCL, this code passes the 5.12 test suite.

;; Return x^n, where n is a nonnegative integer and x is a nonzero double
float.
;; This function uses the Russian peasant algorithm. Because x is nonzero,
we
;; don't have to trap 0.0^0. The first argument to n-power could be a CL
complex
;; double float (but this doesn't happen).

(defun n-power (x n)
  (let ((r) (z))
    (cond ((= 0 n) 1.0)
   (t
    (multiple-value-setq (n r) (floor n 2))
    (setq z (n-power x n))
    (* z z (if (= r 0) 1 x))))))

;; Return x^n, where x is either a double float or an integer. When x and n
are
;; both integers, use expt; when x is a double float, use the Russian
peasant algorithm.
;; We don't worry about 0^0 or 0.0^0 --- Maxima catches these cases before
getting
;; to exptb. For big floats, it seems that Maxima doesn't use this
function. But
;; maybe it should.

(defun exptb (x n)
  (if (integerp x) (cl-rat-to-maxima (expt x n))
    (if (< n 0) (n-power (/ 1 x) (- n)) (n-power x n))))

Barton

-----maxima-bounces at math.utexas.edu wrote: -----

>To: maxima <maxima at math.utexas.edu>
>From: "Robert Dodier" <robert.dodier at gmail.com>
>Sent by: maxima-bounces at math.utexas.edu
>Date: 07/29/2007 09:03PM
>Subject: computing expontential as repeated multiplication
>for  small integer exponent
>
>Hello,
>
>At present there is test suite case which fails (#69 in rtest8)
>for some Lisp implementations, apparently due to different
>results returned by the Lisp built-in function EXPT.
>
>I find that the attached patch for Maxima's EXPTB makes the
>result of rtest8 #69 come out the same for the Lisp
>implementations I have tested (SBCL, GCL, Clisp).
>Before, GCL and Clisp agreed and SBCL was different.
>
>This is somewhat cheesy so I thought I would gather 2nd
>opinions before committing. Comments? Maybe there is
>a better way to go about resolving this problem?
>
>best
>Robert Dodier
>
>PS.
>Index: src/simp.lisp
>===================================================================
>RCS file: /cvsroot/maxima/maxima/src/simp.lisp,v
>retrieving revision 1.41
>diff -u -r1.41 simp.lisp
>--- src/simp.lisp   4 Jul 2007 02:37:00 -0000   1.41
>+++ src/simp.lisp   30 Jul 2007 01:52:02 -0000
>@@ -652,14 +652,25 @@
>
> ;; I (rtoy) think EXPTB is meant to compute a^b, where b is an
> ;; integer.
>+
>+(defvar *exptb-max-integer-exponent* 16)
>+
> (defun exptb (a b)
>   (cond ((equal a %e-val)
>     ;; Make B a float so we'll get double-precision result.
>     (exp (float b)))
>+
>+    ;; Compute a^b as a repeated product when a is a float and b is a
>small positive integer.
>+    ;; This makes results for different Lisp implementations agree
>(otherwise EXPT differs).
>+    ((and (floatp a) (integerp b) (not (minusp b)) (< b
>*exptb-max-integer-exponent*))
>+     (let ((product 1))
>+       (dotimes (i b) (setq product (* product a)))
>+       product))
>+
>    ((or (floatp a) (not (minusp b)))
>     (expt a b))
>    (t
>-    (setq b (expt a (- b)))
>+    (setq b (exptb a (- b)))
>     (*red 1 b))))
>
> (defmfun simplus (x w z)       ; W must be 1
>_______________________________________________
>Maxima mailing list
>Maxima at math.utexas.edu
>http://www.math.utexas.edu/mailman/listinfo/maxima