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 211Here Maxima solves for z. There are more than one solution for z, but Maxima returns 212just one solution. A user might think that there is one solution or that 213the equations are inconsistent. 214 215(2) eliminate([x+y-1,x+y-8],[x,y]) -> Argument to `last' is empty. 216 217Here the equations are inconsistent, but we get an error (from solve) instead 218of a clear message about what happened. 219 220(3) eliminate([a],[x]) -> Can't eliminate from only one equation -- an error. 221but eliminate([a,a],[x]) -> [a,a]. This is silly. 222 223(4) eliminate([x],[]) -> Can't eliminate from only one equation. 224but 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 227but 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 228a long time (I never let it run to completion). Unlike 'eliminate,' the function 'elim' 229makes some guesses about which polynomial to use as the pivot and which variable 230to 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 (let ((g)) 598 (cond 599 ((eq type '$integer) 600 (setq g (gentemp *integer-gentemp-prefix*)) 601 (setf (get g 'integer-gentemp) t) 602 (mfuncall '$declare g '$integer)) 603 604 ((eq type '$natural_number) 605 (setq g (gentemp *natural-gentemp-prefix*)) 606 (setf (get g 'natural-gentemp) t) 607 (mfuncall '$declare g '$integer) 608 (mfuncall '$assume (take '(mgeqp) g 0))) 609 610 ((eq type '$real) 611 (setq g (gentemp *real-gentemp-prefix*)) 612 (setf (get g 'real-gentemp) t)) 613 614 ((eq type '$complex) 615 (setq g (gentemp *complex-gentemp-prefix*)) 616 (setf (get g 'complex-gentemp) t) 617 (mfuncall '$declare g '$complex)) 618 619 (t 620 (setq g (gentemp *general-gentemp-prefix*)) 621 (setf (get g 'general-gentemp) t))) 622 623 g)) 624 625;; Find all the new-gentemp variables in an expression e and re-index them starting from 0. 626 627(defun $nicedummies (e) 628 (let (($listconstvars nil) 629 (z-vars nil) (z-k 0) 630 (n-vars nil) (n-k 0) 631 (r-vars nil) (r-k 0) 632 (c-vars nil) (c-k 0) 633 (g-vars nil) (g-k 0)) 634 635 (mapcar #'(lambda (s) (let ((a)) 636 (cond (($subvarp s)) 637 ((get s 'integer-gentemp) 638 (setq a (symbolconc *integer-gentemp-prefix* z-k)) 639 (mfuncall '$declare a '$integer) 640 (incf z-k) 641 (push `((mequal) ,s ,a) z-vars)) 642 643 ((get s 'natural-gentemp) 644 (setq a (symbolconc *natural-gentemp-prefix* n-k)) 645 (mfuncall '$declare a '$integer) 646 (mfuncall '$assume (take '(mgeqp) a 0)) 647 (incf n-k) 648 (push `((mequal) ,s ,a) n-vars)) 649 650 ((get s 'real-gentemp) 651 (setq a (symbolconc *real-gentemp-prefix* r-k)) 652 (incf r-k) 653 (push `((mequal) ,s ,a) r-vars)) 654 655 ((get s 'complex-gentemp) 656 (setq a (symbolconc *complex-gentemp-prefix* c-k)) 657 (mfuncall '$declare a '$complex) 658 (incf c-k) 659 (push `((mequal) ,s ,a) c-vars)) 660 661 ((get s 'general-gentemp) 662 (setq a (symbolconc *general-gentemp-prefix* g-k)) 663 (incf g-k) 664 (push `((mequal) ,s ,a) g-vars))))) 665 (margs ($sort ($listofvars e)))) 666 667 (push '(mlist) z-vars) 668 (push '(mlist) n-vars) 669 (push '(mlist) r-vars) 670 (push '(mlist) c-vars) 671 (push '(mlist) g-vars) 672 ($sublis g-vars ($sublis c-vars ($sublis r-vars ($sublis n-vars ($sublis z-vars e))))))) 673 674(defun $complex_number_p (e) 675 (complex-number-p e #'$numberp)) 676