(in-package :maxima)

;;; This file over-writes a bunch of programs from various
;;; files to implement 1/0, 0/-1, 0/0   also  0/rat(0) etc.
;;; see paper c:/papers/infinity.tex

(defvar $ratinf t)
(defvar *ratsimpind '(rat simp))
(setf (get 'rat 'msimpind) *ratsimpind)  ;; make them all the same
(defvar *ratinf `(,*ratsimpind 1 0))
(defvar *ratind `(,*ratsimpind 0 0)) ;; indeterminate, unknown, ???

(defun $ratinf_p (x)(alike1 x *ratinf))
(defun $ratind_p (x)(alike1 x *ratind))



;; from simp.lisp
(defun zerop1 (x)
  (or (and (integerp x) (= 0 x))
      (and (floatp x) (= 0.0 x))
      (and ($bfloatp x) (= 0 (second x)))
      (and (ratnump x)(= (cadr x) 0) (not(= (caddr x) 0))) ;; added 1/2016 RJF for -0 as ((rat) 0 -1)
      ))


(defmfun simpexpt (x y z)
  (prog (gr pot check res rulesw w mlpgr mlppot)
     (setq check x)
     (cond (z (setq gr (cadr x) pot (caddr x)) (go cont)))
     (twoargcheck x)
     (setq gr (simplifya (cadr x) nil))
     (setq pot
           (let (($%enumer $numer))
             ;; Switch $%enumer on, when $numer is TRUE to allow 
             ;; simplification of $%e to its numerical value.
             (simplifya (if $ratsimpexpons ($ratsimp (caddr x)) (caddr x))
                        nil)))
  cont  
     (cond (($ratp pot)
            (setq pot (ratdisrep pot))
            (go cont))
           (($ratp gr)
            (cond ((member 'trunc (car gr) :test #'eq)
                   (return (srf (list '(mexpt) gr pot))))
                  ((integerp pot)
                   (let ((varlist (caddar gr)) (genvar (cadddr (car gr))))
                     (return (ratrep* (list '(mexpt) gr pot)))))
                  (t
                   (setq gr (ratdisrep gr))
                   (go cont))))
           ((or (setq mlpgr (mxorlistp gr))
                (setq mlppot (mxorlistp pot)))
            (go matrix))
           ((onep1 pot) (go atgr))
           ((or (zerop1 pot) (onep1 gr)) (go retno))
           
           ;; This code tries to handle 0^a more complete.
           ;; If the sign of realpart(a) is not known return an unsimplified
           ;; expression. The handling of the flag *zexptsimp? is not changed.
           ;; Reverting the return of an unsimplified 0^a, because timesin
           ;; can not handle such expressions. (DK 02/2010)
           ((zerop1 gr)
            (cond ((or (member (setq z ($csign pot)) '($neg $nz))
                       (and *zexptsimp? (eq ($asksign pot) '$neg)))
		   ;; A negative exponent. Maxima error.
		   
		   
                   (cond 
		    ($ratinf 
				 (return *ratinf))  ;; NEW LINE 12/2015 RJF
		    ((not errorsw) (merror (intl:gettext "expt: undefined: 0 to a negative exponent.")))
			 
			 
			 
                         (t (throw 'errorsw t))))
                  ((and (member z '($complex $imaginary))
                        ;; A complex exponent. Look at the sign of the realpart.
                        (member (setq z ($sign ($realpart pot))) 
                                '($neg $nz $zero)))
                   (cond ((not errorsw)
                          (merror (intl:gettext "expt: undefined: 0 to a complex exponent.")))
                         (t (throw 'errorsw t))))
                  ((and *zexptsimp? (eq ($asksign pot) '$zero))
                   (cond ((not errorsw)
                          (merror (intl:gettext "expt: undefined: 0^0")))
                         (t (throw 'errorsw t))))
                  ((not (member z '($pos $pz)))
                   ;; The sign of realpart(pot) is not known. We can not return
                   ;; an unsimplified 0^a expression, because timesin can not
                   ;; handle it. We return ZERO. That is the old behavior.
                   ;; Look for the imaginary symbol to be consistent with 
                   ;; old code.
                   (cond ((not (free pot '$%i))
                          (cond ((not errorsw)
                                 (merror (intl:gettext "expt: undefined: 0 to a complex exponent.")))
                                (t (throw 'errorsw t))))
                         (t
                          ;; Return ZERO and not an unsimplified expression.
                          (return (zerores gr pot)))))
                  (t (return (zerores gr pot)))))
           
           ((and (mnump gr)
                 (mnump pot)
                 (or (not (ratnump gr)) (not (ratnump pot))))
            (return (eqtest (exptrl gr pot) check)))
           ;; Check for numerical evaluation of the sqrt.
           ((and (alike1 pot '((rat) 1 2))
                 (or (setq res (flonum-eval '%sqrt gr))
                     (and (not (member 'simp (car x) :test #'eq))
                          (setq res (big-float-eval '%sqrt gr)))))
            (return res))
           ((eq gr '$%i)
            (return (%itopot pot)))
           ((and (realp gr) (minusp gr) (mevenp pot))
            (setq gr (- gr))
            (go cont))
           ((and (realp gr) (minusp gr) (moddp pot))
            (return (mul2 -1 (power (- gr) pot))))
           ((and (equal gr -1) (maxima-integerp pot) (mminusp pot))
            (setq pot (neg pot))
            (go cont))
           ((and (equal gr -1)
                 (maxima-integerp pot)
                 (mtimesp pot)
                 (= (length pot) 3)
                 (fixnump (cadr pot))
                 (oddp (cadr pot))
                 (maxima-integerp (caddr pot)))
            (setq pot (caddr pot))
            (go cont))
           ((atom gr) (go atgr))
           ((and (eq (caar gr) 'mabs)
                 (evnump pot)
                 (or (and (eq $domain '$real) (not (decl-complexp (cadr gr))))
                     (and (eq $domain '$complex) (decl-realp (cadr gr)))))
            (return (power (cadr gr) pot)))
           ((and (eq (caar gr) 'mabs)
                 (integerp pot)
                 (oddp pot)
                 (not (equal pot -1))
                 (or (and (eq $domain '$real) (not (decl-complexp (cadr gr))))
                     (and (eq $domain '$complex) (decl-realp (cadr gr)))))
            ;; abs(x)^(2*n+1) -> abs(x)*x^(2*n), n an integer number
            (if (plusp pot)
                (return (mul (power (cadr gr) (add pot -1))
                             gr))
                (return (mul (power (cadr gr) (add pot 1))
                             (inv gr)))))
           ((eq (caar gr) 'mequal)
            (return (eqtest (list (ncons (caar gr))
                                  (power (cadr gr) pot)
                                  (power (caddr gr) pot))
                            gr)))
           ((symbolp pot) (go opp))
           ((eq (caar gr) 'mexpt) (go e1))
           ((and (eq (caar gr) '%sum)
                 $sumexpand
                 (integerp pot)
                 (signp g pot)
                 (< pot $maxposex))
            (return (do ((i (1- pot) (1- i))
                         (an gr (simptimes (list '(mtimes) an gr) 1 t)))
                        ((signp e i) an))))
           ((equal pot -1) 
            (return (eqtest (testt (tms gr pot nil)) check)))
           ((fixnump pot)
            (return (eqtest (cond ((and (mplusp gr)
                                        (not (or (> pot $expop)
                                                 (> (- pot) $expon))))
                                   (expandexpt gr pot))
                                  (t (simplifya (tms gr pot nil) t)))
                            check))))
     
  opp
     (cond ((eq (caar gr) 'mexpt) (go e1))
           ((eq (caar gr) 'rat)
            (return (mul2 (power (cadr gr) pot)
                          (power (caddr gr) (mul2 -1 pot)))))
           ((not (eq (caar gr) 'mtimes)) (go up))
           ((or (eq $radexpand '$all) (and $radexpand (simplexpon pot)))
            (setq res (list 1))
            (go start))
           ((and (or (not (numberp (cadr gr)))
                     (equal (cadr gr) -1))
                 (equal -1 ($num gr)) ; only for -1
                 ;; Do not simplify for a complex base.
                 (not (member ($csign gr) '($complex $imaginary)))
                 (and (eq $domain '$real) $radexpand))
            ;; (-1/x)^a -> 1/(-x)^a for x negative
            ;; For all other cases (-1)^a/x^a
            (if (eq ($csign (setq w ($denom gr))) '$neg)
                (return (inv (power (neg w) pot)))
                (return (div (power -1 pot)
                             (power w pot)))))
           ((or (eq $domain '$complex) (not $radexpand)) (go up)))
     (return (do ((l (cdr gr) (cdr l)) (res (ncons 1)) (rad))
                 ((null l)
                  (cond ((equal res '(1))
                         (eqtest (list '(mexpt) gr pot) check))
                        ((null rad) 
                         (testt (cons '(mtimes simp) res)))
                        (t
                         (setq rad (power* ; RADEXPAND=()?
                                     (cons '(mtimes) (nreverse rad)) pot))
                         (cond ((not (onep1 rad))
                                (setq rad
                                      (testt (tms rad 1 (cons '(mtimes) res))))
                                (cond (rulesw
                                       (setq rulesw nil res (cdr rad))))))
                         (eqtest (testt (cons '(mtimes) res)) check))))
               ;; Check with $csign to be more complete. This prevents wrong 
               ;; simplifications like sqrt(-z^2)->%i*sqrt(z^2) for z complex.
               (setq z ($csign (car l)))
               (if (member z '($complex $imaginary))
                   (setq z '$pnz)) ; if appears complex, unknown sign
               (setq w (cond ((member z '($neg $nz) :test #'eq)
                              (setq rad (cons -1 rad))
                              (mult -1 (car l)))
                             (t (car l))))
               (cond ((onep1 w))
                     ((alike1 w gr) (return (list '(mexpt simp) gr pot)))
                     ((member z '($pn $pnz) :test #'eq)
                      (setq rad (cons w rad)))
                     (t
                      (setq w (testt (tms (simplifya (list '(mexpt) w pot) t)
                                          1 (cons '(mtimes) res))))))
               (cond (rulesw (setq rulesw nil res (cdr w))))))
     
  start
     (cond ((and (cdr res) (onep1 (car res)) (ratnump (cadr res)))
            (setq res (cdr res))))
     (cond ((null (setq gr (cdr gr)))
            (return (eqtest (testt (cons '(mtimes) res)) check)))
           ((mexptp (car gr))
            (setq y (list (caar gr) (cadar gr) (mult (caddar gr) pot))))
           ((eq (car gr) '$%i)
            (setq y (%itopot pot)))
           ((mnump (car gr))
            (setq y (list '(mexpt) (car gr) pot)))
           (t (setq y (list '(mexpt simp) (car gr) pot))))
     (setq w (testt (tms (simplifya y t) 1 (cons '(mtimes) res))))
     (cond (rulesw (setq rulesw nil res (cdr w))))
     (go start)
     
  retno
     (return (exptrl gr pot))
     
  atgr
     (cond ((zerop1 pot) (go retno))
           ((onep1 pot)
            (let ((y (mget gr '$numer)))
              (if (and y (floatp y) (or $numer (not (equal pot 1))))
                  ;; A numeric constant like %e, %pi, ... and 
                  ;; exponent is a float or bigfloat value.
                  (return (if (and (member gr *builtin-numeric-constants*)
                                   (equal pot bigfloatone))
                              ;; Return a bigfloat value.
                              ($bfloat gr)
                              ;; Return a float value.
                              y))
                  ;; In all other cases exptrl simplifies accordingly.
                  (return (exptrl gr pot)))))
           ((eq gr '$%e)
	    
	    (when (and $ratinf (ratnump pot)(= 0 (caddr pot))) ;; exp(0/0), 1/0, -1/0  RJF
	     (return (case (cadr pot) (0 *ratind)(-1 0)(1 *ratinf))))
            ;; Numerically evaluate if the power is a flonum.
            (when $%emode
              (let ((val (flonum-eval '%exp pot)))
                (when val
                  (return val)))
              ;; Numerically evaluate if the power is a (complex)
              ;; big-float.  (This is basically the guts of
              ;; big-float-eval, but we can't use big-float-eval.)
              (when (and (not (member 'simp (car x) :test #'eq))
                         (complex-number-p pot 'bigfloat-or-number-p))
                (let ((x ($realpart pot))
                      (y ($imagpart pot)))
                  (cond ((and ($bfloatp x) (like 0 y))
                         (return ($bfloat `((mexpt simp) $%e ,pot))))
                        ((or ($bfloatp x) ($bfloatp y))
                         (let ((z (add ($bfloat x) (mul '$%i ($bfloat y)))))
                           (setq z ($rectform `((mexpt simp) $%e ,z)))
                           (return ($bfloat z))))))))
            (cond ((and $logsimp (among '%log pot)) (return (%etolog pot)))
                  ((and $demoivre (setq z (demoivre pot))) (return z))
                  ((and $%emode
                        (among '$%i pot)
                        (among '$%pi pot)
                        ;; Exponent contains %i and %pi and %emode is TRUE:
                        ;; Check simplification of exp(%i*%pi*p/q*x)
                        (setq z (%especial pot)))
                   (return z))
                  (($taylorp (third x))
                   ;; taylorize %e^taylor(...)
                   (return ($taylor x)))))
           (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)))
                   (return (exptrl y pot))))))

  up 
     (return (eqtest (list '(mexpt) gr pot) check))

  matrix
     (cond ((zerop1 pot)
            (cond ((mxorlistp1 gr) (return (constmx (addk 1 pot) gr)))
                  (t (go retno))))
           ((onep1 pot) (return gr))
           ((or $doallmxops $doscmxops $domxexpt)
            (cond ((or (and mlpgr
                            (or (not ($listp gr)) $listarith)
                            (scalar-or-constant-p pot $assumescalar))
                       (and $domxexpt
                            mlppot
                            (or (not ($listp pot)) $listarith)
                            (scalar-or-constant-p gr $assumescalar)))
                   (return (simplifya (outermap1 'mexpt gr pot) t)))
                  (t (go up))))
           ((and $domxmxops (member pot '(-1 -1.0) :test #'equal))
            (return (simplifya (outermap1 'mexpt gr pot) t)))
           (t (go up)))
  e1 
     ;; At this point we have an expression: (z^a)^b with gr = z^a and pot = b
     (cond ((or (eq $radexpand '$all)
                ;; b is an integer or an odd rational
                (simplexpon pot)
                (and (eq $domain '$complex)
                     (not (member ($csign (caddr gr)) '($complex $imaginary)))
                         ;; z >= 0 and a not a complex
                     (or (member ($csign (cadr gr)) '($pos $pz $zero))
                         ;; -1 < a <= 1
                         (and (mnump (caddr gr))
                              (eq ($sign (sub 1 (take '(mabs) (caddr gr))))
                                  '$pos))))
                (and (eq $domain '$real)
                     (member ($csign (cadr gr)) '($pos $pz $zero)))
                ;; (1/z)^a -> 1/z^a when z a constant complex
                (and (eql (caddr gr) -1)
                     (or (and $radexpand
                              (eq $domain '$real))
                         (and (eq ($csign (cadr gr)) '$complex)
                              ($constantp (cadr gr)))))
                ;; This does (1/z)^a -> 1/z^a. This is in general wrong.
                ;; We switch this type of simplification on, when
                ;; $ratsimpexpons is T. E.g. radcan sets this flag to T.
                ;; radcan hangs for expressions like sqrt(1/(1+x)) without
                ;; this simplification.
                (and $ratsimpexpons
                     (equal (caddr gr) -1))
                (and $radexpand
                     (eq $domain '$real)
                     (odnump (caddr gr))))
            ;; Simplify (z^a)^b -> z^(a*b)
            (setq pot (mul pot (caddr gr))
                  gr (cadr gr)))
           ((and (eq $domain '$real)
                 (free gr '$%i)
                 $radexpand
                 (not (decl-complexp (cadr gr)))
                 (evnump (caddr gr)))
            ;; Simplify (x^a)^b -> abs(x)^(a*b)
            (setq pot (mul pot (caddr gr))
                  gr (radmabs (cadr gr))))
           ((and $radexpand
                 (eq $domain '$real)
                 (mminusp (caddr gr)))
            ;; Simplify (1/z^a)^b -> 1/(z^a)^b
            (setq pot (neg pot)
                  gr (power (cadr gr) (neg (caddr gr)))))
           (t (go up)))
    (go cont)))


(defun addk (x y)
  (cond ((eql x 0) y)
	((eql y 0) x)
	((and (numberp x) (numberp y))
         (cl-or-bfloat-binary-op mplus + x y))
	((or ($bfloatp x) ($bfloatp y)) ($bfloat (list '(mplus) x y)))
	(t (prog (g a b)
	      (cond ((numberp x)
		     (cond ((floatp x) (return (+ x (fpcofrat y))))
			   (t (setq x (list '(rat) x 1)))))
		    ((numberp y)
		     (cond ((floatp y) (return (+ y (fpcofrat x))))
			   (t (setq y (list '(rat) y 1))))))
	     (setq g (gcd (caddr x) (caddr y)))
	     
	     (cond ((= (caddr x) 0)
		    (cond((= (cadr x) 0)(return x))
			 ((= (caddr y) 0)
			  (return (if (equal x y) x ; oo+oo or -oo-oo
				    *ratind))))) ;oo-oo is und
		   ((= (caddr y) 0) ; return oo or -oo or und		   
		    (return y)))
	     
	    ;; above clause only change 12/2015 rjf
	     
	     (setq a (truncate (caddr x) g)
		   b (truncate (caddr y) g))

	     
	      (return (timeskl (list '(rat) 1 g)
			       (list '(rat)
				     (+ (* (cadr x) b)
					   (* (cadr y) a))
				     (* a b))))))))

(defun partition_ex_lisp   ;; hacked, simplified for this purpose
    (pred e)
  (let ((yes nil)(no nil))
    (map nil
	     #'(lambda(r) 
		 (if (funcall pred r)
                      (setf yes (cons r yes))
                                ; r fails predicate
                          (setf no (cons r no))))
               e)
           (cons yes no)))


(defun simptimes_ratinf (x)
  ;; here we combine floats, integers, ((rat) n d), bfloats, ((rat) n 0)
  ;; maybe common lisp  rats, complexes, intervals.
  ;; into one item. e.g.  ((mtimes) 0 ((rat) 1 0)  $r $s)
  ;; becomes ((mtimes)((rat) 0 0)  $r $s)
  ;; phew 
  ;;first pass, look for ANY  ((rat *) * 0)
  (if (notany #'(lambda(s) ;;(format t"~%s=~s" s)
		 (and (consp s)
		      (or (and(eq (caar s) 'rat)
			      (= (caddr s) 0)) ; denom is zero
			  (and(eq (caar s) 'mrat) ;; ((mrat) n .  0)
			      (equal (cddr s) 0)))
		      
		      ))
	      (cdr x))
     x ;; no rational infinities. quickest check, usually comes out here
    (let ((h (partition_ex_lisp  #'(lambda(r)(or (mnump r)
						 (and (consp r)
						      (eq (caar r)'mrat))))
				 (cdr x))))
   ; (format t "~% (cdar h) is ~s" (cdar h))
    (if (null(cdar h)) x		; zero or only one item in number part
					; just return same value
      (infprod_hackit (car h)(cdr h))))))

;; nums is a list of 2 or more items.
;; look for at least one 0/0  --> 0.
;; look for at least one n/0 * 0 --> 0/0.
;;  (-1/0)^n  is  1/0 if n is even.  but  -3* (1/0) is -1/0 also.
;; need to check sign of all numbers.

(defun ratinfp(s)(and (consp s)
		   (eq (caar s) 'rat)
		   (equal (caddr s) 0)))

(defun infprod_hackit (nums rest) 
  (let ((poscount 0)(indefcount 0) (zerocount 0) (sign 0) (mrat nil))
    ;;counts of +oo -oo  0/0 0
    ;; mrat = t if the item encountered is mrat encoded not rat.
    ;; mrat allows for encoding -0, as one bonus.
    ;;sign
  (loop for s in nums do
	(cond((ratinfp s) ;denom is zero, look at numerator
	      (case (cadr s) ; numerator of ((rat) n d) is cadr
		(1 (incf poscount))
		(-1 (incf sign))
		(0 (incf indefcount))))
	     ;; the mrat case
	     ((and (consp s)
		   (eq (caar s) 'mrat)
		   (equal (cddr s) 0)) ;denom is zero, look at numerator
	      (setq mrat t)
	      (case (cadr s) ;; the numerator of ((mrat) . ( n . d))
		(1 (incf poscount))
		(-1 (incf sign))
		(0 (incf indefcount))))
	     
	     ;; not an infinity or indefinite
	     ((zerop1 s)		; 0 or 0.0 or 0.0b0
	      (incf zerocount))
	     ;leftover numbers, possibly bigfloat
	     ((and (mnump s)(mgrp 0 s)) (incf sign))))
  
;; (format t "~% pos=~s  indef=~s zer=~s sign=~s"	  poscount indefcount zerocount sign)
;;  here we make it an mrat if we encountered an mrat inf or 0/0
  (if mrat
     (list '(mtimes) (cons '(mrat simp nil nil) ; return mrat of oo -oo or 0/0
	    (cons
					; new numerator is..
	     (if (or (> indefcount 0)(> zerocount 0)) ;0/0 or 0* infinity
		 0 
	       (expt -1 sign))
					; (-1/0)^negcount
	     0)))

      (cons '(mtimes) ; return product: (rat of oo -oo or 0/0)* rest
	(cons
	 (cond ((or (> indefcount 0)(> zerocount 0)) ;0/0 or 0* infinity
		*ratind)
	       (t (list *ratsimpind (expt -1 sign) 0))) ; (-1/0)^negcount
	 rest)))))

  
(defmfun simptimes (x w z)			; W must be 1
  (prog (res check eqnflag matrixflag sumflag)
     (if (null (cdr x)) (return 1))
    (setq check x)
    
    (if $ratinf (setf x (simptimes_ratinf x)))
    
  start
    (setq x (cdr x))
    
    (cond ;;added 1/2016 .
    
     #+ignore
     ((and(consp  res) 
	;  (prog nil (format t "~%start res=~s x=~s" res x) t)
	  (member 0 (cdr res) :test #'(lambda(r s) ; car of res is (mtimes)
			   (and (consp s)
				(eq (caar s) 'rat)
				(= (caddr s) 0) ; denom is zero
				)))
	  (member 0 x))
     
    ;  (format t "~%DONE! res=~s x=~s" res x)
      (if (member 0 x) (return *ratind)  ;; 0*(1/0).  ;;; not there yet
	))
     
     
     
     ((zerop1 res) 
	   ;;; added 1/2016 RJF
      #+ignore
      (if (member 0 x :test #'(lambda(r s)
			   (and (consp s)
				(eq (caar s) 'rat)
				(= (caddr s) 0) ; denom is zero
				))) (return *ratind))

	    (cond ($mx0simp
		   (cond ((and matrixflag (mxorlistp1 matrixflag))
			  (return (constmx res matrixflag)))
			 (eqnflag (return (list '(mequal simp)
						(mul2 res (cadr eqnflag))
						(mul2 res (caddr eqnflag)))))
		         (t
		          (dolist (u x)
			    (cond ((mxorlistp u)
				   (return (setq res (constmx res u))))
				  ((and (mexptp u)
					(mxorlistp1 (cadr u))
					($numberp (caddr u)))
				   (return (setq res (constmx res (cadr u)))))
				  ((mequalp u)
				   (return
				     (setq res 
				           (list '(mequal simp)
						 (mul2 res (cadr u))
						 (mul2 res (caddr u))))))))))))
	    (return res))
	   ((null x) (go end)))
     (setq w (if z (car x) (simplifya (car x) nil)))
  st1
     (cond ((atom w) nil)
	   ((eq (caar w) 'mrat)
	    (cond ((or eqnflag matrixflag
	               (and sumflag
	                    (not (member 'trunc (cdar w) :test #'eq)))
		       (spsimpcases (cdr x) w))
	           (setq w (ratdisrep w))
	           (go st1))
	          (t
	           (return 
	             (ratf (cons '(mtimes)
			         (nconc (mapcar #'simplify (cons w (cdr x)))
					(cdr res))))))))
	   ((eq (caar w) 'mequal)
	    (setq eqnflag
		  (if (not eqnflag)
		      w
		      (list (car eqnflag)
			    (mul2 (cadr eqnflag) (cadr w))
			    (mul2 (caddr eqnflag) (caddr w)))))
	    (go start))
	   ((member (caar w) '(mlist $matrix) :test #'eq)
	    (setq matrixflag
		  (cond ((not matrixflag) w)
			((and (or $doallmxops $domxmxops $domxtimes)
			      (or (not (eq (caar w) 'mlist)) $listarith)
			      (not (eq *inv* '$detout)))
			 (stimex matrixflag w))
			(t (setq res (tms w 1 res)) matrixflag)))
	    (go start))
	   ((and (eq (caar w) '%sum) $sumexpand)
	    (setq sumflag (sumtimes sumflag w))
	    (go start)))
     (setq res (tms w 1 res))
     (go start)
  end
     (cond ((mtimesp res) (setq res (testt res))))
     (cond (sumflag (setq res (cond ((or (null res) (equal res 1)) sumflag)
				    ((not (mtimesp res))
				     (list '(mtimes) res sumflag))
				    (t (nconc res (list sumflag)))))))
     (cond ((or (atom res)
		(not (member (caar res) '(mexpt mtimes) :test #'eq))
		(and (zerop $expop) (zerop $expon))
		expandflag))
	   ((eq (caar res) 'mtimes) (setq res (expandtimes res)))
	   ((and (mplusp (cadr res))
		 (fixnump (caddr res))
		 (not (or (> (caddr res) $expop)
			  (> (- (caddr res)) $expon))))
	    (setq res (expandexpt (cadr res) (caddr res)))))
     (cond (matrixflag
            (setq res
                  (cond ((null res) matrixflag)
                        ((and (or ($listp matrixflag)
                                  $doallmxops
			          (and $doscmxops
			               (not (member res '(-1 -1.0) :test #'equal)))
			          ;; RES should only be -1 here (not = 1)
                                  (and $domxmxops
                                       (member res '(-1 -1.0) :test #'equal)))
                              (or (not ($listp matrixflag)) $listarith))
                         (mxtimesc res matrixflag))
			(t (testt (tms matrixflag 1 (tms res 1 nil))))))))
     (if res (setq res (eqtest res check)))
     (return (cond (eqnflag
		    (if (null res) (setq res 1))
		    (list (car eqnflag)
			  (mul2 (cadr eqnflag) res)
			  (mul2 (caddr eqnflag) res)))
		   (t res)))))


;; we do   need this change.??
(defun *red (n d)
  (cond ((zerop n)(if $ratinf (list *ratsimpind
				    0 (if (numberp d)(cond ((< d 0) -1)((= d 0) 0)(t 1))
				    0)) 0))	; 0/x is 0. 0/-5 is -0 1/2015 rjf
	((equal d 1) n)
	((equal d 0) (list *ratsimpind (signum n) d)) ;; changed 12/2015 rjf
	(t (let ((u (gcd n d)))
	     (setq n (truncate n u)
	           d (truncate d u))
	     (if (minusp d) (setq n (- n) d (- d)))
	     (cond ((equal d 1) n)
		   ($float (fpcofrat1 n d))
		   (t (list *ratsimpind n d)))))))

(defun simpmin (x vestigial z)
  (declare (ignore vestigial))
  (cond ((null (cdr x)) 0)
	((equal 0 (cadr x)) `(,*ratsimpind 0 -1))  ;rjf 1/2016  -0
        ((null (cddr x))
	 (mul -1 (simplifya (cadr x) z)))
        (t
         ;; ((mminus) a b ...) -> ((mplus) a ((mtimes) -1 b) ...)
         (sub (simplifya (cadr x) z) (addn (cddr x) z)))))

(defmfun simpquot (x y z)
  (twoargcheck x)
  (cond ((and (integerp (cadr x)) (integerp (caddr x)) (not (zerop (caddr x))))
	 (*red (cadr x) (caddr x)))
	((and (numberp (cadr x)) (numberp (caddr x)) (not (zerop (caddr x))))
         (cl-or-bfloat-binary-op mquotient / (cadr x) (caddr x)))
	(t (setq y (simplifya (cadr x) z))
	   (setq x (simplifya (list '(mexpt) (caddr x) -1) z))
	   (if (equal y 1) x (simplifya (list '(mtimes) y x) t)))))


(defmfun exptrl (r1 r2)  ;compute r1^r2, both numbers of some sort
  (cond ((equal r2 1) r1)
        ((equal r2 1.0) 
         (cond ((mnump r1) (addk 0.0 r1))
               ;; Do not simplify the type of the number away.
               (t (list '(mexpt simp) r1 1.0))))
        ((equal r2 bigfloatone)
         (cond ((mnump r1) ($bfloat r1))
               ;; Do not simplify the type of the number away.
               (t (list '(mexpt simp) r1 bigfloatone))))
	((zerop1 r1)
	 
	 (cond 
	  ;; changed 1/2016 RJF
	  ((or (zerop1 r2) (mnegp r2))
	   (cond ($ratinf *ratind) ;; this is RJF change
		 ((not errorsw)
		    (merror (intl:gettext "expt: undefined: ~M") (list '(mexpt) r1 r2)))
		   (t (throw 'errorsw t))))
	  (t (zerores r1 r2))))
	((alike1 r2 *ratinf) (mul (signum r1) *ratinf)) ;new rjf
	((alike1 r2 *ratind) *ratind)	;new rjf
	((alike1 r1 *ratind) *ratind)	;new rjf

	((or (zerop1 r2) (onep1 r1))
	 (cond ((or ($bfloatp r1) ($bfloatp r2)) bigfloatone)
	       ((or (floatp r1) (floatp r2)) 1.0)
	       (t 1)))
	((or ($bfloatp r1) ($bfloatp r2)) ($bfloat (list '(mexpt) r1 r2)))
	((and (numberp r1) (integerp r2)) (exptb r1 r2))
	((and (numberp r1) (floatp r2) (equal r2 (float (floor r2))))
	 (exptb (float r1) (floor r2)))
	((or $numer (and (floatp r2) (or (plusp (num1 r1)) $numer_pbranch)))
	 (let (y  #+kcl(r1 r1) #+kcl(r2 r2))
	   (cond ((mnegp (setq r1 (addk 0.0 r1)))
		  (cond ((or $numer_pbranch (eq $domain '$complex))
		         ;; for R1<0:
		         ;; R1^R2 = (-R1)^R2*cos(pi*R2) + i*(-R1)^R2*sin(pi*R2)
			 (setq r2 (addk 0.0 r2))
			 (setq y (exptrl (- r1) r2) r2 (* %pi-val r2))
			 (add2 (* y (cos r2))
			       (list '(mtimes simp) (* y (sin r2)) '$%i)))
			(t (setq y (let ($numer $float $keepfloat $ratprint)
				     (power -1 r2)))
			   (mul2 y (exptrl (- r1) r2)))))
	         ((equal (setq r2 (addk 0.0 r2)) (float (floor r2)))
	          (exptb r1 (floor r2)))
	         ((and (equal (setq y (* 2.0 r2)) (float (floor y)))
	               (not (equal r1 %e-val)))
		  (exptb (sqrt r1) (floor y)))
                 ;; If bfloat conversion turned r1 or r2 into a bigfloat, we
                 ;; can't call cl:exp or cl:log here, but we *can* safely
                 ;; recurse (which will usually give a noun form that then gets
                 ;; simplified)
                 ((or ($bfloatp r1) ($bfloatp r2))
                  (exptrl r1 r2))
		 (t
                  (exp (* r2 (log r1)))))))
	((floatp r2) (list '(mexpt simp) r1 r2))
	((integerp r2) ;; here we are left with rat^power RJF
	 (cond ((minusp r2)
		;; this looks dubious too RJF. 1/2016. I'd like (-1/0)^-1
		;; to come out as -0. and (-0)^-1 to be (-1/0).i.e. -oo.
		;; some lisps may have (* -1 0) be -0. which would be right.
		;; but may never happen in others. most others. RJF
		(if

		    (and $ratinf 
		      (= 0 (caddr r1))) ;; (n/0) ^ -1
			 #+ignore
		  '((rat) 0 0)		;rjf
		  (list '(rat)0 (cadr r1) )  ; 0/n
		
	        (exptrl (cond ((equal (abs (cadr r1)) 1)
	                       (* (cadr r1) (caddr r1)))
	                       ;; We set the simp flag at this place. This
	                       ;; changes nothing for an exponent r2 # -1.
	                       ;; exptrl is called again and does not look at
	                       ;; the simp flag. For the case r2 = -1 exptrl
	                       ;; is called with an exponent 1. For this case
	                       ;; the base is immediately returned. Now the
	                       ;; base has the correct simp flag. (DK 02/2010)
			      ((minusp (cadr r1))
			       (list *ratsimpind (- (caddr r1)) (- (cadr r1))))
			      (t (list *ratsimpind (caddr r1) (cadr r1))))
			(- r2))))
	       (t (list *ratsimpind (exptb (cadr r1) r2) (exptb (caddr r1) r2)))))

	((and (floatp r1) (alike1 r2 '((rat) 1 2)))
	 ;; this all looks dubious-- RJF.  There are 2 square roots. 1/2016
	 (if (minusp r1)
	     (list '(mtimes simp) (sqrt (- r1)) '$%i)
	     (sqrt r1)))
	((and (floatp r1) (alike1 r2 '((rat) -1 2)))
	 (if (minusp r1)
	     (list '(mtimes simp) (/ -1.0 (sqrt (- r1))) '$%i)
	     (/ (sqrt r1))))
	((floatp r1)
	 (if (plusp r1)
	     (exptrl r1 (fpcofrat r2))
	     (mul2 (exptrl -1 r2) ;; (-4.5)^(1/4) -> (4.5)^(1/4) * (-1)^(1/4)
		   (exptrl (- r1) r2))))
	(exptrlsw (list '(mexpt simp) r1 r2))
	(t
	 (let ((exptrlsw t))
	   (simptimes (list '(mtimes)
			    (exptrl r1 (truncate (cadr r2) (caddr r2)))
			    (let ((y (let ($keepfloat $ratprint)
				       (simpnrt r1 (caddr r2))))
				  (z (rem (cadr r2) (caddr r2))))
			      (if (mexptp y)
				  (list (car y) (cadr y) (mul2 (caddr y) z))
				  (power y z))))
		      1 t)))))















;;from nforma.lisp

#+ignore
(setf (get 'rat 'formatter) 
  #'(lambda(form)(cond ((minusp (cadr form))
			(list '(mminus) (list '(rat) (- (cadr form)) (caddr form))))
		       (t (cons '(rat) (cdr form))))))
(setf (get 'rat 'formatter) 
  #'(lambda(form)(cond 

		       
		       ;; next clauses 1/2016, Rjf
		       ((and (= 0 (cadr form))(= -1 (caddr form))) ;-0
			(list '(mminus) 0))
		       ((= 0 (caddr form)) ; inf? indef?
			(case (cadr form)
			  (1 '$oo)
			  (-1 '$\-oo)
			  (0 '$indef)))
		       ((minusp (cadr form)) ;neg numerator
			(list '(mminus) (list '(rat) (- (cadr form)) (caddr form))))
			
		       (t (cons '(rat) (cdr form))))))




;;; from trigi.lisp
;;; example of how to do ratinf
(defmfun simp-%sin (form y z) 
  (oneargcheck form)
  (setq y (simpcheck (cadr form) z))
  (cond ((flonum-eval (mop form) y))
	((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
	((taylorize (mop form) (second form)))
	((and $%piargs (cond ((zerop1 y) 0)
			     ((has-const-or-int-term y '$%pi) (%piargs-sin/cos y)))))
	((and $%iargs (multiplep y '$%i)) (mul '$%i (cons-exp '%sinh (coeff y '$%i 1))))
	((and $triginverses (not (atom y))
	      (cond ((eq '%asin (setq z (caar y))) (cadr y))
		    ((eq '%acos z) (sqrt1-x^2 (cadr y)))
		    ((eq '%atan z) (div (cadr y) (sqrt1+x^2 (cadr y))))
		    ((eq '%acot z) (div 1 (sqrt1+x^2 (cadr y))))
		    ((eq '%asec z) (div (sqrtx^2-1 (cadr y)) (cadr y)))
		    ((eq '%acsc z) (div 1 (cadr y)))
		    ((eq '$atan2 z) (div (cadr y) (sq-sumsq (cadr y) (caddr y)))))))
	((and $trigexpand (trigexpand '%sin y)))
	($exponentialize (exponentialize '%sin y))
	((and $halfangles (halfangle '%sin y)))
	((apply-reflection-simp (mop form) y $trigsign))
					;((and $trigsign (mminusp* y)) (neg (cons-exp '%sin (neg y))))
	((and $ratinf (ratnump y)(= 0 (caddr y)))  ;; added for 0/0 etc  RJF 1/2016
	 '((rat) 0 0))		; denominator is zero. numerator 0,1,-1 -- all same value of sin()

	(t (eqtest (list '(%sin) y) form))))


;;; from compar.lisp

(defmfun infsimp* (e) ;; redone entirely RJF 1/2016
  
  (cond
       ((and $ratinf (ratnump e)(= (caddr e) 0)) (cadr e)
	;(case (cadr e)(0 '$und)(1 '$pos)(-1 '$neg))
	)
       ((or (atom e) (and (free e '$inf) (free e '$minf)))
	e)
    (t (infsimp e))))

(defmfun simpsignum (e y z)
  (declare (ignore y))
  (oneargcheck e)
  (if(and $ratinf (ratnump (cadr e))
		(= 0 (caddr (cadr e)))) ;; rjf sign of n/0.  
	   ;;sign of 0/0 is ??
	   (if (= 0 (cadr (cadr e))) e
	     (cadr (cadr e)))
  (let ((x (simpcheck (second e) z)) (sgn))
    (cond 

	  ((complex-number-p x #'mnump)
		    (if (complex-number-p x #'$ratnump) ;; nonfloat complex
		        (if (zerop1 x) 0 ($rectform (div x ($cabs x))))
		      (maxima::to (bigfloat::signum (bigfloat::to x)))))
		   
	  ;; idempotent: signum(signum(z)) = signum(z).
	  ((and (consp x) (consp (car x)) (eq '%signum (mop x))) x)
		   
	  (t
	   (setq sgn ($csign x))
	   (cond ((eq sgn '$neg) -1)
		 ((eq sgn '$zero) 0)
		 ((eq sgn '$pos) 1)

		 ;; multiplicative: signum(ab) = signum(a) * signum(b).
		 ((mtimesp x)
		  (muln (mapcar #'(lambda (s) (take '(%signum) s)) (margs x)) t))

		 ;; Reflection rule: signum(-x) --> -signum(x).
		 ((great (neg x) x) (neg (take '(%signum) (neg x))))
	
		 ;; nounform return
		 (t (eqtest (list '(%signum) x) e))))))))

;;; the following section is an attempt to
;;; over-wrap the $limit and $integrate programs to
;;; use these conventions.
;;; It would be nice if solve(1/x=inf,x) got x=0,
;;; but it doesn't so there is not much point in wrapping it.

;; fixing up limit,hackish.

(defvar limit-redefined nil)		; do this only once
(defvar integrate-redefined nil)	; do this only once

(unless limit-redefined 
  (setf (symbol-function '$old-limit)
    (symbol-function '$limit))
  (setf limit-redefined t))

(unless integrate-redefined  ;; nicer would be to store old def on proplist?
  (setf (symbol-function '$old-integrate)
    (symbol-function '$integrate))
  (setf (symbol-function 'old-simpinteg)
      (symbol-function 'simpinteg))
  (setf integrate-redefined t))

(defun $limit(&rest z)
  (let ((oldans (apply '$old-limit z))
	(newans nil))
;;    (format t "~%old-limit returns ~s" oldans)
    (setf newans
      (sublis
     '(
       ($inf . *ratinf)		;should cover minf too
       ($und . *ratind)
       ($infinity . *ratind)	; "complex infinity"?
       ($ind . *ratind)) ; bounded but unknown?? maybe interval(-oo,oo)
     oldans))
;;    (format t "~%newans is ~s" newans)
    (simplifya     newans     nil)))

(defun $integrate(&rest z)
  (let ((oldans (apply '$old-integrate z))
	(newans nil))
    (setf newans
	  (sublis
	   '(
	     ($inf . *ratinf)	;should cover minf too
	     ($und . *ratind)
	     ($infinity . *ratind)	; "complex infinity"?
	     ($ind . *ratind)) ; bounded but unknown?? maybe interval(-oo,oo)
	   oldans))
    ;;  (format t "~%newans is ~s" newans)
    (simplifya     newans     nil)))

(defun simpinteg(&rest z)
  (let ((oldans (apply 'old-simpinteg z))
	(newans nil))
    (setf newans
	  (sublis
	   '(
	     ($inf . *ratinf)	;should cover minf too
	     ($und . *ratind)
	     ($infinity . *ratind)	; "complex infinity"?
	     ($ind . *ratind)) ; bounded but unknown?? maybe interval(-oo,oo)
	   oldans))
    (simplifya     newans     nil)))

	      
;; from simp.lisp

(defun addk (x y)
  (cond ((eql x 0) y)
	((eql y 0) x)
	((and (numberp x) (numberp y))
         (cl-or-bfloat-binary-op mplus + x y))
	;; new, for infinities
	((ratinfp x) (if (alike1 x y) x (if (ratinfp y) *ratind x )))
	((ratinfp y) y)
	((or ($bfloatp x) ($bfloatp y)) ($bfloat (list '(mplus) x y)))
	
	(t (prog (g a b)
	      (cond ((numberp x)
		     (cond ((floatp x) (return (+ x (fpcofrat y))))
			   (t (setq x (list *ratsimpind x 1)))))
		    ((numberp y)
		     (cond ((floatp y) (return (+ y (fpcofrat x))))
			   (t (setq y (list *ratsimpind y 1))))))
	      (setq g (gcd (caddr x) (caddr y)))
	      (setq a (truncate (caddr x) g)
	            b (truncate (caddr y) g))
	      (return (timeskl (list *ratsimpind 1 g)
			       (list *ratsimpind
				     (+ (* (cadr x) b)
					   (* (cadr y) a))
				     (* a b))))))))









;;; rat (x/0) requires fix  to ratreduce and ratinvert from rat3b.lisp

(defmfun ratreduce (x y &aux b)
  (cond					;((pzerop y) (rat-error "quotient by zero"))
   ((pzerop y) 
    (if (pzerop x) '( 0 . 0) '(1 . 0))) ;; change! rjf 2/2016
   ((pzerop x) (rzero))
	((equal y 1) (cons x 1))
	((and $keepfloat (pcoefp y) (or $float (floatp y) (pfloatp x)))
	 (cons (pctimes (quotient 1.0 y) x) 1))
	(t (setq b (pgcdcofacts x y))
	   (setq b (ratalgdenom (rplacd (cdr b) (caddr b))))
	   (cond ((and modulus (pcoefp (cdr b)))
		  (cons (pctimes (crecip (cdr b)) (car b)) 1))
		 ((pminusp (cdr b))
		  (cons (pminus (car b)) (pminus (cdr b))))
		 (t b)))))


(defun ratinvert (y)
  (ratalgdenom
   (cond ;;((pzerop (car y)) (rat-error "`quotient' by `zero'"))
    ((pzerop (car y)) '(1 . 0)) ;; changed 2/2016 RJF
	 ((and modulus (pcoefp (car y)))
	  (cons (pctimes (crecip (car y)) (cdr y)) 1))
	 ((and $keepfloat (floatp (car y)))
	  (cons (pctimes (/ (car y)) (cdr y)) 1))
	 ((pminusp (car y)) (cons (pminus (cdr y)) (pminus (car y))))
	 (t (cons (cdr y) (car y))))))

#+ignore
(defmfun rattimes (x y gcdsw) ; no change,, yet
  (cond ($ratfac (facrtimes x y gcdsw))
	((and $algebraic gcdsw (ralgp x) (ralgp y))
	 (let ((w  (rattimes x y nil)))
	   (ratreduce (car w) (cdr w))))
	((equal 1 (cdr x))
	 (cond ((equal 1 (cdr y)) (cons (ptimes* (car x) (car y)) 1))
	       (t (cond (gcdsw (rattimes (ratreduce (car x) (cdr y))
					 (cons (car y) 1) nil))
			(t (cons (ptimes* (car x) (car y)) (cdr y)))))))
	((equal 1 (cdr y)) (rattimes y x gcdsw))
	(t (cond (gcdsw (rattimes (ratreduce (car x) (cdr y))
				  (ratreduce (car y) (cdr x)) nil))
		 (t (cons (ptimes* (car x) (car y))
			  (ptimes* (cdr x) (cdr y))))))))

;; must fix ratdisprep or nformat so  ((mrat ..) 0 . 0)
;; prints better .. not as 0  !!

;; from rat3e

(defun cdisrep (x &aux n d sign)
  (cond 
   ((pzerop (cdr x)) (if (pzerop (car x)) *ratind *ratinf))
   ((pzerop (car x)) 0)
	((or (equal 1 (cdr x)) (floatp (cdr x))) (pdisrep (car x)))
	(t (setq sign (cond ($ratexpand (setq n (pdisrep (car x))) 1)
			    ((pminusp (car x))
			     (setq n (pdisrep (pminus (car x)))) -1)
			    (t (setq n (pdisrep (car x))) 1)))
	   (setq d (pdisrep (cdr x)))
	   (cond ((and (numberp n) (numberp d))
		  (list *ratsimpind (* sign n) d))
		 ((and $ratdenomdivide $ratexpand
		       (not (atom n))
		       (eq (caar n) 'mplus))
		  (fancydis n d))
		 ((numberp d)
		  (list '(mtimes ratsimp)
			(list *ratsimpind sign d) n))
		 ((equal sign -1)
		  (cons '(mtimes ratsimp)
			(cond ((numberp n)
			       (list (* n -1)
				     (list '(mexpt ratsimp) d -1)))
			      (t (list sign n (list '(mexpt ratsimp) d -1))))))
		 ((equal n 1)
		  (list '(mexpt ratsimp) d -1))
		 (t (list '(mtimes ratsimp) n
			  (list '(mexpt ratsimp) d -1)))))))

;;;;fixed BUGS 4/13/16
#|  3.0^(1/0),  3.0^(0/0)  fixed
(0/0)^0  comes out indef  as does (0^0)/(0^0)

;; bug or feature:?? float(0/0) error
|#