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