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