1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module trigi
)
15 (load-macsyma-macros mrgmac
)
17 (declare-top (special errorsw $demoivre
1//2 -
1//2))
21 (defmvar $triginverses t
)
22 (defmvar $trigexpand nil
)
23 (defmvar $trigexpandplus t
)
24 (defmvar $trigexpandtimes t
)
26 (defmvar $exponentialize nil
)
28 (defmvar $halfangles nil
)
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
)
35 ((mexpt simp
) 3 ((rat simp
) 1 2))))
36 (defvar -sqrt3
//2 '((mtimes simp
)
38 ((mexpt simp
) 3 ((rat simp
) 1 2))))
40 ;;; Arithmetic utilities.
43 (power (sub 1 (power x
2)) 1//2))
46 (power (add 1 (power x
2)) 1//2))
49 (power (add (power x
2) -
1) 1//2))
52 (power (add (power x
2) (power y
2)) 1//2))
55 (member func
'(%sin %cos %tan %csc %sec %cot %sinh %cosh %tanh %csch %sech %coth
)
59 (member func
'(%asin %acos %atan %acsc %asec %acot %asinh %acosh %atanh %acsch %asech %acoth
)
62 ;;; The trigonometric functions distribute of lists, matrices and equations.
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
)))
70 (defun domain-error (x f
)
71 (merror (intl:gettext
"~A: argument ~:M isn't in the domain of ~A.") f
(complexify x
) f
))
73 ;; Build hash tables '*flonum-op*' and '*big-float-op*' that map Maxima
74 ;; function names to their corresponding Lisp functions.
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
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.")
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.
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.
100 ;; Let's rewrite it, where abs(x)>1
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))]
110 ;; Now, look at log part. If x > 0, we have
112 ;; log(x*(1+sqrt(1-1/x^2)))
114 ;; which is just fine. For x < 0, we have
116 ;; log(abs(x)*(1-sqrt(1-1/x^2))).
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))
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)))
132 ;; asin(x) = -%pi/2+%i*log(abs(x)*(1+sqrt(1-1/x^2)))
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
140 (let* ((absx (abs x
))
142 (result (complex (/ #.
(float pi
) 2)
144 (1+ (* (sqrt (+ 1 recip
))
145 (sqrt (- 1 recip
))))))))))
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
))
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))
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)
173 ;; Fill the hash table.
174 (macrolet ((frob (mfun dfun
) `(setf (gethash ',mfun
*flonum-op
*) ,dfun
)))
184 (frob %sec
#'(lambda (x)
185 (let ((y (ignore-errors (/ 1 (cl:cos x
)))))
186 (if y y
(domain-error x
'sec
)))))
188 (frob %csc
#'(lambda (x)
189 (let ((y (ignore-errors (/ 1 (cl:sin x
)))))
190 (if y y
(domain-error x
'csc
)))))
192 (frob %cot
#'(lambda (x)
193 (let ((y (ignore-errors (/ 1 (cl:tan x
)))))
194 (if y y
(domain-error x
'cot
)))))
196 (frob %acos
#'maxima-branch-acos
)
197 (frob %asin
#'maxima-branch-asin
)
199 (frob %atan
#'cl
:atan
)
201 (frob %asec
#'(lambda (x)
202 (let ((y (ignore-errors (maxima-branch-acos (/ 1 x
)))))
203 (if y y
(domain-error x
'asec
)))))
205 (frob %acsc
#'(lambda (x)
206 (let ((y (ignore-errors (maxima-branch-asin (/ 1 x
)))))
207 (if y y
(domain-error x
'acsc
)))))
209 (frob %acot
#'(lambda (x)
210 (let ((y (ignore-errors (maxima-branch-acot x
))))
211 (if y y
(domain-error x
'acot
)))))
213 (frob %cosh
#'cl
:cosh
)
214 (frob %sinh
#'cl
:sinh
)
215 (frob %tanh
#'cl
:tanh
)
217 (frob %sech
#'(lambda (x)
218 (let ((y (ignore-errors (/ 1 (cl:cosh x
)))))
219 (if y y
(domain-error x
'sech
)))))
221 (frob %csch
#'(lambda (x)
222 (let ((y (ignore-errors (/ 1 (cl:sinh x
)))))
223 (if y y
(domain-error x
'csch
)))))
225 (frob %coth
#'(lambda (x)
226 (let ((y (ignore-errors (/ 1 (cl:tanh x
)))))
227 (if y y
(domain-error x
'coth
)))))
229 (frob %acosh
#'cl
:acosh
)
230 (frob %asinh
#'cl
:asinh
)
232 (frob %atanh
#'maxima-branch-atanh
)
234 (frob %asech
#'(lambda (x)
235 (let ((y (ignore-errors (cl:acosh
(/ 1 x
)))))
236 (if y y
(domain-error x
'asech
)))))
238 (frob %acsch
#'(lambda (x)
239 (let ((y (ignore-errors (cl:asinh
(/ 1 x
)))))
240 (if y y
(domain-error x
'acsch
)))))
242 (frob %acoth
#'(lambda (x)
243 (let ((y (ignore-errors (maxima-branch-atanh (/ 1 x
)))))
244 (if y y
(domain-error x
'acoth
)))))
248 (frob mexpt
#'cl
:expt
)
249 (frob %sqrt
#'cl
:sqrt
)
250 (frob %log
#'(lambda (x)
251 (let ((y (ignore-errors (cl:log x
))))
252 (if y y
(domain-error x
'log
)))))
254 (frob %plog
#'(lambda (x)
255 (let ((y (ignore-errors (cl:log x
))))
256 (if y y
(domain-error x
'log
)))))
258 (frob $conjugate
#'cl
:conjugate
)
259 (frob $floor
#'cl
:ffloor
)
260 (frob $ceiling
#'cl
:fceiling
)
261 (frob $realpart
#'cl
:realpart
)
262 (frob $imagpart
#'cl
:imagpart
)
265 (frob %signum
#'cl
:signum
)
266 (frob $atan2
#'cl
:atan
))
268 (macrolet ((frob (mfun dfun
) `(setf (gethash ',mfun
*big-float-op
*) ,dfun
)))
269 ;; All big-float implementation functions MUST support a required x
270 ;; arg and an optional y arg for the real and imaginary parts. The
271 ;; imaginary part does not have to be given.
272 (frob %asin
#'big-float-asin
)
273 (frob %sinh
#'big-float-sinh
)
274 (frob %asinh
#'big-float-asinh
)
275 (frob %tanh
#'big-float-tanh
)
276 (frob %atanh
#'big-float-atanh
)
277 (frob %acos
'big-float-acos
)
278 (frob %log
'big-float-log
)
279 (frob %sqrt
'big-float-sqrt
))
281 ;; Here is a general scheme for defining and applying reflection rules. A
282 ;; reflection rule is something like f(-x) --> f(x), or f(-x) --> %pi - f(x).
284 ;; We define functions for the two most common reflection rules; these
285 ;; are the odd function rule (f(-x) --> -f(x)) and the even function rule
286 ;; (f(-x) --> f(x)). A reflection rule takes two arguments (the operator and
289 (defun odd-function-reflect (op x
)
290 (neg (take (list op
) (neg x
))))
292 (defun even-function-reflect (op x
)
293 (take (list op
) (neg x
)))
295 ;; Put the reflection rule on the property list of the exponential-like
298 (setf (get '%cos
'reflection-rule
) #'even-function-reflect
)
299 (setf (get '%sin
'reflection-rule
) #'odd-function-reflect
)
300 (setf (get '%tan
'reflection-rule
) #'odd-function-reflect
)
301 (setf (get '%sec
'reflection-rule
) #'even-function-reflect
)
302 (setf (get '%csc
'reflection-rule
) #'odd-function-reflect
)
303 (setf (get '%cot
'reflection-rule
) #'odd-function-reflect
)
305 ;; See A&S 4.4.14--4.4.19
307 (setf (get '%acos
'reflection-rule
) #'(lambda (op x
) (sub '$%pi
(take (list op
) (neg x
)))))
308 (setf (get '%asin
'reflection-rule
) #'odd-function-reflect
)
309 (setf (get '%atan
'reflection-rule
) #'odd-function-reflect
)
310 (setf (get '%asec
'reflection-rule
) #'(lambda (op x
) (sub '$%pi
(take (list op
) (neg x
)))))
311 (setf (get '%acsc
'reflection-rule
) #'odd-function-reflect
)
312 (setf (get '%acot
'reflection-rule
) #'odd-function-reflect
)
314 (setf (get '%cosh
'reflection-rule
) #'even-function-reflect
)
315 (setf (get '%sinh
'reflection-rule
) #'odd-function-reflect
)
316 (setf (get '%tanh
'reflection-rule
) #'odd-function-reflect
)
317 (setf (get '%sech
'reflection-rule
) #'even-function-reflect
)
318 (setf (get '%csch
'reflection-rule
) #'odd-function-reflect
)
319 (setf (get '%coth
'reflection-rule
) #'odd-function-reflect
)
321 (setf (get '%asinh
'reflection-rule
) #'odd-function-reflect
)
322 (setf (get '%atanh
'reflection-rule
) #'odd-function-reflect
)
323 (setf (get '%asech
'reflection-rule
) #'even-function-reflect
)
324 (setf (get '%acsch
'reflection-rule
) #'odd-function-reflect
)
325 (setf (get '%acoth
'reflection-rule
) #'odd-function-reflect
)
327 ;; When b is nil, do not apply the reflection rule. For trigonometric like
328 ;; functions, b is $trigsign. This function uses 'great' to decide when to
329 ;; apply the rule. Another possibility is to apply the rule when (mminusp* x)
330 ;; evaluates to true. Maxima <= 5.9.3 uses this scheme; with this method, we have
331 ;; assume(z < 0), cos(z) --> cos(-z). I (Barton Willis) think this goofy.
333 ;; The function 'great' is non-transitive. I don't think this bug will cause
334 ;; trouble for this function. If there is an expression such that both
335 ;; (great (neg x) x) and (great x (neg x)) evaluate to true, this function
336 ;; could cause an infinite loop. I could protect against this possibility with
337 ;; (and b f (great (neg x) x) (not (great x (neg x))).
339 (defun apply-reflection-simp (op x
&optional
(b t
))
340 (let ((f (get op
'reflection-rule
)))
341 (if (and b f
(great (neg x
) x
)) (funcall f op x
) nil
)))
343 (defun taylorize (op x
)
345 (mfuncall '$apply
'$taylor
`((mlist) ((,op
) ,($ratdisrep x
)) ,@(cdr ($taylorinfo x
)))) nil
))
347 (defun float-or-rational-p (x)
348 (or (floatp x
) ($ratnump x
)))
350 (defun bigfloat-or-number-p (x)
351 (or ($bfloatp x
) (numberp x
) ($ratnump x
)))
353 ;; When z is a Maxima complex float or when 'numer' is true and z is a
354 ;; Maxima complex number, evaluate (op z) by applying the mapping from
355 ;; the Maxima operator 'op' to the operator in the hash table
356 ;; '*flonum-op*'. When z isn't a Maxima complex number, return
359 (defun flonum-eval (op z
)
360 (let ((op (gethash op
*flonum-op
*)))
362 (multiple-value-bind (bool R I
)
363 (complex-number-p z
#'float-or-rational-p
)
364 (when (and bool
(or $numer
(floatp R
) (floatp I
)))
367 (complexify (funcall op
(if (zerop I
) R
(complex R I
)))))))))
369 ;; For now, big float evaluation of trig-like functions for complex
370 ;; big floats uses rectform. I suspect that for some functions (not
371 ;; all of them) rectform generates expressions that are poorly suited
372 ;; for numerical evaluation. For better accuracy, these functions
373 ;; (maybe acosh, for one) may need to be special cased. If they are
374 ;; special-cased, the *big-float-op* hash table contains the special
377 (defun big-float-eval (op z
)
378 (when (complex-number-p z
'bigfloat-or-number-p
)
379 (let ((x ($realpart z
))
381 (bop (gethash op
*big-float-op
*)))
382 ;; If bop is non-NIL, we want to try that first. If bop
383 ;; declines (by returning NIL), we silently give up and use the
385 (cond ((and ($bfloatp x
) (like 0 y
))
386 (or (and bop
(funcall bop x
))
387 ($bfloat
`((,op simp
) ,x
))))
388 ((or ($bfloatp x
) ($bfloatp y
))
389 (or (and bop
(funcall bop
($bfloat x
) ($bfloat y
)))
390 (let ((z (add ($bfloat x
) (mul '$%i
($bfloat y
)))))
391 (setq z
($rectform
`((,op simp
) ,z
)))
394 ;; For complex big float evaluation, it's important to check the
395 ;; simp flag -- otherwise Maxima can get stuck in an infinite loop:
396 ;; asin(1.23b0 + %i * 4.56b0) ---> (simp-%asin ((%asin) ...) -->
397 ;; (big-float-eval ((%asin) ...) --> (risplit ((%asin simp) ...) -->
398 ;; (simp-%asin ((%asin simp) ...). If the simp flag is ignored, we've
401 (def-simplifier sin
(y)
403 (cond ((flonum-eval (mop form
) y
))
404 ((and (not (member 'simp
(car form
))) (big-float-eval (mop form
) y
)))
405 ((taylorize (mop form
) (second form
)))
406 ((and $%piargs
(cond ((zerop1 y
) 0)
407 ((has-const-or-int-term y
'$%pi
) (%piargs-sin
/cos y
)))))
408 ((and $%iargs
(multiplep y
'$%i
)) (mul '$%i
(ftake* '%sinh
(coeff y
'$%i
1))))
409 ((and $triginverses
(not (atom y
))
410 (cond ((eq '%asin
(setq z
(caar y
))) (cadr y
))
411 ((eq '%acos z
) (sqrt1-x^
2 (cadr y
)))
412 ((eq '%atan z
) (div (cadr y
) (sqrt1+x^
2 (cadr y
))))
413 ((eq '%acot z
) (div 1 (sqrt1+x^
2 (cadr y
))))
414 ((eq '%asec z
) (div (sqrtx^
2-
1 (cadr y
)) (cadr y
)))
415 ((eq '%acsc z
) (div 1 (cadr y
)))
416 ((eq '$atan2 z
) (div (cadr y
) (sq-sumsq (cadr y
) (caddr y
)))))))
417 ((and $trigexpand
(trigexpand '%sin y
)))
418 ($exponentialize
(exponentialize '%sin y
))
419 ((and $halfangles
(halfangle '%sin y
)))
420 ((apply-reflection-simp (mop form
) y $trigsign
))
421 ;((and $trigsign (mminusp* y)) (neg (ftake* '%sin (neg y))))
424 (def-simplifier cos
(y)
426 (cond ((flonum-eval (mop form
) y
))
427 ((and (not (member 'simp
(car form
))) (big-float-eval (mop form
) y
)))
428 ((taylorize (mop form
) (second form
)))
429 ((and $%piargs
(cond ((zerop1 y
) 1)
430 ((has-const-or-int-term y
'$%pi
)
431 (%piargs-sin
/cos
(add %pi
//2 y
))))))
432 ((and $%iargs
(multiplep y
'$%i
)) (ftake* '%cosh
(coeff y
'$%i
1)))
433 ((and $triginverses
(not (atom y
))
434 (cond ((eq '%acos
(setq z
(caar y
))) (cadr y
))
435 ((eq '%asin z
) (sqrt1-x^
2 (cadr y
)))
436 ((eq '%atan z
) (div 1 (sqrt1+x^
2 (cadr y
))))
437 ((eq '%acot z
) (div (cadr y
) (sqrt1+x^
2 (cadr y
))))
438 ((eq '%asec z
) (div 1 (cadr y
)))
439 ((eq '%acsc z
) (div (sqrtx^
2-
1 (cadr y
)) (cadr y
)))
440 ((eq '$atan2 z
) (div (caddr y
) (sq-sumsq (cadr y
) (caddr y
)))))))
441 ((and $trigexpand
(trigexpand '%cos y
)))
442 ($exponentialize
(exponentialize '%cos y
))
443 ((and $halfangles
(halfangle '%cos y
)))
444 ((apply-reflection-simp (mop form
) y $trigsign
))
445 ;((and $trigsign (mminusp* y)) (ftake* '%cos (neg y)))
448 (defun %piargs-sin
/cos
(x)
449 (let ($float coeff ratcoeff zl-rem
)
450 (setq ratcoeff
(get-const-or-int-terms x
'$%pi
)
451 coeff
(linearize ratcoeff
)
452 zl-rem
(get-not-const-or-int-terms x
'$%pi
))
453 (cond ((zerop1 zl-rem
) (%piargs coeff ratcoeff
))
454 ((not (mevenp (car coeff
))) nil
)
455 ((equal 0 (setq x
(mmod (cdr coeff
) 2))) (ftake* '%sin zl-rem
))
456 ((equal 1 x
) (neg (ftake* '%sin zl-rem
)))
457 ((alike1 1//2 x
) (ftake* '%cos zl-rem
))
458 ((alike1 '((rat) 3 2) x
) (neg (ftake* '%cos zl-rem
))))))
461 (defun filter-sum (pred form simp-flag
)
462 "Takes form to be a sum and a sum of the summands for which pred is
463 true. Passes simp-flag through to addn if there is more than one
468 (when (funcall pred term
) (list term
))) (cdr form
))
470 (if (funcall pred form
) form
0)))
472 ;; collect terms of form A*var where A is a constant or integer.
473 ;; returns sum of all such A.
474 ;; does not expand form, so does not find constant term in (x+1)*var.
475 ;; thus we cannot simplify sin(2*%pi*(1+x)) => sin(2*%pi*x) unless
476 ;; the user calls expand. this could be extended to look a little
477 ;; more deeply into the expression, but we don't want to call expand
478 ;; in the core simplifier for reasons of speed and predictability.
479 (defun get-const-or-int-terms (form var
)
481 (filter-sum (lambda (term)
482 (let ((coeff (coeff term var
1)))
483 (and (not (zerop1 coeff
))
484 (or ($constantp coeff
)
485 (maxima-integerp coeff
)))))
490 ;; collect terms skipped by get-const-or-int-terms
491 (defun get-not-const-or-int-terms (form var
)
492 (filter-sum (lambda (term)
493 (let ((coeff (coeff term var
1)))
494 (not (and (not (zerop1 coeff
))
495 (or ($constantp coeff
)
496 (maxima-integerp coeff
))))))
500 (defun has-const-or-int-term (form var
)
501 "Tests whether form has at least some term of the form a*var where a
502 is constant or integer"
503 (not (zerop1 (get-const-or-int-terms form var
))))
505 (def-simplifier tan
(y)
507 (cond ((flonum-eval (mop form
) y
))
508 ((and (not (member 'simp
(car form
))) (big-float-eval (mop form
) y
)))
509 ((taylorize (mop form
) (second form
)))
510 ((and $%piargs
(cond ((zerop1 y
) 0)
511 ((has-const-or-int-term y
'$%pi
) (%piargs-tan
/cot y
)))))
512 ((and $%iargs
(multiplep y
'$%i
)) (mul '$%i
(ftake* '%tanh
(coeff y
'$%i
1))))
513 ((and $triginverses
(not (atom y
))
514 (cond ((eq '%atan
(setq z
(caar y
))) (cadr y
))
515 ((eq '%asin z
) (div (cadr y
) (sqrt1-x^
2 (cadr y
))))
516 ((eq '%acos z
) (div (sqrt1-x^
2 (cadr y
)) (cadr y
)))
517 ((eq '%acot z
) (div 1 (cadr y
)))
518 ((eq '%asec z
) (sqrtx^
2-
1 (cadr y
)))
519 ((eq '%acsc z
) (div 1 (sqrtx^
2-
1 (cadr y
))))
520 ((eq '$atan2 z
) (div (cadr y
) (caddr y
))))))
521 ((and $trigexpand
(trigexpand '%tan y
)))
522 ($exponentialize
(exponentialize '%tan y
))
523 ((and $halfangles
(halfangle '%tan y
)))
524 ((apply-reflection-simp (mop form
) y $trigsign
))
525 ;((and $trigsign (mminusp* y)) (neg (ftake* '%tan (neg y))))
528 (def-simplifier cot
(y)
530 (cond ((flonum-eval (mop form
) y
))
531 ((and (not (member 'simp
(car form
))) (big-float-eval (mop form
) y
)))
532 ((taylorize (mop form
) (second form
)))
533 ((and $%piargs
(cond ((zerop1 y
) (domain-error y
'cot
))
534 ((and (has-const-or-int-term y
'$%pi
)
535 (setq z
(%piargs-tan
/cot
(add %pi
//2 y
))))
537 ((and $%iargs
(multiplep y
'$%i
)) (mul -
1 '$%i
(ftake* '%coth
(coeff y
'$%i
1))))
538 ((and $triginverses
(not (atom y
))
539 (cond ((eq '%acot
(setq z
(caar y
))) (cadr y
))
540 ((eq '%asin z
) (div (sqrt1-x^
2 (cadr y
)) (cadr y
)))
541 ((eq '%acos z
) (div (cadr y
) (sqrt1-x^
2 (cadr y
))))
542 ((eq '%atan z
) (div 1 (cadr y
)))
543 ((eq '%asec z
) (div 1 (sqrtx^
2-
1 (cadr y
))))
544 ((eq '%acsc z
) (sqrtx^
2-
1 (cadr y
)))
545 ((eq '$atan2 z
) (div (caddr y
) (cadr y
))))))
546 ((and $trigexpand
(trigexpand '%cot y
)))
547 ($exponentialize
(exponentialize '%cot y
))
548 ((and $halfangles
(halfangle '%cot y
)))
549 ((apply-reflection-simp (mop form
) y $trigsign
))
550 ;((and $trigsign (mminusp* y)) (neg (ftake* '%cot (neg y))))
553 (defun %piargs-tan
/cot
(x)
554 "If x is of the form tan(u) where u has a nonzero constant linear
555 term in %pi, then %piargs-tan/cot returns a simplified version of x
556 without this constant term."
557 ;; Set coeff to be the coefficient of $%pi collecting terms with no
558 ;; other atoms, so given %pi(x+1/2), coeff = 1/2. Let zl-rem be the
559 ;; remainder (TODO: computing zl-rem could probably be prettier.)
560 (let* ((nice-terms (get-const-or-int-terms x
'$%pi
))
561 (coeff (linearize nice-terms
))
562 (zl-rem (get-not-const-or-int-terms x
'$%pi
))
566 ;; sin-of-coeff-pi and cos-of-coeff-pi are only non-nil if they
567 ;; are constants that %piargs-offset could compute, and we just
568 ;; checked that cos-of-coeff-pi was nonzero. Thus we can just
569 ;; return their quotient.
570 ((and (zerop1 zl-rem
)
571 (setq sin-of-coeff-pi
572 (%piargs coeff nil
)))
573 (setq cos-of-coeff-pi
574 (%piargs
(cons (car coeff
)
575 (rplus 1//2 (cdr coeff
))) nil
))
576 (cond ((zerop1 sin-of-coeff-pi
)
577 0) ;; tan(integer*%pi)
578 ((zerop1 cos-of-coeff-pi
)
579 (merror (intl:gettext
"tan: ~M isn't in the domain of tan.") x
))
581 (div sin-of-coeff-pi cos-of-coeff-pi
))))
583 ;; This expression sets x to the coeff of %pi (mod 1) as a side
584 ;; effect and then, if this is zero, returns tan of the
585 ;; rest, because tan has periodicity %pi.
586 ((zerop1 (setq x
(mmod (cdr coeff
) 1)))
587 (ftake* '%tan zl-rem
))
589 ;; Similarly, if x = 1/2 then return -cot(x).
591 (neg (ftake* '%cot zl-rem
))))))
593 (def-simplifier csc
(y)
595 (cond ((flonum-eval (mop form
) y
))
596 ((and (not (member 'simp
(car form
))) (big-float-eval (mop form
) y
)))
597 ((taylorize (mop form
) (second form
)))
598 ((and $%piargs
(cond ((zerop1 y
) (domain-error y
'csc
))
599 ((has-const-or-int-term y
'$%pi
) (%piargs-csc
/sec y
)))))
600 ((and $%iargs
(multiplep y
'$%i
)) (mul -
1 '$%i
(ftake* '%csch
(coeff y
'$%i
1))))
601 ((and $triginverses
(not (atom y
))
602 (cond ((eq '%acsc
(setq z
(caar y
))) (cadr y
))
603 ((eq '%asin z
) (div 1 (cadr y
)))
604 ((eq '%acos z
) (div 1 (sqrt1-x^
2 (cadr y
))))
605 ((eq '%atan z
) (div (sqrt1+x^
2 (cadr y
)) (cadr y
)))
606 ((eq '%acot z
) (sqrt1+x^
2 (cadr y
)))
607 ((eq '%asec z
) (div (cadr y
) (sqrtx^
2-
1 (cadr y
))))
608 ((eq '$atan2 z
) (div (sq-sumsq (cadr y
) (caddr y
)) (cadr y
))))))
609 ((and $trigexpand
(trigexpand '%csc y
)))
610 ($exponentialize
(exponentialize '%csc y
))
611 ((and $halfangles
(halfangle '%csc y
)))
612 ((apply-reflection-simp (mop form
) y $trigsign
))
613 ;((and $trigsign (mminusp* y)) (neg (ftake* '%csc (neg y))))
617 (def-simplifier sec
(y)
619 (cond ((flonum-eval (mop form
) y
))
620 ((and (not (member 'simp
(car form
))) (big-float-eval (mop form
) y
)))
621 ((taylorize (mop form
) (second form
)))
622 ((and $%piargs
(cond ((zerop1 y
) 1)
623 ((has-const-or-int-term y
'$%pi
) (%piargs-csc
/sec
(add %pi
//2 y
))))))
624 ((and $%iargs
(multiplep y
'$%i
)) (ftake* '%sech
(coeff y
'$%i
1)))
625 ((and $triginverses
(not (atom y
))
626 (cond ((eq '%asec
(setq z
(caar y
))) (cadr y
))
627 ((eq '%asin z
) (div 1 (sqrt1-x^
2 (cadr y
))))
628 ((eq '%acos z
) (div 1 (cadr y
)))
629 ((eq '%atan z
) (sqrt1+x^
2 (cadr y
)))
630 ((eq '%acot z
) (div (sqrt1+x^
2 (cadr y
)) (cadr y
)))
631 ((eq '%acsc z
) (div (cadr y
) (sqrtx^
2-
1 (cadr y
))))
632 ((eq '$atan2 z
) (div (sq-sumsq (cadr y
) (caddr y
)) (caddr y
))))))
633 ((and $trigexpand
(trigexpand '%sec y
)))
634 ($exponentialize
(exponentialize '%sec y
))
635 ((and $halfangles
(halfangle '%sec y
)))
636 ((apply-reflection-simp (mop form
) y $trigsign
))
637 ;((and $trigsign (mminusp* y)) (ftake* '%sec (neg y)))
641 (defun %piargs-csc
/sec
(x)
642 (prog ($float coeff ratcoeff zl-rem
)
643 (setq ratcoeff
(get-const-or-int-terms x
'$%pi
)
644 coeff
(linearize ratcoeff
)
645 zl-rem
(get-not-const-or-int-terms x
'$%pi
))
646 (return (cond ((and (zerop1 zl-rem
) (setq zl-rem
(%piargs coeff nil
))) (div 1 zl-rem
))
647 ((not (mevenp (car coeff
))) nil
)
648 ((equal 0 (setq x
(mmod (cdr coeff
) 2))) (ftake* '%csc zl-rem
))
649 ((equal 1 x
) (neg (ftake* '%csc zl-rem
)))
650 ((alike1 1//2 x
) (ftake* '%sec zl-rem
))
651 ((alike1 '((rat) 3 2) x
) (neg (ftake* '%sec zl-rem
)))))))
653 (def-simplifier atan
(y)
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 ;; Simplification for special values
659 ((or (eq y
'$inf
) (alike1 y
'((mtimes) -
1 $minf
)))
661 ((or (eq y
'$minf
) (alike1 y
'((mtimes) -
1 $inf
)))
664 ;; Recognize more special values
665 (cond ((equal 1 y
) (div '$%pi
4))
666 ((equal -
1 y
) (div '$%pi -
4))
668 ((alike1 y
'((mexpt) 3 ((rat) 1 2)))
671 ((alike1 y
'((mtimes) -
1 ((mexpt) 3 ((rat) 1 2))))
674 ((alike1 y
'((mexpt) 3 ((rat) -
1 2)))
677 ((alike1 y
'((mtimes) -
1 ((mexpt) 3 ((rat) -
1 2))))
679 ((alike1 y
'((mplus) -
1 ((mexpt) 2 ((rat) 1 2))))
681 ((alike1 y
'((mplus) 1 ((mexpt) 2 ((rat) 1 2))))
682 (mul 3 (div '$%pi
8))))))
683 ((and $%iargs
(multiplep y
'$%i
))
684 ;; atan(%i*y) -> %i*atanh(y)
685 (mul '$%i
(take '(%atanh
) (coeff y
'$%i
1))))
686 ((and (not (atom y
)) (member (caar y
) '(%cot %tan
))
687 (if ($constantp
(cadr y
))
688 (let ((y-val (mfuncall '$mod
689 (if (eq (caar y
) '%tan
)
691 (sub %pi
//2 (cadr y
)))
693 (cond ((eq (mlsp y-val %pi
//2) t
) y-val
)
694 ((eq (mlsp y-val
'$%pi
) t
) (sub y-val
'$%pi
)))))))
695 ((and (eq $triginverses
'$all
) (not (atom y
))
696 (if (eq (caar y
) '%tan
) (cadr y
))))
697 ((and (eq $triginverses t
) (not (atom y
)) (eq (caar y
) '%tan
)
698 ;; Check if y in [-%pi/2, %pi/2]
699 (if (and (member (csign (sub (cadr y
) %pi
//2)) '($nz $neg
) :test
#'eq
)
700 (member (csign (add (cadr y
) %pi
//2)) '($pz $pos
) :test
#'eq
))
702 ($logarc
(logarc '%atan y
))
703 ((apply-reflection-simp (mop form
) y $trigsign
))
706 (defun %piargs
(x ratcoeff
)
708 (cond ((and (integerp (car x
)) (integerp (cdr x
))) 0)
709 ((not (mevenp (car x
)))
710 (cond ((null ratcoeff
) nil
)
711 ((and (integerp (car x
))
712 (setq offset-result
(%piargs-offset
(cdr x
))))
713 (mul (power -
1 (sub ratcoeff
(cdr x
)))
715 ((%piargs-offset
(mmod (cdr x
) 2))))))
717 ; simplifies sin(%pi * x) where x is between 0 and 1
718 ; returns nil if can't simplify
719 (defun %piargs-offset
(x)
720 (cond ((or (alike1 '((rat) 1 6) x
) (alike1 '((rat) 5 6) x
)) 1//2)
721 ((or (alike1 '((rat) 1 4) x
) (alike1 '((rat) 3 4) x
)) (div (power 2 1//2) 2))
722 ((or (alike1 '((rat) 1 3) x
) (alike1 '((rat) 2 3) x
)) (div (power 3 1//2) 2))
724 ((or (alike1 '((rat) 7 6) x
) (alike1 '((rat) 11 6) x
)) -
1//2)
725 ((or (alike1 '((rat) 4 3) x
) (alike1 '((rat) 5 3) x
)) (div (power 3 1//2) -
2))
726 ((or (alike1 '((rat) 5 4) x
) (alike1 '((rat) 7 4) x
)) (mul -
1//2 (power 2 1//2)))
727 ((alike1 '((rat) 3 2) x
) -
1)))
729 ;; identifies integer part of form
730 ;; returns (X . Y) if form can be written as X*some_integer + Y
731 ;; returns nil otherwise
732 (defun linearize (form)
733 (cond ((integerp form
) (cons 0 form
))
737 (cond ((setq dum
(evod form
))
738 (if (eq '$even dum
) '(2 .
0) '(2 .
1)))
739 ((maxima-integerp form
) '(1 .
0)))))
740 ((eq 'rat
(caar form
)) (cons 0 form
))
741 ((eq 'mplus
(caar form
)) (lin-mplus form
))
742 ((eq 'mtimes
(caar form
)) (lin-mtimes form
))
743 ((eq 'mexpt
(caar form
)) (lin-mexpt form
))))
745 (defun lin-mplus (form)
746 (do ((tl (cdr form
) (cdr tl
)) (dummy) (coeff 0) (zl-rem 0))
747 ((null tl
) (cons coeff
(mmod zl-rem coeff
)))
748 (setq dummy
(linearize (car tl
)))
749 (if (null dummy
) (return nil
)
750 (setq coeff
(rgcd (car dummy
) coeff
) zl-rem
(rplus (cdr dummy
) zl-rem
)))))
752 (defun lin-mtimes (form)
753 (do ((fl (cdr form
) (cdr fl
)) (dummy) (coeff 0) (zl-rem 1))
754 ((null fl
) (cons coeff
(mmod zl-rem coeff
)))
755 (setq dummy
(linearize (car fl
)))
756 (cond ((null dummy
) (return nil
))
757 (t (setq coeff
(rgcd (rtimes coeff
(car dummy
))
758 (rgcd (rtimes coeff
(cdr dummy
)) (rtimes zl-rem
(car dummy
))))
759 zl-rem
(rtimes (cdr dummy
) zl-rem
))))))
761 (defun lin-mexpt (form)
763 (cond ((and (integerp (caddr form
)) (not (minusp (caddr form
)))
764 (not (null (setq dummy
(linearize (cadr form
))))))
765 (return (cons (car dummy
) (mmod (cdr dummy
) (caddr form
))))))))
769 (cond ((integerp y
) (gcd x y
))
770 (t (list '(rat) (gcd x
(cadr y
)) (caddr y
)))))
771 ((integerp y
) (list '(rat) (gcd (cadr x
) y
) (caddr x
)))
772 (t (list '(rat) (gcd (cadr x
) (cadr y
)) (lcm (caddr x
) (caddr y
))))))
774 (defun maxima-reduce (x y
)
776 (setq gcd
(gcd x y
) x
(truncate x gcd
) y
(truncate y gcd
))
777 (if (minusp y
) (setq x
(- x
) y
(- y
)))
778 (return (if (eql y
1) x
(list '(rat simp
) x y
)))))
780 ;; The following four functions are generated in code by TRANSL. - JPG 2/1/81
782 (defun rplus (x y
) (addk x y
))
784 (defun rdifference (x y
) (addk x
(timesk -
1 y
)))
786 (defun rtimes (x y
) (timesk x y
))
788 (defun rremainder (x y
)
789 (cond ((equal 0 y
) (dbz-err))
791 (cond ((integerp y
) (maxima-reduce x y
))
792 (t (maxima-reduce (* x
(caddr y
)) (cadr y
)))))
793 ((integerp y
) (maxima-reduce (cadr x
) (* (caddr x
) y
)))
794 (t (maxima-reduce (* (cadr x
) (caddr y
)) (* (caddr x
) (cadr y
))))))
796 (defmfun $exponentialize
(exp)
798 (cond ((atom exp
) exp
)
800 (exponentialize (caar exp
) ($exponentialize
(cadr exp
))))
801 (t (recur-apply #'$exponentialize exp
)))))
803 (defun exponentialize (op arg
)
805 (div (sub (power '$%e
(mul '$%i arg
)) (power '$%e
(mul -
1 '$%i arg
)))
808 (div (add (power '$%e
(mul '$%i arg
)) (power '$%e
(mul -
1 '$%i arg
))) 2))
810 (div (sub (power '$%e
(mul '$%i arg
)) (power '$%e
(mul -
1 '$%i arg
)))
811 (mul '$%i
(add (power '$%e
(mul '$%i arg
))
812 (power '$%e
(mul -
1 '$%i arg
))))))
814 (div (mul '$%i
(add (power '$%e
(mul '$%i arg
))
815 (power '$%e
(mul -
1 '$%i arg
))))
816 (sub (power '$%e
(mul '$%i arg
)) (power '$%e
(mul -
1 '$%i arg
)))))
819 (sub (power '$%e
(mul '$%i arg
)) (power '$%e
(mul -
1 '$%i arg
)))))
821 (div 2 (add (power '$%e
(mul '$%i arg
)) (power '$%e
(mul -
1 '$%i arg
)))))
823 (div (sub (power '$%e arg
) (power '$%e
(neg arg
))) 2))
825 (div (add (power '$%e arg
) (power '$%e
(mul -
1 arg
))) 2))
827 (div (sub (power '$%e arg
) (power '$%e
(neg arg
)))
828 (add (power '$%e arg
) (power '$%e
(mul -
1 arg
)))))
830 (div (add (power '$%e arg
) (power '$%e
(mul -
1 arg
)))
831 (sub (power '$%e arg
) (power '$%e
(neg arg
)))))
833 (div 2 (sub (power '$%e arg
) (power '$%e
(neg arg
)))))
835 (div 2 (add (power '$%e arg
) (power '$%e
(mul -
1 arg
)))))))
837 (defun coefficient (exp var pow
)
841 (cond ((and (integerp x
) (integerp mod
))
842 (if (minusp (if (zerop mod
) x
(setq x
(- x
(* mod
(truncate x mod
))))))
845 ((and ($ratnump x
) ($ratnump mod
))
847 ((d (lcm ($denom x
) ($denom mod
))))
849 (setq mod
(mul* d mod
))
850 (div (mod x mod
) d
)))
853 (defun multiplep (exp var
)
854 (and (not (zerop1 exp
)) (zerop1 (sub exp
(mul var
(coeff exp var
1))))))
856 (defun linearp (exp var
)
857 (and (setq exp
(islinear exp var
)) (not (equal (car exp
) 0))))
864 (setq sign
(csign x
))
865 (or (member sign
'($neg $nz
) :test
#'eq
)
866 (and (mminusp x
) (not (member sign
'($pos $pz
) :test
#'eq
))))))
868 ;; This should give more information somehow.
871 (cond ((not errorsw
) (merror (intl:gettext
"Division by zero attempted.")))
872 (t (throw 'errorsw t
))))
874 (defun dbz-err1 (func)
875 (cond ((not errorsw
) (merror (intl:gettext
"~A: division by zero attempted.") func
))
876 (t (throw 'errorsw t
))))