src/trigi.lisp

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