share/to_poly_solve/to_poly.lisp

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