share/to_poly_solve/to_poly.lisp

   1 ;;  Author Barton Willis
   2 ;;  University of Nebraska at Kearney
   3 ;;  Copyright (C) 2006, 2007, 2008, 2009 Barton Willis
   4
   5 ;;  This program is free software; you can redistribute it and/or modify
   6 ;;  it under the terms of the GNU General Public License as published by
   7 ;;  the Free Software Foundation; either version 2 of the License, or
   8 ;;  (at your option) any later version.
   9
  10 ;;  This program is distributed in the hope that it will be useful,
  11 ;;  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ;;  GNU General Public License for more details.
  14
  15
  16 ($put '$to_poly 2 '$version)
  17 ($load '$polynomialp)
  18
  19 (defmacro opapply (op args)
  20   `(simplify (cons (list ,op) ,args)))
  21
  22 ;; The next three functions convert max and min to abs functions.
  23
  24 (defun max-to-abs (e)
  25   (reduce #'(lambda (a b) (div (add (add a b) (take '(mabs) (sub a b))) 2)) e))
  26
  27 (defun min-to-abs (e)
  28   (reduce #'(lambda (a b) (div (sub (add a b) (take '(mabs) (sub a b))) 2)) e))
  29
  30 (defun convert-from-max-min-to-abs (e)
  31   (cond (($mapatom e) e)
  32         ((op-equalp e '$max) (max-to-abs (mapcar 'convert-from-max-min-to-abs (margs e))))
  33         ((op-equalp e '$min) (min-to-abs (mapcar 'convert-from-max-min-to-abs (margs e))))
  34         (t (opapply (mop e) (mapcar 'convert-from-max-min-to-abs (margs e))))))
  35
  36 (defun maxima-variable-p (e)
  37   (or (symbolp e) ($subvarp e)))
  38
  39 (defun list-subst (l p)
  40   (if (null l) p (list-subst (rest l) ($substitute (first l) p))))
  41
  42 ;; to_poly(p,vars) returns a polynomial in the variables 'vars' that has a zero whenever
  43 ;; p has a zero. When 1 is a member of vars, constant terms, such as sqrt(5) also get
  44 ;; converted to polynomial form. The value of vars defaults to all variables including
  45 ;; constants.
  46
  47
  48 (defun $to_poly (p &optional (vars 'convert-all-vars))
  49   (let (($listconstvars t) (q) (convert-cnst nil) (proviso nil) (non-alg) (qk) (pp nil) (subs))
  50
  51     (if (eq vars 'convert-all-vars) (setq vars ($cons 1 ($listofvars p))))
  52
  53     (if (not ($listp vars))
  54         (merror "The second argument to 'to_poly' must be a list"))
  55
  56     (cond (($member 1 vars)
  57            (setq convert-cnst t)
  58            (setq vars ($delete 1 vars))))
  59
  60     ;; If p is a list or a set, set p to the members in p; otherwise (list p)
  61     (setq p (if (or ($listp p) ($setp p)) (margs p) (list p)))
  62
  63     ;; Convert each member of p that is an equation to a nonequation.
  64     ;; Thus transform a = b into a - b.
  65     (setq p (mapcar 'meqhk p))
  66
  67     ;; Extract the deomominators of p and require them to not vanish.
  68     ;; Replace the list p by a list of the numerators.
  69     (setq q (mapcar '$ratdenom p))
  70     (setq p (mapcar '$ratnumer p))
  71
  72     (setq proviso (delete t (mapcar (lambda (s) (mnqp s 0)) q)))
  73     (setq proviso (mapcar #'(lambda (s) (maxima-substitute 'mnotequal '$notequal s)) proviso))
  74     (setq p (mapcar #'sratsimp p))
  75     ;;(multiple-value-setq (p g-vars) (non-algebraic-subst-list p vars))
  76     (setq p (non-algebraic-subst-list p vars))
  77     (setq non-alg ($second p))
  78     (setq p (margs ($first p)))
  79     ;; It's OK to send every expression through convert-from-max-min-to-abs.
  80     ;; I put in the conditional to skip the ratsimp for expressions that don't
  81     ;; involve max or min.
  82
  83     (setq pp nil)
  84     (setq subs nil)
  85     (dolist (pk p)
  86       (setq pk (if ($freeof '$max '$min pk) pk (sratsimp (convert-from-max-min-to-abs pk))))
  87       (setq pk (to-polynomial pk vars convert-cnst))
  88
  89       ;; After conversion, the members of pp can be rational expressions, not polynomials.
  90       ;; This happens, for example, when there is a 2 cos(x) --> exp(%i x) + 1/exp(%i x)
  91       ;; conversion. So we need to extract numerators of the members of p.
  92
  93       (setq proviso (append proviso (mapcar #'sratsimp (third pk))))
  94       (setq subs (append subs (mapcar #'sratsimp (second pk))))
  95       (setq pk (sratsimp (first pk)))
  96       (setq qk (sratsimp ($ratdenom pk)))
  97       (if (not (eq t (mnqp qk 0))) (push (take '(mnotequal) qk 0) proviso))
  98       (setq pk (sratsimp ($ratnumer pk)))
  99       (push pk pp))
 100
 101     (setq pp (append pp subs))
 102     (push '(mlist) pp)
 103     (push '(mlist) proviso)
 104     (list '(mlist) pp proviso non-alg)))
 105
 106 (defun to-polynomial (p vars convert-cnst)
 107   (let ((n) (b) (nv) (acc nil) (subs nil) (pk) (q) (inequal) (np-subs))
 108     (cond ((or (maxima-variable-p p)
 109                (mnump p)
 110                (and ($emptyp vars) (not convert-cnst))
 111                (and (not ($constantp p)) ($lfreeof vars p))
 112                (and ($constantp p) (not convert-cnst))
 113                (and (op-equalp p '$conjugate) (maxima-variable-p (second p))))
 114            (list p nil nil nil))
 115
 116           ((mexptp p)
 117
 118            (setq n (nth 2 p))
 119            (setq b (nth 1 p))
 120            (cond ((and (integerp n) (> n 0))
 121                   (list p nil nil nil))
 122
 123                  (($ratnump n)
 124                   (setq b (to-polynomial b vars convert-cnst))
 125                   (setq subs (second b))
 126                   (setq inequal (third b))
 127                   (setq np-subs (fourth b))
 128                   (setq b (first b))
 129                   (setq nv (new-gentemp 'general))
 130                   (cond ((or (mgrp n 0) (mnump b))
 131                          (let ((q) (r) (rr))
 132                            (setq q (take '($floor) n))
 133                            (setq r (sub n q))
 134                            (setq rr ($denom r))
 135                            (push (take '(mleqp) (take '($parg) nv) (div '$%pi rr)) inequal)
 136                            (push (take '(mlessp) (div '$%pi (neg rr)) (take '($parg) nv)) inequal)
 137                            (push (take '(mequal) (power b ($num r)) (power nv ($denom r))) subs)
 138                            (push (take '(mequal) p nv) np-subs)
 139                            (list (mul nv (power b q))  subs inequal np-subs)))
 140
 141                         ;   (push (take '(mleqp) (take '($parg) nv) (mul n '$%pi)) inequal)
 142                         ;   (push (take '(mlessp) (mul -1 '$%pi n) (take '($parg) nv)) inequal)
 143                         ;   (push (take '(mequal) (power b ($num n)) (power nv ($denom n))) subs)
 144                         ;   (push (take '(mequal) p nv) np-subs)
 145                         ;   (list nv subs inequal np-subs)))
 146
 147                         (t
 148                          (setq n (neg n))
 149                          (push (take '(mequal) 1 (mul (power nv ($denom n)) (power b ($num n)))) subs)
 150                          (push (take '(mlessp) (take '($parg) nv) (mul n '$%pi)) inequal)
 151                          (push (take '(mleqp) (mul -1 '$%pi n) (take '($parg) nv)) inequal)
 152                          (push (take '(mnotequal) nv 0) inequal)
 153                          (list nv subs inequal np-subs))))
 154
 155                  (t (merror "Nonalgebraic argument given to 'to_poly'"))))
 156
 157           ((op-equalp p 'mabs)
 158            (setq b (to-polynomial (first (margs p)) vars convert-cnst))
 159            (setq acc (second b))
 160            (setq inequal (third b))
 161            (setq np-subs (fourth b))
 162            (setq b (first b))
 163            (setq nv (new-gentemp '$general))
 164            ;; Ok, I'm indecisive.
 165            ;(push (take '(mgeqp) nv 0) inequal)
 166            ;(push (take '(mequal) (take '($parg) nv) 0) inequal)
 167            (push (take '($isnonnegative_p) nv) inequal)
 168            (list nv (cons (take '(mequal) (mul b (take '($conjugate) b)) (mul nv (take '($conjugate) nv))) acc)
 169                  inequal np-subs))
 170
 171           ((mtimesp p)
 172            (setq acc 1)
 173            (setq p (margs p))
 174            (while p
 175              (setq pk (first p))
 176              (setq p (rest p))
 177              (setq q (to-polynomial pk vars convert-cnst))
 178              (setq acc (mul acc (first q)))
 179              (setq subs (append (second q) subs))
 180              (setq inequal (append inequal (third q)))
 181              (setq np-subs (append np-subs (fourth q)))
 182              (setq vars ($append vars ($listofvars `((mlist) ,@subs))))
 183
 184              (setq p (mapcar #'(lambda (s) (list-subst np-subs s)) p)))
 185            (list acc subs inequal np-subs))
 186
 187           ((mplusp p)
 188            (setq acc 0)
 189            (setq p (margs p))
 190            (while p
 191              (setq pk (first p))
 192              (setq p (rest p))
 193              (setq q (to-polynomial pk vars convert-cnst))
 194              (setq acc (add acc (first q)))
 195              (setq subs (append (second q) subs))
 196              (setq inequal (append (third q) inequal))
 197              (setq np-subs (append (fourth q) np-subs))
 198              (setq vars ($append vars ($listofvars `((mlist) ,@subs))))
 199              (setq p (mapcar #'(lambda (s) (list-subst np-subs s)) p)))
 200            (list acc subs inequal np-subs))
 201
 202
 203           (t (merror "Nonalgebraic argument given to 'to_poly'")))))
 204
 205
 206 #|
 207   Things I don't like about eliminate:
 208
 209 (1)  eliminate([x + x*y + z-4, x+y+z-12, x^2 + y^2 + z^2=7],[x,y,z]) -> [4]
 210
 211 Here Maxima solves for z. There are more than one solution for z, but Maxima returns
 212 just one solution. A user might think that there is one solution or that
 213 the equations are inconsistent.
 214
 215 (2)  eliminate([x+y-1,x+y-8],[x,y]) -> Argument to `last' is empty.
 216
 217 Here the equations are inconsistent, but we get an error (from solve) instead
 218 of a clear message about what happened.
 219
 220 (3) eliminate([a],[x]) -> Can't eliminate from only one equation -- an error.
 221 but eliminate([a,a],[x]) -> [a,a]. This is silly.
 222
 223 (4) eliminate([x],[]) -> Can't eliminate from only one equation.
 224 but eliminate([x,x],[]) -> [x,x]. Again, this is silly.
 225
 226 (5) elim([x^3-y^2,x^7 + y+z*p,q*q-23+ x+y,p-q,x+y+z-5],[x,y,z,p]) takes 0.3 second
 227 but eliminate([x^3-y^2,x^7 + y+z*p,q*q-23+ x+y,p-q,x+y+z-5],[x,y,z,p]) takes
 228 a long time (I never let it run to completion). Unlike 'eliminate,' the function 'elim'
 229 makes some guesses about which polynomial to use as the pivot and which variable
 230 to eliminate.
 231 |#
 232
 233 (defun require-maxima-variable (x context-string)
 234   (setq x (ratdisrep x))
 235   (if (maxima-variable-p x) x
 236     (merror "Function ~:M expects a symbol, instead found ~:M" context-string x)))
 237
 238 ;; Simplify a polynomial equation p = 0 by
 239
 240 ;;   (1) nonzero constant * p --> p,
 241 ;;   (2) p^n --> p.
 242
 243 ;; If you want to use $factor instead of $sqfr, go ahead. But if you do that, you might want to
 244 ;; locally set factorflag to false.
 245
 246 (defun suppress-multiple-zeros (q)
 247   (setq q ($sqfr q))
 248   (cond ((mnump q) (if (zerop1 q) 0 1))
 249         ((mtimesp q) (muln (mapcar 'suppress-multiple-zeros (margs q)) t))
 250         ((and (mexptp q) (integerp (third q)) (> (third q) 0)) (second q))
 251         (($constantp q) 1) ; takes care of things like (1 + sqrt(5)) * x --> x.
 252         (t q)))
 253
 254 ;; Using eq as the "pivot," eliminate x from the list or set of equations eqs.
 255
 256 (defun $eliminate_using (eqs eq x)
 257   (if (or ($listp eqs) ($setp eqs))
 258       (progn
 259         (setq eqs (mapcar #'(lambda (s) ($ratexpand (meqhk s))) (margs eqs)))
 260         (setq eqs (margs (opapply '$set eqs))))
 261     (merror "The first argument to 'eliminate_using' must be a list or a set"))
 262
 263   (setq x (require-maxima-variable x "$eliminate_using"))
 264
 265   (if (not (every #'(lambda (s) ($polynomialp s `((mlist) ,x)
 266                                               `((lambda) ((mlist) s) (($freeof) ,x s)))) eqs))
 267       (merror "The first argument to 'eliminate_using' must be a set or list of polynomials"))
 268
 269   (setq eq ($ratexpand (meqhk eq)))
 270   (if (not ($polynomialp eq `((mlist) ,x) `((lambda) ((mlist) s) (($freeof) ,x s))))
 271       (merror "The second argument to 'eliminate_using' must be a polynomial"))
 272
 273   (setq eqs (mapcar #'suppress-multiple-zeros eqs))
 274   (setq eq (suppress-multiple-zeros eq))
 275   (opapply '$set (eliminate-using eqs eq x)))
 276
 277 (defun eliminate-using (eqs eq x)
 278   (delete 0 (mapcar #'(lambda (s) (suppress-multiple-zeros ($resultant s eq x))) eqs)))
 279
 280 ;; Return an upper bound for the total degree of the polynomial p in the variables vars.
 281 ;; When p is fully expanded, the bound is exact.
 282
 283 (defun degree-upper-bound (p vars)
 284   (cond ((maxima-variable-p p) (if (member p vars :test #'equal) 1 0))
 285
 286         ((mnump p) 0)
 287
 288         ((and (mexptp p) (integerp (third p)) (> (third p) 0))
 289          (* (degree-upper-bound (nth 1 p) vars) (nth 2 p)))
 290
 291         ((mtimesp p)
 292          (reduce '+ (mapcar #'(lambda (s) (degree-upper-bound s vars)) (margs p))))
 293
 294         ((mplusp p)
 295          (simplify `(($max) ,@(mapcar #'(lambda (s) (degree-upper-bound s vars)) (margs p)))))
 296
 297         ((apply '$freeof (append vars (list p))) 0)
 298         (t (merror "Nonpolynomial argument given to degree-upper-bound"))))
 299
 300 (defun unk-eliminate (eqs vars &optional (pivots nil))
 301    (let ((ni) (n-min nil) (xeqs `(($set))) (pivot-var) (pivot-eq) (acc `(($set))) ($ratfac nil))
 302      (cond ((or (null vars) (null eqs)) (list eqs pivots))
 303            (t
 304
 305             ;; The pivot is a nonconstant member of eqs with minimal total degree.
 306             ;; The  constant members of eqs get adjoined into acc -- all other  members get
 307             ;; adjoined into xeqs. Since each member of eqs has been ratexpanded,
 308             ;; the degree bound is exact.
 309
 310             (dolist (ei eqs)
 311               (setq ei ($ratexpand (suppress-multiple-zeros ei)))
 312               (setq ni (degree-upper-bound ei vars))
 313
 314               (if (and (or (null n-min) (< ni n-min)) (> ni 0))
 315                   (setq n-min ni pivot-eq ei))
 316
 317               (if (and (> ni 0) (not (equal 0 ei))) (setq xeqs ($adjoin ei xeqs)) (setq acc ($adjoin ei acc))))
 318
 319             (setq xeqs (margs xeqs))
 320             (setq acc (margs acc))
 321             ;; Now we'll decide which variable to eliminate. The pivot variable
 322             ;; is the variable that has the least (but nonzero) degree in pivot-eq.
 323
 324             (setq n-min nil)
 325             (dolist (vi vars)
 326               (setq ni (degree-upper-bound pivot-eq (list vi)))
 327               (if (and (or (null n-min) (< ni n-min)) (> ni 0))
 328                   (setq pivot-var vi n-min ni)))
 329
 330             (if (null pivot-var) (list eqs pivots)
 331               (unk-eliminate (append acc (eliminate-using xeqs pivot-eq pivot-var)) (delete pivot-var vars)
 332                              (cons pivot-eq pivots)))))))
 333
 334 (defun $elim (eqs x)
 335   (if (or ($listp eqs) ($setp eqs))
 336       (progn
 337         (setq eqs (mapcar #'(lambda (s) ($ratexpand (suppress-multiple-zeros (meqhk s)))) (margs eqs)))
 338         (setq eqs (margs (opapply '$set eqs))))
 339     (merror "The first argument to 'elim' must be a list or a set"))
 340
 341   (setq x (margs (cond (($listp x) ($setify x))
 342                        (($setp x) x)
 343                        (t (merror "The second argument to 'elim' must be a list or a set")))))
 344
 345   (setq x (mapcar #'(lambda (s) (require-maxima-variable s "$elim")) x))
 346
 347   (setq x (opapply 'mlist x))
 348   (if (not (every #'(lambda (s) ($polynomialp s x `((lambda) ((mlist) s) (($lfreeof) ,x s)))) eqs))
 349       (merror "Each member of the first argument to 'elim' must be a polynomial"))
 350
 351   (setq x (margs x))
 352   (opapply 'mlist (mapcar #'(lambda (s) (opapply 'mlist s)) (unk-eliminate eqs x))))
 353
 354 (defun $elim_allbut (eqs x)
 355   (let (($listconstvars nil) (v))
 356     (setq v ($listofvars eqs))
 357     (setq x (cond (($listp x) ($setify x))
 358                   (($setp x) x)
 359                   (t (take '($set) x))))
 360     (setq v ($setdifference ($setify v) x))
 361     ($elim eqs ($listify v))))
 362
 363 ;; Return a list of the arguments to the operator op in the expression e.
 364 ;; This function only looks inside + and * and it doesn't recurse inside
 365 ;; + and *. So
 366
 367 ;;  gather-args(log(x) + log(a + log(b)), log) --> (x, a + log(b)),
 368 ;;  gather-args(cos(log(x)), log) --> (),
 369 ;;  gather-args(x^log(x) + log(s)^x, log) --> ().
 370
 371 (defun gather-args (e op vars)
 372   (cond ((and (consp e) (consp (car e)) (eq (caar e) 'mlist))
 373          (apply 'append (mapcar #'(lambda (s) (gather-args s op vars)) (margs e))))
 374         (($mapatom e) nil)
 375         ((and (eq (mop e) op) (not ($lfreeof vars e))) (margs e))
 376         ((memq (mop e) '(mplus mtimes))
 377          (let ((acc nil))
 378            (setq e (margs e))
 379            (dolist (ek e acc)
 380              (setq acc (append acc (gather-args ek op vars))))))
 381         (t nil)))
 382
 383 (defun gather-nonrational-powers (e vars)
 384   (cond ((and (consp e) (consp (car e)) (eq (caar e) 'mlist))
 385          (apply 'append (mapcar #'(lambda (s) (gather-nonrational-powers s vars)) (margs e))))
 386         (($mapatom e) nil)
 387         ((and (eq (mop e) 'mexpt) (memq (second e) vars) (not ($ratnump (third e)))) (list e))
 388         ((memq (mop e) '(mplus mtimes))
 389          (mapcan #'(lambda (s) (gather-nonrational-powers s vars)) (margs e)))
 390         (t nil)))
 391
 392 (defun gather-exp-args (e vars)
 393   (cond (($mapatom e) nil)
 394         ((eq (mop e) 'mexpt)
 395          (if (and (eq (second e) '$%e) (not ($lfreeof vars (third e)))) (list (third e))
 396            (gather-exp-args (second e) vars)))
 397         (t (mapcan #'(lambda (s) (gather-exp-args s vars)) (margs e)))))
 398
 399 ;; Return true iff the set {a,b} is linearly dependent. Thus return true iff
 400 ;; either a = 0 or b = 0, or a/b is a number.
 401
 402 (defun $linearly_dependent_p (a b)
 403   (linearly-dependent-p a b))
 404
 405 (defun linearly-dependent-p (a b)
 406   (if (or (=0 a)
 407           (=0 b)
 408           (mnump (sratsimp (div a b)))) t nil))
 409
 410 (defun exponentialize-if (e vars)
 411   (cond (($mapatom e) e)
 412         ((and (trigp (mop e)) (not ($lfreeof vars e))) ($exponentialize e))
 413         (t (opapply (mop e) (mapcar #'(lambda (s) (exponentialize-if s vars)) (margs e))))))
 414
 415 (defun logarc-if (e vars)
 416   (cond (($mapatom e) e)
 417         ((and (arcp (mop e)) (not ($lfreeof vars e))) ($logarc e))
 418         (t (opapply (mop e) (mapcar #'(lambda (s) (logarc-if s vars)) (margs e))))))
 419
 420 (defun non-algebraic-subst (e vars)
 421   (let ((log-args nil) (exp-args nil) (mexpt-args) (ee) (s) (g) (p) (q) (na-subs nil) (g-vars nil) (sz))
 422
 423     (setq e ($ratexpand (exponentialize-if (logarc-if e vars) vars)))
 424     (setq log-args (gather-args e '%log vars))
 425     (setq log-args (margs (opapply '$set log-args)))
 426     (setq exp-args (gather-exp-args e vars))
 427     (setq exp-args (opapply '$set exp-args))
 428     (setq s (margs (simplify ($equiv_classes exp-args 'linearly-dependent-p))))
 429
 430     (dolist (sk s)
 431       (setq g (new-gentemp '$general))
 432       (push g g-vars)
 433       (setq sz ($first (apply '$ezgcd (margs sk))))
 434       (setq p (power '$%e sz))
 435       (push (take '(mequal) g p) na-subs)
 436       (setq sk (margs sk))
 437       (dolist (sl sk)
 438         (setq q (sratsimp (div sl sz)))
 439         (setq e ($ratexpand ($substitute (power g  q) (power '$%e sl) e)))))
 440
 441     (dolist (sk log-args)
 442       (setq g (new-gentemp '$general))
 443       (push g g-vars)
 444       (setq p (take '(%log) sk))
 445       (setq e ($substitute g p e))
 446       (push (take '(mequal) g p) na-subs))
 447
 448     ;; Attempt the substitution %g = x^a, where x is in vars and a is not a rational number.
 449     ;; Accept the substitution if it eliminates x from e and it changes e, otherwise reject it.
 450
 451     (setq mexpt-args (gather-nonrational-powers e (margs vars)))
 452     (dolist (sk mexpt-args)
 453       (setq g (new-gentemp '$general))
 454       (setq ee ($ratsubst g sk e))
 455       (if (and (freeof (second sk) ee) (not (like e ee)))
 456           (progn
 457             (push g g-vars)
 458             (setq e ee)
 459             (push (take '(mequal) g sk) na-subs))))
 460
 461     (values `((mlist) ,e ((mlist) ,@na-subs)) g-vars)))
 462
 463 (defun non-algebraic-subst-list (l vars)
 464   (let ((log-args nil) (exp-args nil) (mexpt-args) (ll) (s) (g) (p) (q) (na-subs nil) (g-vars nil) (sz))
 465
 466     (setq l (mapcar #'(lambda (s) ($ratexpand (exponentialize-if (logarc-if s vars) vars))) l))
 467     (push '(mlist) l)
 468     (setq log-args (gather-args l '%log vars))
 469     (setq log-args (margs (opapply '$set log-args)))
 470
 471     (setq exp-args (gather-exp-args l vars))
 472     (setq exp-args (opapply '$set exp-args))
 473     (setq s (margs (simplify ($equiv_classes exp-args 'linearly-dependent-p))))
 474
 475     (dolist (sk s)
 476       (setq g (new-gentemp '$general))
 477       (push g g-vars)
 478       (setq sz ($first (apply '$ezgcd (margs sk))))
 479       (setq p (power '$%e sz))
 480       (push (take '(mequal) g p) na-subs)
 481       (setq sk (margs sk))
 482       (dolist (sl sk)
 483         (setq q (sratsimp (div sl sz)))
 484         (setq l ($ratexpand ($substitute (power g  q) (power '$%e sl) l)))))
 485
 486     (dolist (sk log-args)
 487       (setq g (new-gentemp '$general))
 488       (push g g-vars)
 489       (setq p (take '(%log) sk))
 490       (setq l ($substitute g p l))
 491       (push (take '(mequal) g p) na-subs))
 492
 493     ;; Attempt the substitution %g = x^a, where x is in vars and a is not a rational number.
 494     ;; Accept the substitution if it eliminates x from e and it changes e, otherwise reject it.
 495
 496     (setq mexpt-args (gather-nonrational-powers l (margs vars)))
 497     (dolist (sk mexpt-args)
 498       (setq g (new-gentemp '$general))
 499       (setq ll ($ratsubst g sk l))
 500       (if (and (freeof (second sk) ll) (not (like l ll)))
 501           (progn
 502             (push g g-vars)
 503             (setq l ll)
 504             (push (take '(mequal) g sk) na-subs))))
 505
 506     (values `((mlist) ,l ((mlist) ,@na-subs)) g-vars)))
 507
 508 ;; A simplifying carg function.
 509
 510 (setf (get '$parg 'operators) 'simp-parg)
 511
 512 (defun simp-parg (e yy z)
 513   (declare (ignore yy))
 514   (let (($domain '$complex) (sgn) (isgn))
 515
 516     (oneargcheck e)
 517     (setq e (simplifya (specrepcheck (cadr e)) z))
 518     (setq e (sratsimp ($exponentialize e)))
 519
 520     (cond ((zerop1 e) e) ;; parg(0) = 0,parg(0.0) = 0.0, parg(0.0b0) = 0.0b0.
 521
 522           ;; For a complex number, use atan2
 523           ((complex-number-p e '$constantp)
 524            (take '($atan2) ($imagpart e) ($realpart e)))
 525
 526           ;; off the negative real axis, parg(conjugate(x)) = -parg(x)
 527           ((and (op-equalp e '$conjugate) (off-negative-real-axisp e))
 528            (neg (take '($parg) (cadr e))))
 529
 530           ;; parg(a * x) = parg(x) provided a > 0.
 531           ((and (mtimesp e) (eq t (mgrp (cadr e) 0)))
 532            (take '($parg) (reduce 'mul (cddr e))))
 533
 534           ;; parg exp(a + %i*b) = %pi-mod(%pi-b,2*%pi)
 535           ((and (mexptp e) (eq '$%e (second e)) (linearp (third e) '$%i))
 536            (sub '$%pi (take '($mod) (sub '$%pi ($imagpart (third e))) (mul 2 '$%pi))))
 537
 538           ;; parg(x^number) = number * parg(x), provided number in (-1,1).
 539           ((and (mexptp e) ($numberp (caddr e)) (eq t (mgrp (caddr e) -1)) (eq t (mgrp 1 (caddr e))))
 540            (mul (caddr e) (take '($parg) (cadr e))))
 541
 542           ;; sign rules parg(x) = %pi, x < 0 and parg(x) = 0, x > 0
 543           ((eq '$neg (setq sgn ($csign e))) '$%pi)
 544           ((eq '$pos sgn) 0)
 545
 546           ;; more sign rules parg(%i * x) = %pi /2, x > 0 and -%pi / 2, x < 0.
 547           ((and (eq '$imaginary sgn) (setq isgn ($csign (div e '$%i))) (eq '$pos isgn))
 548            (div '$%pi 2))
 549
 550           ((and (eq '$imaginary sgn) (eq '$neg isgn))
 551            (div '$%pi -2))
 552
 553           ;; nounform return
 554           (t `(($parg simp) ,e)))))
 555
 556 (setf (get '$isreal_p 'operators) 'simp-isreal-p)
 557
 558 (defun simp-isreal-p (e yy z)
 559   (declare (ignore yy))
 560   (oneargcheck e)
 561   (let ((ec) ($domain '$complex))
 562
 563     (setq e (simplifya (specrepcheck (cadr e)) z))
 564
 565     ;; Simplifications:
 566
 567     ;;  (1) if r is real then isreal_p(r + x) --> isreal_p(x),
 568     ;;  (2) if r is real then isreal_p(r * x) --> isreal_p(x),
 569     ;;  (3) if e = conjugate(e) then true,
 570     ;;  (4) if is(notequal(e,conjugate(e))) then false.
 571
 572     (cond ((mtimesp e)
 573            (setq e (muln (mapcar #'(lambda (s) (if (eq t (take '($isreal_p) s)) 1 s)) (margs e)) t)))
 574           ((mplusp e)
 575            (setq e (addn (mapcar #'(lambda (s) (if (eq t (take '($isreal_p) s)) 0 s)) (margs e)) t))))
 576
 577     ;; If e is constant, apply rectform.
 578     (if ($constantp e) (setq e ($rectform e)))
 579     (setq ec (take '($conjugate) e))
 580     (cond ((eq t (meqp e ec)) t)
 581           ((eq t (mnqp e ec)) nil)
 582           (t `(($isreal_p simp) ,e)))))
 583
 584 (defvar *integer-gentemp-prefix* "$%Z")
 585 (defvar *natural-gentemp-prefix* "$%N")
 586 (defvar *real-gentemp-prefix*    "$%R")
 587 (defvar *complex-gentemp-prefix* "$%C")
 588 (defvar *general-gentemp-prefix* "$%G")
 589
 590 ;; Return a new gentemp with a prefix that is determined by 'type.' Also put
 591 ;; an identifying tag on the symbol property of each new gentemp.
 592
 593 (defun $new_variable (type)
 594   (new-gentemp type))
 595
 596 (defun new-gentemp (type)
 597   (flet ((%gentemp (prefix)
 598            (intern (symbol-name (gensym prefix)) :maxima)))
 599     (let ((g))
 600       (cond
 601        ((eq type '$integer)
 602         (setq g (%gentemp *integer-gentemp-prefix*))
 603         (setf (get g 'integer-gentemp) t)
 604         (mfuncall '$declare g '$integer))
 605
 606        ((eq type '$natural_number)
 607         (setq g (%gentemp *natural-gentemp-prefix*))
 608         (setf (get g 'natural-gentemp) t)
 609         (mfuncall '$declare g '$integer)
 610         (mfuncall '$assume (take '(mgeqp) g 0)))
 611
 612        ((eq type '$real)
 613         (setq g (%gentemp *real-gentemp-prefix*))
 614         (setf (get g 'real-gentemp) t))
 615
 616        ((eq type '$complex)
 617         (setq g (%gentemp *complex-gentemp-prefix*))
 618         (setf (get g 'complex-gentemp) t)
 619         (mfuncall '$declare g '$complex))
 620
 621        (t
 622         (setq g (%gentemp *general-gentemp-prefix*))
 623         (setf (get g 'general-gentemp) t)))
 624
 625       g)))
 626
 627 ;; Find all the new-gentemp variables in an expression e and re-index them starting from 0.
 628
 629 (defun $nicedummies (e)
 630   (let (($listconstvars nil)
 631         (z-vars nil) (z-k 0)
 632         (n-vars nil) (n-k 0)
 633         (r-vars nil) (r-k 0)
 634         (c-vars nil) (c-k 0)
 635         (g-vars nil) (g-k 0))
 636
 637     (mapcar #'(lambda (s) (let ((a))
 638                             (cond (($subvarp s))
 639                                   ((get s 'integer-gentemp)
 640                                    (setq a (symbolconc *integer-gentemp-prefix* z-k))
 641                                    (mfuncall '$declare a '$integer)
 642                                    (incf z-k)
 643                                    (push `((mequal) ,s ,a) z-vars))
 644
 645                                   ((get s 'natural-gentemp)
 646                                    (setq a (symbolconc *natural-gentemp-prefix* n-k))
 647                                    (mfuncall '$declare a '$integer)
 648                                    (mfuncall '$assume (take '(mgeqp) a 0))
 649                                    (incf n-k)
 650                                    (push `((mequal) ,s ,a) n-vars))
 651
 652                                   ((get s 'real-gentemp)
 653                                    (setq a (symbolconc *real-gentemp-prefix* r-k))
 654                                    (incf r-k)
 655                                    (push `((mequal) ,s ,a) r-vars))
 656
 657                                   ((get s 'complex-gentemp)
 658                                    (setq a (symbolconc *complex-gentemp-prefix* c-k))
 659                                    (mfuncall '$declare a '$complex)
 660                                    (incf c-k)
 661                                    (push `((mequal) ,s ,a) c-vars))
 662
 663                                   ((get s 'general-gentemp)
 664                                    (setq a (symbolconc *general-gentemp-prefix* g-k))
 665                                    (incf g-k)
 666                                    (push `((mequal) ,s ,a) g-vars)))))
 667             (margs ($sort ($listofvars e))))
 668
 669     (push '(mlist) z-vars)
 670     (push '(mlist) n-vars)
 671     (push '(mlist) r-vars)
 672     (push '(mlist) c-vars)
 673     (push '(mlist) g-vars)
 674     ($sublis g-vars ($sublis c-vars ($sublis r-vars ($sublis n-vars ($sublis z-vars e)))))))
 675
 676 (defun $complex_number_p (e)
 677   (complex-number-p e #'$numberp))