src/trigi.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter, University of Texas  ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module trigi)
  14
  15 (load-macsyma-macros mrgmac)
  16
  17 ;;; Arithmetic utilities.
  18
  19 (defun sqrt1-x^2 (x)
  20   (power (sub 1 (power x 2)) 1//2))
  21
  22 (defun sqrt1+x^2 (x)
  23   (power (add 1 (power x 2)) 1//2))
  24
  25 (defun sqrtx^2-1 (x)
  26   (power (add (power x 2) -1) 1//2))
  27
  28 (defun sq-sumsq (x y)
  29   (power (add (power x 2) (power y 2)) 1//2))
  30
  31 (defun trigp (func)
  32   (member func '(%sin %cos %tan %csc %sec %cot %sinh %cosh %tanh %csch %sech %coth)
  33           :test #'eq))
  34
  35 (defun arcp (func)
  36   (member func '(%asin %acos %atan %acsc %asec %acot %asinh %acosh %atanh %acsch %asech %acoth)
  37           :test #'eq))
  38
  39 ;;; The trigonometric functions distribute of lists, matrices and equations.
  40
  41 (dolist (x '(%sin   %cos   %tan   %cot   %csc   %sec
  42              %sinh  %cosh  %tanh  %coth  %csch  %sech
  43              %asin  %acos  %atan  %acot  %acsc  %asec
  44              %asinh %acosh %atanh %acoth %acsch %asech))
  45   (setf (get x 'distribute_over) '(mlist $matrix mequal)))
  46
  47 (defun domain-error (x f)
  48   (merror (intl:gettext "~A: argument ~:M isn't in the domain of ~A.") f (complexify x) f))
  49
  50 ;; Some Lisp implementations goof up branch cuts for ASIN, ACOS, and/or ATANH.
  51 ;; Here are definitions which have the right branch cuts
  52 ;; (assuming LOG, PHASE, and SQRT have the right branch cuts).
  53 ;; Don't bother trying to sort out which implementations get it right or wrong;
  54 ;; we'll make all implementations use these functions.
  55
  56 ;; Apply formula from CLHS if X falls on a branch cut.
  57 ;; Otherwise punt to CL:ASIN.
  58 (defun maxima-branch-asin (x)
  59   ;; Test for (IMAGPART X) is EQUAL because signed zero is EQUAL to zero.
  60   (if (and (> (abs (realpart x)) 1.0) (equal (imagpart x) 0.0))
  61       ;; The formula from CLHS is asin(x) = -%i*log(%i*x+sqrt(1-x^2)).
  62       ;; This has problems with overflow for large x.
  63       ;;
  64       ;; Let's rewrite it, where abs(x)>1
  65       ;;
  66       ;; asin(x) = -%i*log(%i*x+abs(x)*sqrt(1-1/x^2))
  67       ;;         = -%i*log(%i*x*(1+abs(x)/x*sqrt(1-1/x^2)))
  68       ;;         = -%i*[log(abs(x)*abs(1+abs(x)/x*sqrt(1-1/x^2)))
  69       ;;                 + %i*arg(%i*x*(1+abs(x)/x*sqrt(1-1/x^2)))]
  70       ;;         = -%i*[log(abs(x)*(1+abs(x)/x*sqrt(1-1/x^2)))
  71       ;;                 + %i*%pi/2*sign(x)]
  72       ;;         = %pi/2*sign(x) - %i*[log(abs(x)*(1+abs(x)/x*sqrt(1-1/x^2))]
  73       ;;
  74       ;; Now, look at log part.  If x > 0, we have
  75       ;;
  76       ;;    log(x*(1+sqrt(1-1/x^2)))
  77       ;;
  78       ;; which is just fine.  For x < 0, we have
  79       ;;
  80       ;;    log(abs(x)*(1-sqrt(1-1/x^2))).
  81       ;;
  82       ;; But
  83       ;;    1-sqrt(1-1/x^2) = (1-sqrt(1-1/x^2))*(1+sqrt(1-1/x^2))/(1+sqrt(1-1/x^2))
  84       ;;                    = (1-(1-1/x^2))/(1+sqrt(1-1/x^2))
  85       ;;                    = 1/x^2/(1+sqrt(1-1/x^2))
  86       ;;
  87       ;; So
  88       ;;
  89       ;;    log(abs(x)*(1-sqrt(1-1/x^2)))
  90       ;;        = log(abs(x)/x^2/(1+sqrt(1-1/x^2)))
  91       ;;        = -log(x^2/abs(x)*(1+sqrt(1-1/x^2))
  92       ;;        = -log(abs(x)*(1+sqrt(1-1/x^2)))
  93       ;;
  94       ;; Thus, for x < 0,
  95       ;;
  96       ;; asin(x) = -%pi/2+%i*log(abs(x)*(1+sqrt(1-1/x^2)))
  97       ;;         = -asin(-x)
  98       ;;
  99       ;; If we had an accurate f(x) = log(1+x) function, we should
 100       ;; probably evaluate log(1+sqrt(1-1/x^2)) via f(x) instead of
 101       ;; log.  One other accuracy change is to evaluate sqrt(1-1/x^2)
 102       ;; as sqrt(1-1/x)*sqrt(1+1/x), because 1/x^2 won't underflow as
 103       ;; soon as 1/x.
 104       (let* ((absx (abs x))
 105              (recip (/ absx))
 106              (result (complex (/ #.(float pi) 2)
 107                               (- (log (* absx
 108                                          (1+ (* (sqrt (+ 1 recip))
 109                                                 (sqrt (- 1 recip))))))))))
 110         (if (minusp x)
 111             (- result)
 112             result))
 113       (cl:asin x)))
 114
 115 ;; Apply formula from CLHS if X falls on a branch cut.
 116 ;; Otherwise punt to CL:ACOS.
 117 (defun maxima-branch-acos (x)
 118   ; Test for (IMAGPART X) is EQUAL because signed zero is EQUAL to zero.
 119   (if (and (> (abs (realpart x)) 1.0) (equal (imagpart x) 0.0))
 120     (- #.(/ (float pi) 2) (maxima-branch-asin x))
 121     (cl:acos x)))
 122
 123 (defun maxima-branch-acot (x)
 124   ;; Allow 0.0 in domain of acot, otherwise use atan(1/x)
 125   (if (and (equal (realpart x) 0.0) (equal (imagpart x) 0.0))
 126     #.(/ (float pi) 2)
 127     (cl:atan (/ 1 x))))
 128
 129 ;; Apply formula from CLHS if X falls on a branch cut.
 130 ;; Otherwise punt to CL:ATANH.
 131 (defun maxima-branch-atanh (x)
 132   ; Test for (IMAGPART X) is EQUAL because signed zero is EQUAL to zero.
 133   (if (and (> (abs (realpart x)) 1.0) (equal (imagpart x) 0.0))
 134     (/ (- (cl:log (+ 1 x)) (cl:log (- 1 x))) 2)
 135     (cl:atanh x)))
 136
 137 ;; Fill the hash table.
 138 (macrolet ((frob (mfun dfun) `(setf (gethash ',mfun *flonum-op*) ,dfun)))
 139   (frob mplus #'+)
 140   (frob mtimes #'*)
 141   (frob mquotient #'/)
 142   (frob mminus #'-)
 143
 144   (frob %cos #'cl:cos)
 145   (frob %sin #'cl:sin)
 146   (frob %tan #'cl:tan)
 147
 148   (frob %sec #'(lambda (x)
 149                  (let ((y (ignore-errors (/ 1 (cl:cos x)))))
 150                    (if y y (domain-error x 'sec)))))
 151
 152   (frob %csc #'(lambda (x)
 153                  (let ((y (ignore-errors (/ 1 (cl:sin x)))))
 154                    (if y y (domain-error x 'csc)))))
 155
 156   (frob %cot #'(lambda (x)
 157                  (let ((y (ignore-errors (/ 1 (cl:tan x)))))
 158                    (if y y (domain-error x 'cot)))))
 159
 160   (frob %acos #'maxima-branch-acos)
 161   (frob %asin #'maxima-branch-asin)
 162
 163   (frob %atan #'cl:atan)
 164
 165   (frob %asec #'(lambda (x)
 166                   (let ((y (ignore-errors (maxima-branch-acos (/ 1 x)))))
 167                     (if y y (domain-error x 'asec)))))
 168
 169   (frob %acsc #'(lambda (x)
 170                   (let ((y (ignore-errors (maxima-branch-asin (/ 1 x)))))
 171                     (if y y (domain-error x 'acsc)))))
 172
 173   (frob %acot #'(lambda (x)
 174                   (let ((y (ignore-errors (maxima-branch-acot x))))
 175                     (if y y (domain-error x 'acot)))))
 176
 177   (frob %cosh #'cl:cosh)
 178   (frob %sinh #'cl:sinh)
 179   (frob %tanh #'cl:tanh)
 180
 181   (frob %sech #'(lambda (x)
 182                   (flet ((sech (x)
 183                            ;; For large x > 0, cosh(x) ~= exp(x)/2.
 184                            ;; Hence, sech(x) ~= 2*exp(-x).  And since
 185                            ;; cosh(x) is even, we only need to deal
 186                            ;; with |x|.  Note also that if |x| >=
 187                            ;; sqrt(most-positive-double-float),
 188                            ;; exp(-x) is basically zero, so we can use
 189                            ;; a threshold of sqrt(most-positive).
 190                            ;;
 191                            ;; Several Lisp's can not compute acosh()
 192                            ;; for very large values, e.g.
 193                            ;; (acosh most-positive-double-float)
 194                            ;; Therefore use the numerical value
 195                            ;; 710.4758600739439d0 = (acosh most-positive-double-float)
 196                            ;; instead of computing the value.
 197                            ;; The most-positive-double-float is standardized (IEEE 754).
 198                            (if (and (floatp x)
 199                                     (>= (abs x) 710.4758600739439d0))
 200                                (* 2 (exp (- (abs x))))
 201                                (/ (cl:cosh x)))))
 202                   (let ((y (ignore-errors (sech x))))
 203                     (if y y (domain-error x 'sech))))))
 204
 205   (frob %csch #'(lambda (x)
 206                   (flet ((csch (x)
 207                            ;; For large x > 0, sinh(x) ~= exp(x)/2.
 208                            ;; Hence csch(x) = 2*exp(-x).  Since
 209                            ;; sinh(x) is odd, we also have csch(x) =
 210                            ;; -2*exp(x) when x < 0 and |x| is large.
 211                            ;;
 212                            ;; Several Lisp's can not compute asinh()
 213                            ;; for very large values, e.g.
 214                            ;; (asinh most-positive-double-float)
 215                            ;; Therefore use the numerical value
 216                            ;; 710.4758600739439d0 = (asinh most-positive-double-float)
 217                            ;; instead of computing the value.
 218                            ;; The most-positive-double-float is standardized (IEEE 754).
 219                            (if (and (floatp x)
 220                                     (>= (abs x) 710.4758600739439d0))
 221                                (float-sign x (* 2 (exp (- (abs x)))))
 222                                (/ (cl:sinh x)))))
 223                   (let ((y (ignore-errors (csch x))))
 224                     (if y y (domain-error x 'csch))))))
 225
 226   (frob %coth #'(lambda (x)
 227                   (let ((y (ignore-errors (/ 1 (cl:tanh x)))))
 228                     (if y y (domain-error x 'coth)))))
 229
 230   (frob %acosh #'cl:acosh)
 231   (frob %asinh #'cl:asinh)
 232
 233   (frob %atanh #'maxima-branch-atanh)
 234
 235   (frob %asech #'(lambda (x)
 236                    (let ((y (ignore-errors (cl:acosh (/ 1 x)))))
 237                      (if y y (domain-error x 'asech)))))
 238
 239   (frob %acsch
 240         #'(lambda (x)
 241             (flet ((acsch (x)
 242                      ;; logarc(acsch(x)) = log(1/x+sqrt(1/x^2+1)).
 243                      ;; Assume x > 0.  Then we can rewrite this as
 244                      ;; log((1+sqrt(1+x^2))/x).  If we choose x such
 245                      ;; that 1+x^2 = 1, then this simplifies to
 246                      ;; log(2/x).  However for very small x, 2/x can
 247                      ;; overflow, so use log(2)-log(x).
 248                      ;;
 249                      ;; 1+x^2 = 1 when x^2 = double-float-epsilon.  So
 250                      ;; x = sqrt(double-float-epsilon).  We'd really
 251                      ;; like to use
 252                      ;; least-positive-normalized-double-float, but
 253                      ;; some lisps like clisp don't have denormals.
 254                      ;; In that case, use sqrt(double-float-epsilon).
 255                      (let ((absx (abs x)))
 256                        (cond ((and (floatp x)
 257                                    (< absx
 258                                       #-clisp
 259                                       least-positive-normalized-double-float
 260                                       #+clisp
 261                                       (sqrt double-float-epsilon)))
 262                               (float-sign x (- (log 2d0) (log (abs x)))))
 263                              (t
 264                               (cl:asinh (/ x)))))))
 265               (let ((y (ignore-errors (acsch x))))
 266                 (if y y (domain-error x 'acsch))))))
 267
 268   (frob %acoth #'(lambda (x)
 269                    (let ((y (ignore-errors (maxima-branch-atanh (/ 1 x)))))
 270                      (if y y (domain-error x 'acoth)))))
 271
 272   (frob mabs #'cl:abs)
 273   (frob %exp #'cl:exp)
 274   (frob mexpt #'cl:expt)
 275   (frob %sqrt #'cl:sqrt)
 276   (frob %log #'(lambda (x)
 277                  (let ((y (ignore-errors (cl:log x))))
 278                    (if y y (domain-error x 'log)))))
 279
 280   (frob %plog #'(lambda (x)
 281                   (let ((y (ignore-errors (cl:log x))))
 282                     (if y y (domain-error x 'log)))))
 283
 284   (frob $conjugate #'cl:conjugate)
 285   (frob $floor #'cl:ffloor)
 286   (frob $ceiling #'cl:fceiling)
 287   (frob $realpart #'cl:realpart)
 288   (frob $imagpart #'cl:imagpart)
 289   (frob $max #'cl:max)
 290   (frob $min #'cl:min)
 291   (frob %signum #'cl:signum)
 292   (frob %atan2 #'cl:atan))
 293
 294 (macrolet ((frob (mfun dfun) `(setf (gethash ',mfun *big-float-op*) ,dfun)))
 295   ;; All big-float implementation functions MUST support a required x
 296   ;; arg and an optional y arg for the real and imaginary parts.  The
 297   ;; imaginary part does not have to be given.
 298   (frob %asin #'big-float-asin)
 299   (frob %sinh #'big-float-sinh)
 300   (frob %asinh #'big-float-asinh)
 301   (frob %tanh #'big-float-tanh)
 302   (frob %atanh #'big-float-atanh)
 303   (frob %acos 'big-float-acos)
 304   (frob %log 'big-float-log)
 305   (frob %sqrt 'big-float-sqrt))
 306
 307 ;; Here is a general scheme for defining and applying reflection rules. A
 308 ;; reflection rule is something like f(-x) --> f(x), or  f(-x) --> %pi - f(x).
 309
 310 ;; We define functions for the two most common reflection rules; these
 311 ;; are the odd function rule (f(-x) --> -f(x)) and the even function rule
 312 ;; (f(-x) --> f(x)). A reflection rule takes two arguments (the operator and
 313 ;; the operand).
 314
 315 (defun odd-function-reflect (op x)
 316   (neg (take (list op) (neg x))))
 317
 318 (defun even-function-reflect (op x)
 319   (take (list op) (neg x)))
 320
 321 ;; Put the reflection rule on the property list of the exponential-like
 322 ;; functions.
 323
 324 (setf (get '%cos 'reflection-rule) 'even-function-reflect)
 325 (setf (get '%sin 'reflection-rule) 'odd-function-reflect)
 326 (setf (get '%tan 'reflection-rule) 'odd-function-reflect)
 327 (setf (get '%sec 'reflection-rule) 'even-function-reflect)
 328 (setf (get '%csc 'reflection-rule) 'odd-function-reflect)
 329 (setf (get '%cot 'reflection-rule) 'odd-function-reflect)
 330
 331 ;; See A&S 4.4.14--4.4.19
 332
 333 (setf (get '%acos 'reflection-rule) #'(lambda (op x) (sub '$%pi (take (list op) (neg x)))))
 334 (setf (get '%asin 'reflection-rule) 'odd-function-reflect)
 335 (setf (get '%atan 'reflection-rule) 'odd-function-reflect)
 336 (setf (get '%asec 'reflection-rule) #'(lambda (op x) (sub '$%pi (take (list op) (neg x)))))
 337 (setf (get '%acsc 'reflection-rule) 'odd-function-reflect)
 338 (setf (get '%acot 'reflection-rule) 'odd-function-reflect)
 339
 340 (setf (get '%cosh 'reflection-rule) 'even-function-reflect)
 341 (setf (get '%sinh 'reflection-rule) 'odd-function-reflect)
 342 (setf (get '%tanh 'reflection-rule) 'odd-function-reflect)
 343 (setf (get '%sech 'reflection-rule) 'even-function-reflect)
 344 (setf (get '%csch 'reflection-rule) 'odd-function-reflect)
 345 (setf (get '%coth 'reflection-rule) 'odd-function-reflect)
 346
 347 (setf (get '%asinh 'reflection-rule) 'odd-function-reflect)
 348 (setf (get '%atanh 'reflection-rule) 'odd-function-reflect)
 349 (setf (get '%asech 'reflection-rule) 'even-function-reflect)
 350 (setf (get '%acsch 'reflection-rule) 'odd-function-reflect)
 351 (setf (get '%acoth 'reflection-rule) 'odd-function-reflect)
 352
 353 ;; When b is nil, do not apply the reflection rule. For trigonometric like
 354 ;; functions, b is $trigsign.  This function uses 'great' to decide when to
 355 ;; apply the rule.  Another possibility is to apply the rule when (mminusp* x)
 356 ;; evaluates to true. Maxima <= 5.9.3 uses this scheme; with this method, we have
 357 ;; assume(z < 0), cos(z) --> cos(-z). I (Barton Willis) think this goofy.
 358
 359 ;; The function 'great' is non-transitive. I don't think this bug will cause
 360 ;; trouble for this function. If there is an expression such that both
 361 ;; (great (neg x) x) and (great x (neg x)) evaluate to true, this function
 362 ;; could cause an infinite loop. I could protect against this possibility with
 363 ;; (and b f (great (neg x) x) (not (great x (neg x))).
 364
 365 (defun apply-reflection-simp (op x &optional (b t))
 366   (let ((f (get op 'reflection-rule)))
 367     (if (and b f (great (neg x) x)) (funcall f op x) nil)))
 368
 369 (defun taylorize (op x)
 370   (if ($taylorp x)
 371       (mfuncall '$apply '$taylor `((mlist) ((,op) ,($ratdisrep x)) ,@(cdr ($taylorinfo x)))) nil))
 372
 373 (defun float-or-rational-p (x)
 374   (or (floatp x) ($ratnump x)))
 375
 376 (defun bigfloat-or-number-p (x)
 377   (or ($bfloatp x) (numberp x) ($ratnump x)))
 378
 379 ;; When z is a Maxima complex float or when 'numer' is true and z is a
 380 ;; Maxima complex number, evaluate (op z) by applying the mapping from
 381 ;; the Maxima operator 'op' to the operator in the hash table
 382 ;; '*flonum-op*'. When z isn't a Maxima complex number, return
 383 ;; nil.
 384
 385 (defun flonum-eval (op z)
 386   (let ((op (gethash op *flonum-op*)))
 387     (when op
 388       (multiple-value-bind (bool R I)
 389         (complex-number-p z #'float-or-rational-p)
 390         (when (and bool (or $numer (floatp R) (floatp I)))
 391           (setq R ($float R))
 392           (setq I ($float I))
 393           (complexify (funcall op (if (zerop I) R (complex R I)))))))))
 394
 395 ;; For now, big float evaluation of trig-like functions for complex
 396 ;; big floats uses rectform.  I suspect that for some functions (not
 397 ;; all of them) rectform generates expressions that are poorly suited
 398 ;; for numerical evaluation. For better accuracy, these functions
 399 ;; (maybe acosh, for one) may need to be special cased.  If they are
 400 ;; special-cased, the *big-float-op* hash table contains the special
 401 ;; cases.
 402
 403 (defun big-float-eval (op z)
 404   (when (complex-number-p z 'bigfloat-or-number-p)
 405     (let ((x ($realpart z))
 406           (y ($imagpart z))
 407           (bop (gethash op *big-float-op*)))
 408       ;; If bop is non-NIL, we want to try that first.  If bop
 409       ;; declines (by returning NIL), we silently give up and use the
 410       ;; rectform version.
 411       (cond ((and ($bfloatp x) (like 0 y))
 412              (or (and bop (funcall bop x))
 413                  ($bfloat `((,op simp) ,x))))
 414             ((or ($bfloatp x) ($bfloatp y))
 415              (or (and bop (funcall bop ($bfloat x) ($bfloat y)))
 416                  (let ((z (add ($bfloat x) (mul '$%i ($bfloat y)))))
 417                    (setq z ($rectform `((,op simp) ,z)))
 418                    ($bfloat z))))))))
 419
 420 ;; For complex big float evaluation, it's important to check the
 421 ;; simp flag -- otherwise Maxima can get stuck in an infinite loop:
 422 ;; asin(1.23b0 + %i * 4.56b0) ---> (simp-%asin ((%asin) ...) -->
 423 ;; (big-float-eval ((%asin) ...) --> (risplit ((%asin simp) ...) -->
 424 ;; (simp-%asin ((%asin simp) ...). If the simp flag is ignored, we've
 425 ;; got trouble.
 426
 427 (def-simplifier sin (y)
 428   (let (z)
 429     (cond ((flonum-eval (mop form) y))
 430           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 431           ((taylorize (mop form) (second form)))
 432           ((and $%piargs (cond ((zerop1 y) 0)
 433                                ((has-const-or-int-term y '$%pi) (%piargs-sin/cos y)))))
 434           ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%sinh (coeff y '$%i 1))))
 435           ((and $triginverses (not (atom y))
 436                 (cond ((eq '%asin (setq z (caar y))) (cadr y))
 437                       ((eq '%acos z) (sqrt1-x^2 (cadr y)))
 438                       ((eq '%atan z) (div (cadr y) (sqrt1+x^2 (cadr y))))
 439                       ((eq '%acot z) (div 1 (mul (sqrt1+x^2 (div 1 (cadr y))) (cadr y))))
 440                       ((eq '%asec z) (sqrt1-x^2 (div 1 (cadr y))))
 441                       ((eq '%acsc z) (div 1 (cadr y)))
 442                       ((eq '%atan2 z) (div (cadr y) (sq-sumsq (cadr y) (caddr y)))))))
 443           ((and $trigexpand (trigexpand '%sin y)))
 444           ($exponentialize (exponentialize '%sin y))
 445           ((and $halfangles (halfangle '%sin y)))
 446           ((apply-reflection-simp (mop form) y $trigsign))
 447                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%sin (neg y))))
 448           (t (give-up)))))
 449
 450 (def-simplifier cos (y)
 451   (let (z)
 452     (cond ((flonum-eval (mop form) y))
 453           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 454           ((taylorize (mop form) (second form)))
 455           ((and $%piargs (cond ((zerop1 y) 1)
 456                                ((has-const-or-int-term y '$%pi)
 457                                 (%piargs-sin/cos (add %pi//2 y))))))
 458           ((and $%iargs (multiplep y '$%i)) (ftake* '%cosh (coeff y '$%i 1)))
 459           ((and $triginverses (not (atom y))
 460                 (cond ((eq '%acos (setq z (caar y))) (cadr y))
 461                       ((eq '%asin z) (sqrt1-x^2 (cadr y)))
 462                       ((eq '%atan z) (div 1 (sqrt1+x^2 (cadr y))))
 463                       ((eq '%acot z) (div 1 (sqrt1+x^2 (div 1 (cadr y)))))
 464                       ((eq '%asec z) (div 1 (cadr y)))
 465                       ((eq '%acsc z) (sqrt1-x^2 (div 1 (cadr y))))
 466                       ((eq '%atan2 z) (div (caddr y) (sq-sumsq (cadr y) (caddr y)))))))
 467           ((and $trigexpand (trigexpand '%cos y)))
 468           ($exponentialize (exponentialize '%cos y))
 469           ((and $halfangles (halfangle '%cos y)))
 470           ((apply-reflection-simp (mop form) y $trigsign))
 471                                         ;((and $trigsign (mminusp* y)) (ftake* '%cos (neg y)))
 472           (t (give-up)))))
 473
 474 (defun %piargs-sin/cos (x)
 475   (let ($float coeff ratcoeff zl-rem)
 476     (setq ratcoeff (get-const-or-int-terms x '$%pi)
 477           coeff (linearize ratcoeff)
 478           zl-rem (get-not-const-or-int-terms x '$%pi))
 479     (cond ((zerop1 zl-rem) (%piargs coeff ratcoeff))
 480           ((not (mevenp (car coeff))) nil)
 481           ((equal 0 (setq x (mmod (cdr coeff) 2))) (ftake* '%sin zl-rem))
 482           ((equal 1 x) (neg (ftake* '%sin zl-rem)))
 483           ((alike1 1//2 x) (ftake* '%cos zl-rem))
 484           ((alike1 '((rat) 3 2) x) (neg (ftake* '%cos zl-rem))))))
 485
 486
 487 (defun filter-sum (pred form simp-flag)
 488   "Takes form to be a sum and a sum of the summands for which pred is
 489    true. Passes simp-flag through to addn if there is more than one
 490    term in the sum."
 491   (if (mplusp form)
 492       (addn (mapcan
 493              #'(lambda (term)
 494                  (when (funcall pred term) (list term))) (cdr form))
 495             simp-flag)
 496     (if (funcall pred form) form 0)))
 497
 498 ;; collect terms of form A*var where A is a constant or integer.
 499 ;; returns sum of all such A.
 500 ;; does not expand form, so does not find constant term in (x+1)*var.
 501 ;; thus we cannot simplify sin(2*%pi*(1+x)) => sin(2*%pi*x) unless
 502 ;;  the user calls expand.  this could be extended to look a little
 503 ;;  more deeply into the expression, but we don't want to call expand
 504 ;;  in the core simplifier for reasons of speed and predictability.
 505 (defun get-const-or-int-terms (form var)
 506   (coeff
 507    (filter-sum (lambda (term)
 508                  (let ((coeff (coeff term var 1)))
 509                    (and (not (zerop1 coeff))
 510                         (or ($constantp coeff)
 511                             (maxima-integerp coeff)))))
 512                form
 513                0)
 514    var 1))
 515
 516 ;; collect terms skipped by get-const-or-int-terms
 517 (defun get-not-const-or-int-terms (form var)
 518   (filter-sum (lambda (term)
 519                 (let ((coeff (coeff term var 1)))
 520                   (not (and (not (zerop1 coeff))
 521                             (or ($constantp coeff)
 522                                 (maxima-integerp coeff))))))
 523               form
 524               0))
 525
 526 (defun has-const-or-int-term (form var)
 527   "Tests whether form has at least some term of the form a*var where a
 528   is constant or integer"
 529   (not (zerop1 (get-const-or-int-terms form var))))
 530
 531 (def-simplifier tan (y)
 532   (let (z)
 533     (cond ((flonum-eval (mop form) y))
 534           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 535           ((taylorize (mop form) (second form)))
 536           ((and $%piargs (cond ((zerop1 y) 0)
 537                                ((has-const-or-int-term y '$%pi) (%piargs-tan/cot y)))))
 538           ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%tanh (coeff y '$%i 1))))
 539           ((and $triginverses (not (atom y))
 540                 (cond ((eq '%atan (setq z (caar y))) (cadr y))
 541                       ((eq '%asin z) (div (cadr y) (sqrt1-x^2 (cadr y))))
 542                       ((eq '%acos z) (div (sqrt1-x^2 (cadr y)) (cadr y)))
 543                       ((eq '%acot z) (div 1 (cadr y)))
 544                       ((eq '%asec z) (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y)))
 545                       ((eq '%acsc z) (div 1 (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y))))
 546                       ((eq '%atan2 z) (div (cadr y) (caddr y))))))
 547           ((and $trigexpand (trigexpand '%tan y)))
 548           ($exponentialize (exponentialize '%tan y))
 549           ((and $halfangles (halfangle '%tan y)))
 550           ((apply-reflection-simp (mop form) y $trigsign))
 551                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%tan (neg y))))
 552           (t (give-up)))))
 553
 554 (def-simplifier cot (y)
 555   (let (z)
 556     (cond ((flonum-eval (mop form) y))
 557           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 558           ((taylorize (mop form) (second form)))
 559           ((and $%piargs (cond ((zerop1 y) (domain-error y 'cot))
 560                                ((and (has-const-or-int-term y '$%pi)
 561                                      (setq z (%piargs-tan/cot (add %pi//2 y))))
 562                                 (neg z)))))
 563           ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%coth (coeff y '$%i 1))))
 564           ((and $triginverses (not (atom y))
 565                 (cond ((eq '%acot (setq z (caar y))) (cadr y))
 566                       ((eq '%asin z) (div (sqrt1-x^2 (cadr y)) (cadr y)))
 567                       ((eq '%acos z) (div (cadr y) (sqrt1-x^2 (cadr y))))
 568                       ((eq '%atan z) (div 1 (cadr y)))
 569                       ((eq '%asec z) (div 1 (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y))))
 570                       ((eq '%acsc z) (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y)))
 571                       ((eq '%atan2 z) (div (caddr y) (cadr y))))))
 572           ((and $trigexpand (trigexpand '%cot y)))
 573           ($exponentialize (exponentialize '%cot y))
 574           ((and $halfangles (halfangle '%cot y)))
 575           ((apply-reflection-simp (mop form) y $trigsign))
 576                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%cot (neg y))))
 577           (t (give-up)))))
 578
 579 (defun %piargs-tan/cot (x)
 580   "If x is of the form tan(u) where u has a nonzero constant linear
 581    term in %pi, then %piargs-tan/cot returns a simplified version of x
 582    without this constant term."
 583   ;; Set coeff to be the coefficient of $%pi collecting terms with no
 584   ;; other atoms, so given %pi(x+1/2), coeff = 1/2. Let zl-rem be the
 585   ;; remainder (TODO: computing zl-rem could probably be prettier.)
 586   (let* ((nice-terms (get-const-or-int-terms x '$%pi))
 587          (coeff (linearize nice-terms))
 588          (zl-rem (get-not-const-or-int-terms x '$%pi))
 589          (sin-of-coeff-pi)
 590          (cos-of-coeff-pi))
 591     (cond
 592      ;; sin-of-coeff-pi and cos-of-coeff-pi are only non-nil if they
 593      ;; are constants that %piargs-offset could compute, and we just
 594      ;; checked that cos-of-coeff-pi was nonzero. Thus we can just
 595      ;; return their quotient.
 596      ((and (zerop1 zl-rem)
 597            (setq sin-of-coeff-pi
 598                  (%piargs coeff nil)))
 599       (setq cos-of-coeff-pi
 600             (%piargs (cons (car coeff)
 601                            (rplus 1//2 (cdr coeff))) nil))
 602       (cond ((zerop1 sin-of-coeff-pi)
 603              0)         ;; tan(integer*%pi)
 604             ((zerop1 cos-of-coeff-pi)
 605              (merror (intl:gettext "tan: ~M isn't in the domain of tan.") x))
 606             (cos-of-coeff-pi
 607              (div sin-of-coeff-pi cos-of-coeff-pi))))
 608
 609      ;; This expression sets x to the coeff of %pi (mod 1) as a side
 610      ;; effect and then, if this is zero, returns tan of the
 611      ;; rest, because tan has periodicity %pi.
 612      ((zerop1 (setq x (mmod (cdr coeff) 1)))
 613       (ftake* '%tan zl-rem))
 614
 615      ;; Similarly, if x = 1/2 then return -cot(x).
 616      ((alike1 1//2 x)
 617         (neg (ftake* '%cot zl-rem))))))
 618
 619 (def-simplifier csc (y)
 620   (let (z)
 621     (cond ((flonum-eval (mop form) y))
 622           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 623           ((taylorize (mop form) (second form)))
 624           ((and $%piargs (cond ((zerop1 y) (domain-error y 'csc))
 625                                ((has-const-or-int-term y '$%pi) (%piargs-csc/sec y)))))
 626           ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%csch (coeff y '$%i 1))))
 627           ((and $triginverses (not (atom y))
 628                 (cond ((eq '%acsc (setq z (caar y))) (cadr y))
 629                       ((eq '%asin z) (div 1 (cadr y)))
 630                       ((eq '%acos z) (div 1 (sqrt1-x^2 (cadr y))))
 631                       ((eq '%atan z) (div (sqrt1+x^2 (cadr y)) (cadr y)))
 632                       ((eq '%acot z) (mul (sqrt1+x^2 (div 1 (cadr y))) (cadr y)))
 633                       ((eq '%asec z) (div 1 (sqrt1-x^2 (div 1 (cadr y)))))
 634                       ((eq '%atan2 z) (div (sq-sumsq (cadr y) (caddr y)) (cadr y))))))
 635           ((and $trigexpand (trigexpand '%csc y)))
 636           ($exponentialize (exponentialize '%csc y))
 637           ((and $halfangles (halfangle '%csc y)))
 638           ((apply-reflection-simp (mop form) y $trigsign))
 639                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%csc (neg y))))
 640
 641           (t (give-up)))))
 642
 643 (def-simplifier sec (y)
 644   (let (z)
 645     (cond ((flonum-eval (mop form) y))
 646           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 647           ((taylorize (mop form) (second form)))
 648           ((and $%piargs (cond ((zerop1 y) 1)
 649                                ((has-const-or-int-term y '$%pi) (%piargs-csc/sec (add %pi//2 y))))))
 650           ((and $%iargs (multiplep y '$%i)) (ftake* '%sech (coeff y '$%i 1)))
 651           ((and $triginverses (not (atom y))
 652                 (cond ((eq '%asec (setq z (caar y))) (cadr y))
 653                       ((eq '%asin z) (div 1 (sqrt1-x^2 (cadr y))))
 654                       ((eq '%acos z) (div 1 (cadr y)))
 655                       ((eq '%atan z) (sqrt1+x^2 (cadr y)))
 656                       ((eq '%acot z) (sqrt1+x^2 (div 1 (cadr y))))
 657                       ((eq '%acsc z) (div 1 (sqrt1-x^2 (div 1 (cadr y)))))
 658                       ((eq '%atan2 z) (div (sq-sumsq (cadr y) (caddr y)) (caddr y))))))
 659           ((and $trigexpand (trigexpand '%sec y)))
 660           ($exponentialize (exponentialize '%sec y))
 661           ((and $halfangles (halfangle '%sec y)))
 662           ((apply-reflection-simp (mop form) y $trigsign))
 663                                         ;((and $trigsign (mminusp* y)) (ftake* '%sec (neg y)))
 664
 665           (t (give-up)))))
 666
 667 (defun %piargs-csc/sec (x)
 668   (prog ($float coeff ratcoeff zl-rem)
 669      (setq ratcoeff (get-const-or-int-terms x '$%pi)
 670            coeff (linearize ratcoeff)
 671            zl-rem (get-not-const-or-int-terms x '$%pi))
 672      (return (cond ((and (zerop1 zl-rem) (setq zl-rem (%piargs coeff nil))) (div 1 zl-rem))
 673                    ((not (mevenp (car coeff))) nil)
 674                    ((equal 0 (setq x (mmod (cdr coeff) 2))) (ftake* '%csc zl-rem))
 675                    ((equal 1 x) (neg (ftake* '%csc zl-rem)))
 676                    ((alike1 1//2 x) (ftake* '%sec zl-rem))
 677                    ((alike1 '((rat) 3 2) x) (neg (ftake* '%sec zl-rem)))))))
 678
 679 (def-simplifier atan (y)
 680   (cond ((flonum-eval (mop form) y))
 681         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 682         ((taylorize (mop form) (second form)))
 683         ;; Simplification for special values
 684         ((zerop1 y) y)
 685         ((or (eq y '$inf) (alike1 y '((mtimes) -1 $minf)))
 686          (div '$%pi 2))
 687         ((or (eq y '$minf) (alike1 y '((mtimes) -1 $inf)))
 688          (div '$%pi -2))
 689         ((and $%piargs
 690               ;; Recognize more special values
 691               (cond ((equal 1 y) (div '$%pi 4))
 692                     ((equal -1 y) (div '$%pi -4))
 693                     ;; sqrt(3)
 694                     ((alike1 y '((mexpt) 3 ((rat) 1 2)))
 695                      (div '$%pi 3))
 696                     ;; -sqrt(3)
 697                     ((alike1 y '((mtimes) -1 ((mexpt) 3 ((rat) 1 2))))
 698                      (div '$%pi -3))
 699                     ;; 1/sqrt(3)
 700                     ((alike1 y '((mexpt) 3 ((rat) -1 2)))
 701                      (div '$%pi 6))
 702                     ;; -1/sqrt(3)
 703                     ((alike1 y '((mtimes) -1 ((mexpt) 3 ((rat) -1 2))))
 704                      (div '$%pi -6))
 705                     ((alike1 y '((mplus) -1 ((mexpt) 2 ((rat) 1 2))))
 706                      (div '$%pi 8))
 707                     ((alike1 y '((mplus) 1 ((mexpt) 2 ((rat) 1 2))))
 708                      (mul 3 (div '$%pi 8))))))
 709         ((and $%iargs (multiplep y '$%i))
 710          ;; atan(%i*y) -> %i*atanh(y)
 711          (mul '$%i (take '(%atanh) (coeff y '$%i 1))))
 712         ((and (not (atom y)) (member (caar y) '(%cot %tan))
 713               (if ($constantp (cadr y))
 714                   (let ((y-val (mfuncall '$mod
 715                                          (if (eq (caar y) '%tan)
 716                                              (cadr y)
 717                                              (sub %pi//2 (cadr y)))
 718                                          '$%pi)))
 719                     (cond ((eq (mlsp y-val %pi//2) t) y-val)
 720                           ((eq (mlsp y-val '$%pi) t) (sub y-val '$%pi)))))))
 721         ((and (eq $triginverses '$all) (not (atom y))
 722               (if (eq (caar y) '%tan) (cadr y))))
 723         ((and (eq $triginverses t) (not (atom y)) (eq (caar y) '%tan)
 724               ;; Check if y in [-%pi/2, %pi/2]
 725               (if (and (member (csign (sub (cadr y) %pi//2)) '($nz $neg) :test #'eq)
 726                        (member (csign (add (cadr y) %pi//2)) '($pz $pos) :test #'eq))
 727                   (cadr y))))
 728         ($logarc (logarc '%atan y))
 729         ((apply-reflection-simp (mop form) y $trigsign))
 730         (t (give-up))))
 731
 732 (defun %piargs (x ratcoeff)
 733   (let (offset-result)
 734     (cond ((and (integerp (car x)) (integerp (cdr x))) 0)
 735           ((not (mevenp (car x)))
 736            (cond ((null ratcoeff) nil)
 737                  ((and (integerp (car x))
 738                        (setq offset-result (%piargs-offset (cdr x))))
 739                   (mul (power -1 (sub ratcoeff (cdr x)))
 740                        offset-result))))
 741           ((%piargs-offset (mmod (cdr x) 2))))))
 742
 743 ; simplifies sin(%pi * x) where x is between 0 and 1
 744 ; returns nil if can't simplify
 745 (defun %piargs-offset (x)
 746   (cond ((or (alike1 '((rat) 1 6) x) (alike1 '((rat) 5 6) x)) 1//2)
 747         ((or (alike1 '((rat) 1 4) x) (alike1 '((rat) 3 4) x)) (div (power 2 1//2) 2))
 748         ((or (alike1 '((rat) 1 3) x) (alike1 '((rat) 2 3) x)) (div (power 3 1//2) 2))
 749         ((alike1 1//2 x) 1)
 750         ((or (alike1 '((rat) 7 6) x) (alike1 '((rat) 11 6) x)) -1//2)
 751         ((or (alike1 '((rat) 4 3) x) (alike1 '((rat) 5 3) x)) (div (power 3 1//2) -2))
 752         ((or (alike1 '((rat) 5 4) x) (alike1 '((rat) 7 4) x)) (mul -1//2 (power 2 1//2)))
 753         ((alike1 '((rat) 3 2) x) -1)))
 754
 755 ;; identifies integer part of form
 756 ;; returns (X . Y) if form can be written as X*some_integer + Y
 757 ;; returns nil otherwise
 758 (defun linearize (form)
 759   (cond ((integerp form) (cons 0 form))
 760         ((numberp form) nil)
 761         ((atom form)
 762          (let (dum)
 763            (cond ((setq dum (evod form))
 764                   (if (eq '$even dum) '(2 . 0) '(2 . 1)))
 765                  ((maxima-integerp form) '(1 . 0)))))
 766         ((eq 'rat (caar form)) (cons 0 form))
 767         ((eq 'mplus (caar form)) (lin-mplus form))
 768         ((eq 'mtimes (caar form)) (lin-mtimes form))
 769         ((eq 'mexpt (caar form)) (lin-mexpt form))))
 770
 771 (defun lin-mplus (form)
 772   (do ((tl (cdr form) (cdr tl)) (dummy) (coeff 0) (zl-rem 0))
 773       ((null tl) (cons coeff (mmod zl-rem coeff)))
 774     (setq dummy (linearize (car tl)))
 775     (if (null dummy) (return nil)
 776         (setq coeff (rgcd (car dummy) coeff) zl-rem (rplus (cdr dummy) zl-rem)))))
 777
 778 (defun lin-mtimes (form)
 779   (do ((fl (cdr form) (cdr fl)) (dummy) (coeff 0) (zl-rem 1))
 780       ((null fl) (cons coeff (mmod zl-rem coeff)))
 781     (setq dummy (linearize (car fl)))
 782     (cond ((null dummy) (return nil))
 783           (t (setq coeff (rgcd (rtimes coeff (car dummy))
 784                                (rgcd (rtimes coeff (cdr dummy)) (rtimes zl-rem (car dummy))))
 785                    zl-rem (rtimes (cdr dummy) zl-rem))))))
 786
 787 (defun lin-mexpt (form)
 788   (prog (dummy)
 789      (cond ((and (integerp (caddr form)) (not (minusp (caddr form)))
 790                  (not (null (setq dummy (linearize (cadr form))))))
 791             (return (cons (car dummy) (mmod (cdr dummy) (caddr form))))))))
 792
 793 (defun rgcd (x y)
 794   (cond ((integerp x)
 795          (cond ((integerp y) (gcd x y))
 796                (t (list '(rat) (gcd x (cadr y)) (caddr y)))))
 797         ((integerp y) (list '(rat) (gcd (cadr x) y) (caddr x)))
 798         (t (list '(rat) (gcd (cadr x) (cadr y)) (lcm (caddr x) (caddr y))))))
 799
 800 (defun maxima-reduce (x y)
 801   (prog (gcd)
 802      (setq gcd (gcd x y) x (truncate x gcd) y (truncate y gcd))
 803      (if (minusp y) (setq x (- x) y (- y)))
 804      (return (if (eql y 1) x (list '(rat simp) x y)))))
 805
 806 ;; The following four functions are generated in code by TRANSL. - JPG 2/1/81
 807
 808 (defun rplus (x y) (addk x y))
 809
 810 (defun rdifference (x y) (addk x (timesk -1 y)))
 811
 812 (defun rtimes (x y) (timesk x y))
 813
 814 (defun rremainder (x y)
 815   (cond ((equal 0 y) (dbz-err))
 816         ((integerp x)
 817          (cond ((integerp y) (maxima-reduce x y))
 818                (t (maxima-reduce (* x (caddr y)) (cadr y)))))
 819         ((integerp y) (maxima-reduce (cadr x) (* (caddr x) y)))
 820         (t (maxima-reduce (* (cadr x) (caddr y)) (* (caddr x) (cadr y))))))
 821
 822 (defmfun $exponentialize (exp)
 823   (let ($demoivre)
 824     (cond ((atom exp) exp)
 825           ((trigp (caar exp))
 826            (exponentialize (caar exp) ($exponentialize (cadr exp))))
 827           (t (recur-apply #'$exponentialize exp)))))
 828
 829 (defun exponentialize (op arg)
 830   (cond ((eq '%sin op)
 831          (div (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))
 832               (mul 2 '$%i)))
 833         ((eq '%cos op)
 834          (div (add (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg))) 2))
 835         ((eq '%tan op)
 836          (div (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))
 837               (mul '$%i (add (power '$%e (mul '$%i arg))
 838                              (power '$%e (mul -1 '$%i arg))))))
 839         ((eq '%cot op)
 840          (div (mul '$%i (add (power '$%e (mul '$%i arg))
 841                              (power '$%e (mul -1 '$%i arg))))
 842               (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))))
 843         ((eq '%csc op)
 844          (div (mul 2 '$%i)
 845               (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))))
 846         ((eq '%sec op)
 847          (div 2 (add (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))))
 848         ((eq '%sinh op)
 849          (div (sub (power '$%e arg) (power '$%e (neg arg))) 2))
 850         ((eq '%cosh op)
 851          (div (add (power '$%e arg) (power '$%e (mul -1 arg))) 2))
 852         ((eq '%tanh op)
 853          (div (sub (power '$%e arg) (power '$%e (neg arg)))
 854               (add (power '$%e arg) (power '$%e (mul -1 arg)))))
 855         ((eq '%coth op)
 856          (div (add (power '$%e arg) (power '$%e (mul -1 arg)))
 857               (sub (power '$%e arg) (power '$%e (neg arg)))))
 858         ((eq '%csch op)
 859          (div 2 (sub (power '$%e arg) (power '$%e (neg arg)))))
 860         ((eq '%sech op)
 861          (div 2 (add (power '$%e arg) (power '$%e (mul -1 arg)))))))
 862
 863 (defun coefficient (exp var pow)
 864   (coeff exp var pow))
 865
 866 (defun mmod (x mod)
 867   (cond ((and (integerp x) (integerp mod))
 868          (if (minusp (if (zerop mod) x (setq x (- x (* mod (truncate x mod))))))
 869              (+ x mod)
 870              x))
 871         ((and ($ratnump x) ($ratnump mod))
 872          (let
 873              ((d (lcm ($denom x) ($denom mod))))
 874            (setq x (mul* d x))
 875            (setq mod (mul* d mod))
 876            (div (mod x mod) d)))
 877         (t nil)))
 878
 879 (defun multiplep (exp var)
 880   (and (not (zerop1 exp)) (zerop1 (sub exp (mul var (coeff exp var 1))))))
 881
 882 (defun linearp (exp var)
 883   (and (setq exp (islinear exp var)) (not (equal (car exp) 0))))
 884
 885 (defun mminusp (x)
 886   (= -1 (signum1 x)))
 887
 888 (defun mminusp* (x)
 889   (let (sign)
 890     (setq sign (csign x))
 891     (or (member sign '($neg $nz) :test #'eq)
 892         (and (mminusp x) (not (member sign '($pos $pz) :test #'eq))))))
 893
 894 ;; This should give more information somehow.
 895
 896 (defun dbz-err ()
 897   (cond ((not errorsw) (merror (intl:gettext "Division by zero attempted.")))
 898         (t (throw 'errorsw t))))
 899
 900 (defun dbz-err1 (func)
 901   (cond ((not errorsw) (merror (intl:gettext "~A: division by zero attempted.") func))
 902         (t (throw 'errorsw t))))