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   (frob %atan 'big-float-atan))
 307
 308 ;; Here is a general scheme for defining and applying reflection rules. A
 309 ;; reflection rule is something like f(-x) --> f(x), or  f(-x) --> %pi - f(x).
 310
 311 ;; We define functions for the two most common reflection rules; these
 312 ;; are the odd function rule (f(-x) --> -f(x)) and the even function rule
 313 ;; (f(-x) --> f(x)). A reflection rule takes two arguments (the operator and
 314 ;; the operand).
 315
 316 (defun odd-function-reflect (op x)
 317   (neg (take (list op) (neg x))))
 318
 319 (defun even-function-reflect (op x)
 320   (take (list op) (neg x)))
 321
 322 ;; Put the reflection rule on the property list of the exponential-like
 323 ;; functions.
 324
 325 (setf (get '%cos 'reflection-rule) 'even-function-reflect)
 326 (setf (get '%sin 'reflection-rule) 'odd-function-reflect)
 327 (setf (get '%tan 'reflection-rule) 'odd-function-reflect)
 328 (setf (get '%sec 'reflection-rule) 'even-function-reflect)
 329 (setf (get '%csc 'reflection-rule) 'odd-function-reflect)
 330 (setf (get '%cot 'reflection-rule) 'odd-function-reflect)
 331
 332 ;; See A&S 4.4.14--4.4.19
 333
 334 (setf (get '%acos 'reflection-rule) #'(lambda (op x) (sub '$%pi (take (list op) (neg x)))))
 335 (setf (get '%asin 'reflection-rule) 'odd-function-reflect)
 336 (setf (get '%atan 'reflection-rule) 'odd-function-reflect)
 337 (setf (get '%asec 'reflection-rule) #'(lambda (op x) (sub '$%pi (take (list op) (neg x)))))
 338 (setf (get '%acsc 'reflection-rule) 'odd-function-reflect)
 339 (setf (get '%acot 'reflection-rule) 'odd-function-reflect)
 340
 341 (setf (get '%cosh 'reflection-rule) 'even-function-reflect)
 342 (setf (get '%sinh 'reflection-rule) 'odd-function-reflect)
 343 (setf (get '%tanh 'reflection-rule) 'odd-function-reflect)
 344 (setf (get '%sech 'reflection-rule) 'even-function-reflect)
 345 (setf (get '%csch 'reflection-rule) 'odd-function-reflect)
 346 (setf (get '%coth 'reflection-rule) 'odd-function-reflect)
 347
 348 (setf (get '%asinh 'reflection-rule) 'odd-function-reflect)
 349 (setf (get '%atanh 'reflection-rule) 'odd-function-reflect)
 350 (setf (get '%asech 'reflection-rule) 'even-function-reflect)
 351 (setf (get '%acsch 'reflection-rule) 'odd-function-reflect)
 352 (setf (get '%acoth 'reflection-rule) 'odd-function-reflect)
 353
 354 ;; When b is nil, do not apply the reflection rule. For trigonometric like
 355 ;; functions, b is $trigsign.  This function uses 'great' to decide when to
 356 ;; apply the rule.  Another possibility is to apply the rule when (mminusp* x)
 357 ;; evaluates to true. Maxima <= 5.9.3 uses this scheme; with this method, we have
 358 ;; assume(z < 0), cos(z) --> cos(-z). I (Barton Willis) think this goofy.
 359
 360 ;; The function 'great' is non-transitive. I don't think this bug will cause
 361 ;; trouble for this function. If there is an expression such that both
 362 ;; (great (neg x) x) and (great x (neg x)) evaluate to true, this function
 363 ;; could cause an infinite loop. I could protect against this possibility with
 364 ;; (and b f (great (neg x) x) (not (great x (neg x))).
 365
 366 (defun apply-reflection-simp (op x &optional (b t))
 367   (let ((f (get op 'reflection-rule)))
 368     (if (and b f (great (neg x) x)) (funcall f op x) nil)))
 369
 370 (defun taylorize (op x)
 371   (if ($taylorp x)
 372       (mfuncall '$apply '$taylor `((mlist) ((,op) ,($ratdisrep x)) ,@(cdr ($taylorinfo x)))) nil))
 373
 374 (defun float-or-rational-p (x)
 375   (or (floatp x) ($ratnump x)))
 376
 377 (defun bigfloat-or-number-p (x)
 378   (or ($bfloatp x) (numberp x) ($ratnump x)))
 379
 380 ;; When z is a Maxima complex float or when 'numer' is true and z is a
 381 ;; Maxima complex number, evaluate (op z) by applying the mapping from
 382 ;; the Maxima operator 'op' to the operator in the hash table
 383 ;; '*flonum-op*'. When z isn't a Maxima complex number, return
 384 ;; nil.
 385
 386 (defun flonum-eval (op z)
 387   (let ((op (gethash op *flonum-op*)))
 388     (when op
 389       (multiple-value-bind (bool R I)
 390         (complex-number-p z #'float-or-rational-p)
 391         (when (and bool (or $numer (floatp R) (floatp I)))
 392           (setq R ($float R))
 393           (setq I ($float I))
 394           (complexify (funcall op (if (zerop I) R (complex R I)))))))))
 395
 396 ;; For now, big float evaluation of trig-like functions for complex
 397 ;; big floats uses rectform.  I suspect that for some functions (not
 398 ;; all of them) rectform generates expressions that are poorly suited
 399 ;; for numerical evaluation. For better accuracy, these functions
 400 ;; (maybe acosh, for one) may need to be special cased.  If they are
 401 ;; special-cased, the *big-float-op* hash table contains the special
 402 ;; cases.
 403
 404 (defun big-float-eval (op z)
 405   (when (complex-number-p z 'bigfloat-or-number-p)
 406     (let ((x ($realpart z))
 407           (y ($imagpart z))
 408           (bop (gethash op *big-float-op*)))
 409       ;; If bop is non-NIL, we want to try that first.  If bop
 410       ;; declines (by returning NIL), we silently give up and use the
 411       ;; rectform version.
 412       (cond ((and ($bfloatp x) (like 0 y))
 413              (or (and bop (funcall bop x))
 414                  ($bfloat `((,op simp) ,x))))
 415             ((or ($bfloatp x) ($bfloatp y))
 416              (or (and bop (funcall bop ($bfloat x) ($bfloat y)))
 417                  (let ((z (add ($bfloat x) (mul '$%i ($bfloat y)))))
 418                    (setq z ($rectform `((,op simp) ,z)))
 419                    ($bfloat z))))))))
 420
 421 ;; For complex big float evaluation, it's important to check the
 422 ;; simp flag -- otherwise Maxima can get stuck in an infinite loop:
 423 ;; asin(1.23b0 + %i * 4.56b0) ---> (simp-%asin ((%asin) ...) -->
 424 ;; (big-float-eval ((%asin) ...) --> (risplit ((%asin simp) ...) -->
 425 ;; (simp-%asin ((%asin simp) ...). If the simp flag is ignored, we've
 426 ;; got trouble.
 427
 428 (def-simplifier sin (y)
 429   (let (z)
 430     (cond ((flonum-eval (mop form) y))
 431           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 432           ((taylorize (mop form) (second form)))
 433           ((and $%piargs (cond ((zerop1 y) 0)
 434                                ((has-const-or-int-term y '$%pi) (%piargs-sin/cos y)))))
 435           ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%sinh (coeff y '$%i 1))))
 436           ((and $triginverses (not (atom y))
 437                 (cond ((eq '%asin (setq z (caar y))) (cadr y))
 438                       ((eq '%acos z) (sqrt1-x^2 (cadr y)))
 439                       ((eq '%atan z) (div (cadr y) (sqrt1+x^2 (cadr y))))
 440                       ((eq '%acot z) (div 1 (mul (sqrt1+x^2 (div 1 (cadr y))) (cadr y))))
 441                       ((eq '%asec z) (sqrt1-x^2 (div 1 (cadr y))))
 442                       ((eq '%acsc z) (div 1 (cadr y)))
 443                       ((eq '%atan2 z) (div (cadr y) (sq-sumsq (cadr y) (caddr y)))))))
 444           ((and $trigexpand (trigexpand '%sin y)))
 445           ($exponentialize (exponentialize '%sin y))
 446           ((and $halfangles (halfangle '%sin y)))
 447           ((apply-reflection-simp (mop form) y $trigsign))
 448                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%sin (neg y))))
 449           (t (give-up)))))
 450
 451 (def-simplifier cos (y)
 452   (let (z)
 453     (cond ((flonum-eval (mop form) y))
 454           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 455           ((taylorize (mop form) (second form)))
 456           ((and $%piargs (cond ((zerop1 y) 1)
 457                                ((has-const-or-int-term y '$%pi)
 458                                 (%piargs-sin/cos (add %pi//2 y))))))
 459           ((and $%iargs (multiplep y '$%i)) (ftake* '%cosh (coeff y '$%i 1)))
 460           ((and $triginverses (not (atom y))
 461                 (cond ((eq '%acos (setq z (caar y))) (cadr y))
 462                       ((eq '%asin z) (sqrt1-x^2 (cadr y)))
 463                       ((eq '%atan z) (div 1 (sqrt1+x^2 (cadr y))))
 464                       ((eq '%acot z) (div 1 (sqrt1+x^2 (div 1 (cadr y)))))
 465                       ((eq '%asec z) (div 1 (cadr y)))
 466                       ((eq '%acsc z) (sqrt1-x^2 (div 1 (cadr y))))
 467                       ((eq '%atan2 z) (div (caddr y) (sq-sumsq (cadr y) (caddr y)))))))
 468           ((and $trigexpand (trigexpand '%cos y)))
 469           ($exponentialize (exponentialize '%cos y))
 470           ((and $halfangles (halfangle '%cos y)))
 471           ((apply-reflection-simp (mop form) y $trigsign))
 472                                         ;((and $trigsign (mminusp* y)) (ftake* '%cos (neg y)))
 473           (t (give-up)))))
 474
 475 (defun %piargs-sin/cos (x)
 476   (let ($float coeff ratcoeff zl-rem)
 477     (setq ratcoeff (get-const-or-int-terms x '$%pi)
 478           coeff (linearize ratcoeff)
 479           zl-rem (get-not-const-or-int-terms x '$%pi))
 480     (cond ((zerop1 zl-rem) (%piargs coeff ratcoeff))
 481           ((not (mevenp (car coeff))) nil)
 482           ((equal 0 (setq x (mmod (cdr coeff) 2))) (ftake* '%sin zl-rem))
 483           ((equal 1 x) (neg (ftake* '%sin zl-rem)))
 484           ((alike1 1//2 x) (ftake* '%cos zl-rem))
 485           ((alike1 '((rat) 3 2) x) (neg (ftake* '%cos zl-rem))))))
 486
 487
 488 (defun filter-sum (pred form simp-flag)
 489   "Takes form to be a sum and a sum of the summands for which pred is
 490    true. Passes simp-flag through to addn if there is more than one
 491    term in the sum."
 492   (if (mplusp form)
 493       (addn (mapcan
 494              #'(lambda (term)
 495                  (when (funcall pred term) (list term))) (cdr form))
 496             simp-flag)
 497     (if (funcall pred form) form 0)))
 498
 499 ;; collect terms of form A*var where A is a constant or integer.
 500 ;; returns sum of all such A.
 501 ;; does not expand form, so does not find constant term in (x+1)*var.
 502 ;; thus we cannot simplify sin(2*%pi*(1+x)) => sin(2*%pi*x) unless
 503 ;;  the user calls expand.  this could be extended to look a little
 504 ;;  more deeply into the expression, but we don't want to call expand
 505 ;;  in the core simplifier for reasons of speed and predictability.
 506 (defun get-const-or-int-terms (form var)
 507   (coeff
 508    (filter-sum (lambda (term)
 509                  (let ((coeff (coeff term var 1)))
 510                    (and (not (zerop1 coeff))
 511                         (or ($constantp coeff)
 512                             (maxima-integerp coeff)))))
 513                form
 514                0)
 515    var 1))
 516
 517 ;; collect terms skipped by get-const-or-int-terms
 518 (defun get-not-const-or-int-terms (form var)
 519   (filter-sum (lambda (term)
 520                 (let ((coeff (coeff term var 1)))
 521                   (not (and (not (zerop1 coeff))
 522                             (or ($constantp coeff)
 523                                 (maxima-integerp coeff))))))
 524               form
 525               0))
 526
 527 (defun has-const-or-int-term (form var)
 528   "Tests whether form has at least some term of the form a*var where a
 529   is constant or integer"
 530   (not (zerop1 (get-const-or-int-terms form var))))
 531
 532 (def-simplifier tan (y)
 533   (let (z)
 534     (cond ((flonum-eval (mop form) y))
 535           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 536           ((taylorize (mop form) (second form)))
 537           ((and $%piargs (cond ((zerop1 y) 0)
 538                                ((has-const-or-int-term y '$%pi) (%piargs-tan/cot y)))))
 539           ((and $%iargs (multiplep y '$%i)) (mul '$%i (ftake* '%tanh (coeff y '$%i 1))))
 540           ((and $triginverses (not (atom y))
 541                 (cond ((eq '%atan (setq z (caar y))) (cadr y))
 542                       ((eq '%asin z) (div (cadr y) (sqrt1-x^2 (cadr y))))
 543                       ((eq '%acos z) (div (sqrt1-x^2 (cadr y)) (cadr y)))
 544                       ((eq '%acot z) (div 1 (cadr y)))
 545                       ((eq '%asec z) (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y)))
 546                       ((eq '%acsc z) (div 1 (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y))))
 547                       ((eq '%atan2 z) (div (cadr y) (caddr y))))))
 548           ((and $trigexpand (trigexpand '%tan y)))
 549           ($exponentialize (exponentialize '%tan y))
 550           ((and $halfangles (halfangle '%tan y)))
 551           ((apply-reflection-simp (mop form) y $trigsign))
 552                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%tan (neg y))))
 553           (t (give-up)))))
 554
 555 (def-simplifier cot (y)
 556   (let (z)
 557     (cond ((flonum-eval (mop form) y))
 558           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 559           ((taylorize (mop form) (second form)))
 560           ((and $%piargs (cond ((zerop1 y) (domain-error y 'cot))
 561                                ((and (has-const-or-int-term y '$%pi)
 562                                      (setq z (%piargs-tan/cot (add %pi//2 y))))
 563                                 (neg z)))))
 564           ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%coth (coeff y '$%i 1))))
 565           ((and $triginverses (not (atom y))
 566                 (cond ((eq '%acot (setq z (caar y))) (cadr y))
 567                       ((eq '%asin z) (div (sqrt1-x^2 (cadr y)) (cadr y)))
 568                       ((eq '%acos z) (div (cadr y) (sqrt1-x^2 (cadr y))))
 569                       ((eq '%atan z) (div 1 (cadr y)))
 570                       ((eq '%asec z) (div 1 (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y))))
 571                       ((eq '%acsc z) (mul (sqrt1-x^2 (div 1 (cadr y))) (cadr y)))
 572                       ((eq '%atan2 z) (div (caddr y) (cadr y))))))
 573           ((and $trigexpand (trigexpand '%cot y)))
 574           ($exponentialize (exponentialize '%cot y))
 575           ((and $halfangles (halfangle '%cot y)))
 576           ((apply-reflection-simp (mop form) y $trigsign))
 577                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%cot (neg y))))
 578           (t (give-up)))))
 579
 580 (defun %piargs-tan/cot (x)
 581   "If x is of the form tan(u) where u has a nonzero constant linear
 582    term in %pi, then %piargs-tan/cot returns a simplified version of x
 583    without this constant term."
 584   ;; Set coeff to be the coefficient of $%pi collecting terms with no
 585   ;; other atoms, so given %pi(x+1/2), coeff = 1/2. Let zl-rem be the
 586   ;; remainder (TODO: computing zl-rem could probably be prettier.)
 587   (let* ((nice-terms (get-const-or-int-terms x '$%pi))
 588          (coeff (linearize nice-terms))
 589          (zl-rem (get-not-const-or-int-terms x '$%pi))
 590          (sin-of-coeff-pi)
 591          (cos-of-coeff-pi))
 592     (cond
 593      ;; sin-of-coeff-pi and cos-of-coeff-pi are only non-nil if they
 594      ;; are constants that %piargs-offset could compute, and we just
 595      ;; checked that cos-of-coeff-pi was nonzero. Thus we can just
 596      ;; return their quotient.
 597      ((and (zerop1 zl-rem)
 598            (setq sin-of-coeff-pi
 599                  (%piargs coeff nil)))
 600       (setq cos-of-coeff-pi
 601             (%piargs (cons (car coeff)
 602                            (rplus 1//2 (cdr coeff))) nil))
 603       (cond ((zerop1 sin-of-coeff-pi)
 604              0)         ;; tan(integer*%pi)
 605             ((zerop1 cos-of-coeff-pi)
 606              (merror (intl:gettext "tan: ~M isn't in the domain of tan.") x))
 607             (cos-of-coeff-pi
 608              (div sin-of-coeff-pi cos-of-coeff-pi))))
 609
 610      ;; This expression sets x to the coeff of %pi (mod 1) as a side
 611      ;; effect and then, if this is zero, returns tan of the
 612      ;; rest, because tan has periodicity %pi.
 613      ((zerop1 (setq x (mmod (cdr coeff) 1)))
 614       (ftake* '%tan zl-rem))
 615
 616      ;; Similarly, if x = 1/2 then return -cot(x).
 617      ((alike1 1//2 x)
 618         (neg (ftake* '%cot zl-rem))))))
 619
 620 (def-simplifier csc (y)
 621   (let (z)
 622     (cond ((flonum-eval (mop form) y))
 623           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 624           ((taylorize (mop form) (second form)))
 625           ((and $%piargs (cond ((zerop1 y) (domain-error y 'csc))
 626                                ((has-const-or-int-term y '$%pi) (%piargs-csc/sec y)))))
 627           ((and $%iargs (multiplep y '$%i)) (mul -1 '$%i (ftake* '%csch (coeff y '$%i 1))))
 628           ((and $triginverses (not (atom y))
 629                 (cond ((eq '%acsc (setq z (caar y))) (cadr y))
 630                       ((eq '%asin z) (div 1 (cadr y)))
 631                       ((eq '%acos z) (div 1 (sqrt1-x^2 (cadr y))))
 632                       ((eq '%atan z) (div (sqrt1+x^2 (cadr y)) (cadr y)))
 633                       ((eq '%acot z) (mul (sqrt1+x^2 (div 1 (cadr y))) (cadr y)))
 634                       ((eq '%asec z) (div 1 (sqrt1-x^2 (div 1 (cadr y)))))
 635                       ((eq '%atan2 z) (div (sq-sumsq (cadr y) (caddr y)) (cadr y))))))
 636           ((and $trigexpand (trigexpand '%csc y)))
 637           ($exponentialize (exponentialize '%csc y))
 638           ((and $halfangles (halfangle '%csc y)))
 639           ((apply-reflection-simp (mop form) y $trigsign))
 640                                         ;((and $trigsign (mminusp* y)) (neg (ftake* '%csc (neg y))))
 641
 642           (t (give-up)))))
 643
 644 (def-simplifier sec (y)
 645   (let (z)
 646     (cond ((flonum-eval (mop form) y))
 647           ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 648           ((taylorize (mop form) (second form)))
 649           ((and $%piargs (cond ((zerop1 y) 1)
 650                                ((has-const-or-int-term y '$%pi) (%piargs-csc/sec (add %pi//2 y))))))
 651           ((and $%iargs (multiplep y '$%i)) (ftake* '%sech (coeff y '$%i 1)))
 652           ((and $triginverses (not (atom y))
 653                 (cond ((eq '%asec (setq z (caar y))) (cadr y))
 654                       ((eq '%asin z) (div 1 (sqrt1-x^2 (cadr y))))
 655                       ((eq '%acos z) (div 1 (cadr y)))
 656                       ((eq '%atan z) (sqrt1+x^2 (cadr y)))
 657                       ((eq '%acot z) (sqrt1+x^2 (div 1 (cadr y))))
 658                       ((eq '%acsc z) (div 1 (sqrt1-x^2 (div 1 (cadr y)))))
 659                       ((eq '%atan2 z) (div (sq-sumsq (cadr y) (caddr y)) (caddr y))))))
 660           ((and $trigexpand (trigexpand '%sec y)))
 661           ($exponentialize (exponentialize '%sec y))
 662           ((and $halfangles (halfangle '%sec y)))
 663           ((apply-reflection-simp (mop form) y $trigsign))
 664                                         ;((and $trigsign (mminusp* y)) (ftake* '%sec (neg y)))
 665
 666           (t (give-up)))))
 667
 668 (defun %piargs-csc/sec (x)
 669   (prog ($float coeff ratcoeff zl-rem)
 670      (setq ratcoeff (get-const-or-int-terms x '$%pi)
 671            coeff (linearize ratcoeff)
 672            zl-rem (get-not-const-or-int-terms x '$%pi))
 673      (return (cond ((and (zerop1 zl-rem) (setq zl-rem (%piargs coeff nil))) (div 1 zl-rem))
 674                    ((not (mevenp (car coeff))) nil)
 675                    ((equal 0 (setq x (mmod (cdr coeff) 2))) (ftake* '%csc zl-rem))
 676                    ((equal 1 x) (neg (ftake* '%csc zl-rem)))
 677                    ((alike1 1//2 x) (ftake* '%sec zl-rem))
 678                    ((alike1 '((rat) 3 2) x) (neg (ftake* '%sec zl-rem)))))))
 679
 680 (def-simplifier atan (y)
 681   (cond ((flonum-eval (mop form) y))
 682         ((and (not (member 'simp (car form))) (big-float-eval (mop form) y)))
 683         ((taylorize (mop form) (second form)))
 684         ;; Simplification for special values
 685         ((zerop1 y) y)
 686         ((or (eq y '$inf) (alike1 y '((mtimes) -1 $minf)))
 687          (div '$%pi 2))
 688         ((or (eq y '$minf) (alike1 y '((mtimes) -1 $inf)))
 689          (div '$%pi -2))
 690         ((and $%piargs
 691               ;; Recognize more special values
 692               (cond ((equal 1 y) (div '$%pi 4))
 693                     ((equal -1 y) (div '$%pi -4))
 694                     ;; sqrt(3)
 695                     ((alike1 y '((mexpt) 3 ((rat) 1 2)))
 696                      (div '$%pi 3))
 697                     ;; -sqrt(3)
 698                     ((alike1 y '((mtimes) -1 ((mexpt) 3 ((rat) 1 2))))
 699                      (div '$%pi -3))
 700                     ;; 1/sqrt(3)
 701                     ((alike1 y '((mexpt) 3 ((rat) -1 2)))
 702                      (div '$%pi 6))
 703                     ;; -1/sqrt(3)
 704                     ((alike1 y '((mtimes) -1 ((mexpt) 3 ((rat) -1 2))))
 705                      (div '$%pi -6))
 706                     ((alike1 y '((mplus) -1 ((mexpt) 2 ((rat) 1 2))))
 707                      (div '$%pi 8))
 708                     ((alike1 y '((mplus) 1 ((mexpt) 2 ((rat) 1 2))))
 709                      (mul 3 (div '$%pi 8))))))
 710         ((and $%iargs (multiplep y '$%i))
 711          ;; atan(%i*y) -> %i*atanh(y)
 712          (mul '$%i (take '(%atanh) (coeff y '$%i 1))))
 713         ((and (not (atom y)) (member (caar y) '(%cot %tan))
 714               (if ($constantp (cadr y))
 715                   (let ((y-val (mfuncall '$mod
 716                                          (if (eq (caar y) '%tan)
 717                                              (cadr y)
 718                                              (sub %pi//2 (cadr y)))
 719                                          '$%pi)))
 720                     (cond ((eq (mlsp y-val %pi//2) t) y-val)
 721                           ((eq (mlsp y-val '$%pi) t) (sub y-val '$%pi)))))))
 722         ((and (eq $triginverses '$all) (not (atom y))
 723               (if (eq (caar y) '%tan) (cadr y))))
 724         ((and (eq $triginverses t) (not (atom y)) (eq (caar y) '%tan)
 725               ;; Check if y in [-%pi/2, %pi/2]
 726               (if (and (member (csign (sub (cadr y) %pi//2)) '($nz $neg) :test #'eq)
 727                        (member (csign (add (cadr y) %pi//2)) '($pz $pos) :test #'eq))
 728                   (cadr y))))
 729         ($logarc (logarc '%atan y))
 730         ((apply-reflection-simp (mop form) y $trigsign))
 731         (t (give-up))))
 732
 733 (defun %piargs (x ratcoeff)
 734   (let (offset-result)
 735     (cond ((and (integerp (car x)) (integerp (cdr x))) 0)
 736           ((not (mevenp (car x)))
 737            (cond ((null ratcoeff) nil)
 738                  ((and (integerp (car x))
 739                        (setq offset-result (%piargs-offset (cdr x))))
 740                   (mul (power -1 (sub ratcoeff (cdr x)))
 741                        offset-result))))
 742           ((%piargs-offset (mmod (cdr x) 2))))))
 743
 744 ; simplifies sin(%pi * x) where x is between 0 and 1
 745 ; returns nil if can't simplify
 746 (defun %piargs-offset (x)
 747   (cond ((or (alike1 '((rat) 1 6) x) (alike1 '((rat) 5 6) x)) 1//2)
 748         ((or (alike1 '((rat) 1 4) x) (alike1 '((rat) 3 4) x)) (div (power 2 1//2) 2))
 749         ((or (alike1 '((rat) 1 3) x) (alike1 '((rat) 2 3) x)) (div (power 3 1//2) 2))
 750         ((alike1 1//2 x) 1)
 751         ((or (alike1 '((rat) 7 6) x) (alike1 '((rat) 11 6) x)) -1//2)
 752         ((or (alike1 '((rat) 4 3) x) (alike1 '((rat) 5 3) x)) (div (power 3 1//2) -2))
 753         ((or (alike1 '((rat) 5 4) x) (alike1 '((rat) 7 4) x)) (mul -1//2 (power 2 1//2)))
 754         ((alike1 '((rat) 3 2) x) -1)))
 755
 756 ;; identifies integer part of form
 757 ;; returns (X . Y) if form can be written as X*some_integer + Y
 758 ;; returns nil otherwise
 759 (defun linearize (form)
 760   (cond ((integerp form) (cons 0 form))
 761         ((numberp form) nil)
 762         ((atom form)
 763          (let (dum)
 764            (cond ((setq dum (evod form))
 765                   (if (eq '$even dum) '(2 . 0) '(2 . 1)))
 766                  ((maxima-integerp form) '(1 . 0)))))
 767         ((eq 'rat (caar form)) (cons 0 form))
 768         ((eq 'mplus (caar form)) (lin-mplus form))
 769         ((eq 'mtimes (caar form)) (lin-mtimes form))
 770         ((eq 'mexpt (caar form)) (lin-mexpt form))))
 771
 772 (defun lin-mplus (form)
 773   (do ((tl (cdr form) (cdr tl)) (dummy) (coeff 0) (zl-rem 0))
 774       ((null tl) (cons coeff (mmod zl-rem coeff)))
 775     (setq dummy (linearize (car tl)))
 776     (if (null dummy) (return nil)
 777         (setq coeff (rgcd (car dummy) coeff) zl-rem (rplus (cdr dummy) zl-rem)))))
 778
 779 (defun lin-mtimes (form)
 780   (do ((fl (cdr form) (cdr fl)) (dummy) (coeff 0) (zl-rem 1))
 781       ((null fl) (cons coeff (mmod zl-rem coeff)))
 782     (setq dummy (linearize (car fl)))
 783     (cond ((null dummy) (return nil))
 784           (t (setq coeff (rgcd (rtimes coeff (car dummy))
 785                                (rgcd (rtimes coeff (cdr dummy)) (rtimes zl-rem (car dummy))))
 786                    zl-rem (rtimes (cdr dummy) zl-rem))))))
 787
 788 (defun lin-mexpt (form)
 789   (prog (dummy)
 790      (cond ((and (integerp (caddr form)) (not (minusp (caddr form)))
 791                  (not (null (setq dummy (linearize (cadr form))))))
 792             (return (cons (car dummy) (mmod (cdr dummy) (caddr form))))))))
 793
 794 (defun rgcd (x y)
 795   (cond ((integerp x)
 796          (cond ((integerp y) (gcd x y))
 797                (t (list '(rat) (gcd x (cadr y)) (caddr y)))))
 798         ((integerp y) (list '(rat) (gcd (cadr x) y) (caddr x)))
 799         (t (list '(rat) (gcd (cadr x) (cadr y)) (lcm (caddr x) (caddr y))))))
 800
 801 (defun maxima-reduce (x y)
 802   (prog (gcd)
 803      (setq gcd (gcd x y) x (truncate x gcd) y (truncate y gcd))
 804      (if (minusp y) (setq x (- x) y (- y)))
 805      (return (if (eql y 1) x (list '(rat simp) x y)))))
 806
 807 ;; The following four functions are generated in code by TRANSL. - JPG 2/1/81
 808
 809 (defun rplus (x y) (addk x y))
 810
 811 (defun rdifference (x y) (addk x (timesk -1 y)))
 812
 813 (defun rtimes (x y) (timesk x y))
 814
 815 (defun rremainder (x y)
 816   (cond ((equal 0 y) (dbz-err))
 817         ((integerp x)
 818          (cond ((integerp y) (maxima-reduce x y))
 819                (t (maxima-reduce (* x (caddr y)) (cadr y)))))
 820         ((integerp y) (maxima-reduce (cadr x) (* (caddr x) y)))
 821         (t (maxima-reduce (* (cadr x) (caddr y)) (* (caddr x) (cadr y))))))
 822
 823 (defmfun $exponentialize (exp)
 824   (let ($demoivre)
 825     (cond ((atom exp) exp)
 826           ((trigp (caar exp))
 827            (exponentialize (caar exp) ($exponentialize (cadr exp))))
 828           (t (recur-apply #'$exponentialize exp)))))
 829
 830 (defun exponentialize (op arg)
 831   (cond ((eq '%sin op)
 832          (div (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))
 833               (mul 2 '$%i)))
 834         ((eq '%cos op)
 835          (div (add (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg))) 2))
 836         ((eq '%tan op)
 837          (div (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))
 838               (mul '$%i (add (power '$%e (mul '$%i arg))
 839                              (power '$%e (mul -1 '$%i arg))))))
 840         ((eq '%cot op)
 841          (div (mul '$%i (add (power '$%e (mul '$%i arg))
 842                              (power '$%e (mul -1 '$%i arg))))
 843               (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))))
 844         ((eq '%csc op)
 845          (div (mul 2 '$%i)
 846               (sub (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))))
 847         ((eq '%sec op)
 848          (div 2 (add (power '$%e (mul '$%i arg)) (power '$%e (mul -1 '$%i arg)))))
 849         ((eq '%sinh op)
 850          (div (sub (power '$%e arg) (power '$%e (neg arg))) 2))
 851         ((eq '%cosh op)
 852          (div (add (power '$%e arg) (power '$%e (mul -1 arg))) 2))
 853         ((eq '%tanh op)
 854          (div (sub (power '$%e arg) (power '$%e (neg arg)))
 855               (add (power '$%e arg) (power '$%e (mul -1 arg)))))
 856         ((eq '%coth op)
 857          (div (add (power '$%e arg) (power '$%e (mul -1 arg)))
 858               (sub (power '$%e arg) (power '$%e (neg arg)))))
 859         ((eq '%csch op)
 860          (div 2 (sub (power '$%e arg) (power '$%e (neg arg)))))
 861         ((eq '%sech op)
 862          (div 2 (add (power '$%e arg) (power '$%e (mul -1 arg)))))))
 863
 864 (defun coefficient (exp var pow)
 865   (coeff exp var pow))
 866
 867 (defun mmod (x mod)
 868   (cond ((and (integerp x) (integerp mod))
 869          (if (minusp (if (zerop mod) x (setq x (- x (* mod (truncate x mod))))))
 870              (+ x mod)
 871              x))
 872         ((and ($ratnump x) ($ratnump mod))
 873          (let
 874              ((d (lcm ($denom x) ($denom mod))))
 875            (setq x (mul* d x))
 876            (setq mod (mul* d mod))
 877            (div (mod x mod) d)))
 878         (t nil)))
 879
 880 (defun multiplep (exp var)
 881   (and (not (zerop1 exp)) (zerop1 (sub exp (mul var (coeff exp var 1))))))
 882
 883 (defun linearp (exp var)
 884   (and (setq exp (islinear exp var)) (not (equal (car exp) 0))))
 885
 886 (defun mminusp (x)
 887   (= -1 (signum1 x)))
 888
 889 (defun mminusp* (x)
 890   (let (sign)
 891     (setq sign (csign x))
 892     (or (member sign '($neg $nz) :test #'eq)
 893         (and (mminusp x) (not (member sign '($pos $pz) :test #'eq))))))
 894
 895 ;; This should give more information somehow.
 896
 897 (defun dbz-err ()
 898   (cond ((not errorsw) (merror (intl:gettext "Division by zero attempted.")))
 899         (t (throw 'errorsw t))))
 900
 901 (defun dbz-err1 (func)
 902   (cond ((not errorsw) (merror (intl:gettext "~A: division by zero attempted.") func))
 903         (t (throw 'errorsw t))))