From ed8dbe4c139f3428abbbda962d9bd243ab5eb777 Mon Sep 17 00:00:00 2001
From: Juan Jose Garcia Ripoll <jjgarcia@users.sourceforge.net>
Date: Fri, 20 Aug 2010 20:19:51 +0200
Subject: [PATCH] (EXPT 2 -2.0d0) is now computed in double precision.

---
 src/CHANGELOG    |  2 ++
 src/c/num_sfun.d | 84 ++++++++++++++++++++++++++++++------------------
 2 files changed, 55 insertions(+), 31 deletions(-)

diff --git a/src/CHANGELOG b/src/CHANGELOG
index 212516fc8..f2b43b6ca 100755
--- a/src/CHANGELOG
+++ b/src/CHANGELOG
@@ -94,6 +94,8 @@ ECL 10.5.1:
  - Function proclamations and declarations are also used to deduce the type
    of their arguments.
 
+ - (EXPT 2 -2.0d0) is now computed in double precision.
+
 ;;; Local Variables: ***
 ;;; mode:text ***
 ;;; fill-column:79 ***
diff --git a/src/c/num_sfun.d b/src/c/num_sfun.d
index 31d49684e..240d765b5 100644
--- a/src/c/num_sfun.d
+++ b/src/c/num_sfun.d
@@ -17,7 +17,8 @@
 
 #define ECL_INCLUDE_MATH_H
 #include <ecl/ecl.h>
-#include "ecl/internal.h"
+#include <ecl/internal.h>
+#include <ecl/ecl-inl.h>
 
 #ifndef HAVE_LOG1P
 double
@@ -160,49 +161,70 @@ cl_expt(cl_object x, cl_object y)
         @(return ecl_expt(x, y));
 }
 
+ecl_def_ct_single_float(singlefloat_one,1,static,const);
+ecl_def_ct_double_float(doublefloat_one,1,static,const);
+#ifdef ECL_LONG_FLOAT
+ecl_def_ct_long_float(longfloat_one,1,static,const);
+#endif
+
+static cl_object
+expt_zero(cl_object x, cl_object y)
+{
+	cl_type ty, tx;
+	cl_object z;
+        ty = type_of(y);
+        tx = type_of(x);
+        if (ecl_unlikely(!ECL_NUMBER_TYPE_P(tx))) {
+                FEwrong_type_nth_arg(@[expt], 1, x, @[number]);
+	}
+        /* INV: The most specific numeric types come first. */
+        switch ((ty > tx)? ty : tx) {
+        case t_fixnum:
+        case t_bignum:
+        case t_ratio:
+                return MAKE_FIXNUM(1);
+        case t_singlefloat:
+                return singlefloat_one;
+        case t_doublefloat:
+                return doublefloat_one;
+#ifdef ECL_LONG_FLOAT
+        case t_longfloat:
+                return longfloat_one;
+#endif
+        case t_complex:
+                z = expt_zero((tx == t_complex)? x->complex.real : x,
+                              (ty == t_complex)? y->complex.real : y);
+                return ecl_make_complex(z, MAKE_FIXNUM(0));
+        default:
+                /* We will never reach this */
+                (void)0;
+        }
+}
+
 cl_object
 ecl_expt(cl_object x, cl_object y)
 {
 	cl_type ty, tx;
 	cl_object z;
-        ty = type_of(y);
-        if (ecl_unlikely(!ECL_NUMBER_TYPE_P(ty))) {
-                FEwrong_type_nth_arg(@[expt], 2, y, @[number]);
+	if (ecl_unlikely(ecl_zerop(y))) {
+                return expt_zero(x, y);
 	}
+        ty = type_of(y);
         tx = type_of(x);
         if (ecl_unlikely(!ECL_NUMBER_TYPE_P(tx))) {
-                FEwrong_type_nth_arg(@[expt], 2, x, @[number]);
+                FEwrong_type_nth_arg(@[expt], 1, x, @[number]);
 	}
-	if (ecl_zerop(y)) {
-		/* INV: The most specific numeric types come first. */
-		switch ((ty > tx)? ty : tx) {
-		case t_fixnum:
-		case t_bignum:
-		case t_ratio:
-			z = MAKE_FIXNUM(1); break;
-		case t_singlefloat:
-			z = ecl_make_singlefloat(1.0); break;
-		case t_doublefloat:
-			z = ecl_make_doublefloat(1.0); break;
-#ifdef ECL_LONG_FLOAT
-		case t_longfloat:
-			z = ecl_make_longfloat(1.0); break;
-#endif
-		case t_complex:
-			z = ecl_expt((tx == t_complex)? x->complex.real : x,
-                                     (ty == t_complex)? y->complex.real : y);
-			z = ecl_make_complex(z, MAKE_FIXNUM(0));
-			break;
-		default:
-			/* We will never reach this */
-			(void)0;
-		}
-	} else if (ecl_zerop(x)) {
+        if (ecl_zerop(x)) {
 		z = ecl_times(x, y);
 		if (!ecl_plusp(ty==t_complex?y->complex.real:y))
 			z = ecl_divide(MAKE_FIXNUM(1), z);
 	} else if (ty != t_fixnum && ty != t_bignum) {
-		z = ecl_log1(x);
+                /* The following could be just
+                   z = ecl_log1(x);
+                   however, Maxima expects EXPT to have double accuracy
+                   when the first argument is integer and the second
+                   is double-float */
+		z = ecl_log1(ecl_times(x, expt_zero(x, y)));
 		z = ecl_times(z, y);
 		z = ecl_exp(z);
 	} else if (ecl_minusp(y)) {