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