Work on simpexpt



I have done some work on simpexpt to overcome some problems:

1. BUG ID: 2825082 "%pi^1.0b0 --> floating point value"

Simpexpt has special code to check the case %e^1.0b0 and returns a
bigfloat value. We can add code for the other numeric constants too.

((onep1 pot)
 ;; Exponent is One.
 (let ((y (mget gr '$numer)))
   (if (and y (floatp y) (or $numer (not (equal pot 1)))
       ;; Base is a numeric constant with a floating point value or 
       ;; $numer is TRUE and the Exponent is not the integer One.
       (return (if (and (member gr *builtin-numeric-constants*)
                        (equal pot bigfloatone))
                        ;; Convert numeric constant to bigfloat value.
                        ($bfloat gr)
                        ;; Return floating point value
                        y))
       ;; Handle other cases in exptrl.
       (return (exptrl gr pot)))))

Now all numeric constants evaluate accordingly:

(%i16) %pi^1.0b0;
(%o16) 3.141592653589793b0
(%i17) %gamma^1.0b0;
(%o17) 5.772156649015329b-1
(%i18) %phi^1.0b0;
(%o18) 1.618033988749895b0


2. BUG ID: 2825092 "%pi^2.0b0 does not evaluate numerically"

Maxima evaluates the numeric constants, if the power is a float number,
but not if a bigfloat number. Simpexpt only evaluates %e^2.0b0
numerically. We can add code to the label ATGR which handles the other
numeric constants too.

(t
 (let ((y (mget gr '$numer)))
   ;; Check for a numeric constant.
   (and y
        (floatp y)
        (or (floatp pot)
            ;; The exponent is a bigfloat. Convert base to bigfloat.
            (and ($bfloatp pot)
                 (member gr *builtin-numeric-constants*)
                 (setq y ($bfloat gr)))
                 (and $numer (integerp pot)))
            ;; The evaluation is done in exprtl.
            (return (exptrl y pot))))))

This are the results:

(%i12) %pi^2.0b0;
(%o12) 9.869604401089359b0
(%i13) %gamma^2.0b0;
(%o13) 3.331779238077187b-1


3. BUG ID: 660948 "simplification of exp(%i*...)"

Most of the problems of this bug report are no longer present. One open
problem is the simplification for float and bigfloat numbers of
expressions like exp(%i*%pi*x).

The simplification is done in the routine %especial. I think, the best
is to call this routine only for integer and rational numbers. This is
the code, which check this:

((and $%emode
      (mtimesp pot)
      (or (not (mnump (cadr pot)))
          (integerp (cadr pot))
          (ratnump (cadr pot)))
      (setq z (%especial pot)))
 (return z))

We no longer will get simplifications for float and bigfloat numbers,
but can evaluate the expression immediately to a float or bigfloat
number:

(%i19) exp(%i*%pi*0.25);
(%o19) %e^(0.25*%i*%pi)

(%i20) exp(%i*%pi*0.25),numer;
(%o20) 0.70710678118655*%i+0.70710678118655

(%i21) exp(%i*%pi*0.25b0);
(%o21) %e^(2.5b-1*%i*%pi)

(%i22) exp(%i*%pi*0.25b0),bfloat;
(%o22) 7.071067811865475b-1*%i+7.071067811865475b-1


4. BUG ID: 1010768 "sqrt(1/z) - 1/sqrt(z) => 0"

simpexpt does in general the wrong simplification:

(%i3) (1/z)^a;
(%o3) 1/z^a

That is for a=1/2 the reported bug. This simplification is true for a an
integer or z a positive real value. The following code in simpexpt is
responsible. In a first step I have removed all code which does the
wrong simplification (1/z)^a:

(cond ((or (eq $radexpand '$all)
           (simplexpon pot)             ; Maxima integer or odd rational
           (noneg (cadr gr))            ; base is positive
;           (equal (caddr gr) -1)       ; that is wrong
           (and (eq $domain '$real)     ; The following in general wrong
                (not (equal (caddr gr) -1)); we exclude only -1
                (odnump (caddr gr))))   ; rational with odd numerator
       (setq pot (mult pot (caddr gr))
                  gr (cadr gr)))
      ((and (eq $domain '$real) 
            (free gr '$%i) 
            $radexpand
            (not (decl-complexp (cadr gr))) 
            (evnump (caddr gr)))
       (setq pot (mult pot (caddr gr)) 
             gr (radmabs (cadr gr))))

; The following is in general wrong.
; The correct case with a positive base is already handled.
;           ((mminusp (caddr gr))
;            ;; Simplify (x^(-a))^b --> 1/(x^a)^b
;            (setq pot (neg pot)
;                  gr (list (car gr) (cadr gr) (neg (caddr gr)))))

           (t (go up)))

With this changes the case (1/z)^a will be correct:

(%i5) (1/z)^a;
(%o5) (1/z)^a

Simplification for an integer:

(%i6) declare(n,integer)$
(%i7) (1/z)^n;
(%o7) 1/z^n

Simplification for x a positive value:

(%i8) assume(x>0)$
(%i9) (1/x)^a;
(%o9) 1/x^a

The testsuite has no problems with the solutions to the problems 1, 2,
and 3. 

With the changes to the last problem 5 expected results have to be
modified, because sqrt(1/x) no longer is equal to 1/sqrt(x). One example
in rtest_taylor will fail completely. 

Further work on simpexpt is possible. But perhaps these are some useful
improvements.

Any suggestions?

Dieter Kaiser