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 1981 Massachusetts Institute of Technology         ;;
9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
10
11(in-package :maxima)
12
13(macsyma-module laplac)
14
15(declare-top (special var $savefactors
16		      checkfactors $ratfac $keepfloat *nounl* *nounsflag*
17                      errcatch $errormsg))
18
19;;; The properties NOUN and VERB give correct linear display
20
21(defprop $laplace %laplace verb)
22(defprop %laplace $laplace noun)
23
24(defprop $ilt %ilt verb)
25(defprop %ilt $ilt noun)
26
27(defun exponentiate (pow)
28       ;;;COMPUTES %E**Z WHERE Z IS AN ARBITRARY EXPRESSION TAKING SOME OF THE WORK AWAY FROM SIMPEXPT
29  (cond ((zerop1 pow) 1)
30	((equal pow 1) '$%e)
31	(t (power '$%e pow))))
32
33(defun fixuprest (rest)
34       ;;;REST IS A PRODUCT WITHOUT THE MTIMES.FIXUPREST PUTS BACK THE MTIMES
35  (cond ((null rest) 1)
36	((cdr rest) (cons '(mtimes) rest))
37	(t (car rest))))
38
39(defun isquadraticp (e x)
40  (let ((b (sdiff e x)))
41    (cond ((zerop1 b) (list 0 0 e))
42	  ((freeof x b) (list 0 b (maxima-substitute 0 x e)))
43	  ((setq b (islinear b x))
44	   (list (div* (car b) 2) (cdr b) (maxima-substitute 0 x e))))))
45
46;;;INITIALIZES SOME GLOBAL VARIABLES THEN CALLS THE DISPATCHING FUNCTION
47
48(defmfun $laplace (fun var parm)
49  (setq fun (mratcheck fun))
50  (cond ((or *nounsflag* (member '%laplace *nounl* :test #'eq))
51         (setq fun (remlaplace fun))))
52  (cond ((and (null (atom fun)) (eq (caar fun) 'mequal))
53	 (list '(mequal)
54	       (laplace (cadr fun) parm)
55	       (laplace (caddr fun) parm)))
56	(t (laplace fun parm))))
57
58;;;LAMBDA BINDS SOME SPECIAL VARIABLES TO NIL AND DISPATCHES
59
60(defun remlaplace (e)
61  (if (atom e)
62      e
63      (cons (delete 'laplace (append (car e) nil) :count 1 :test #'eq)
64	    (mapcar #'remlaplace (cdr e)))))
65
66(defun laplace (fun parm &optional (dvar nil))
67  (let ()
68;;; Handles easy cases and calls appropriate function on others.
69    (cond ((equal fun 0) 0)
70	  ((equal fun 1)
71	   (cond ((zerop1 parm) (simplify (list '($delta) 0)))
72		 (t (power parm -1))))
73	  ((alike1 fun var) (power parm -2))
74	  ((or (atom fun) (freeof var fun))
75	   (cond ((zerop1 parm) (mul2 fun (simplify (list '($delta) 0))))
76		 (t (mul2 fun (power parm -1)))))
77	  (t
78           (let ((op (caar fun)))
79             (let ((result ; We store the result of laplace for further work.
80	       (cond ((eq op 'mplus)
81		      (laplus fun parm))
82		     ((eq op 'mtimes)
83		      (laptimes (cdr fun) parm))
84		     ((eq op 'mexpt)
85		      (lapexpt fun nil parm))
86		     ((eq op '%sin)
87		      (lapsin fun nil nil parm))
88		     ((eq op '%cos)
89		      (lapsin fun nil t parm))
90		     ((eq op '%sinh)
91		      (lapsinh fun nil nil parm))
92		     ((eq op '%cosh)
93		      (lapsinh fun nil t parm))
94		     ((eq op '%log)
95		      (laplog fun parm))
96		     ((eq op '%derivative)
97		      (lapdiff fun parm))
98		     ((eq op '%integrate)
99		      (lapint fun parm dvar))
100		     ((eq op '%sum)
101		      (list '(%sum)
102			    (laplace (cadr fun) parm)
103			    (caddr fun)
104			    (cadddr fun)
105			    (car (cddddr fun))))
106		     ((eq op '%erf)
107		      (laperf fun parm))
108		     ((and (eq op '%ilt)(eq (cadddr fun) var))
109		      (cond ((eq parm (caddr fun))(cadr fun))
110			    (t (subst parm (caddr fun)(cadr fun)))))
111                     ((eq op '$delta)
112		      (lapdelta fun nil parm))
113		     ((setq op ($get op '$laplace))
114		      (mcall op fun var parm))
115		     (t (lapdefint fun parm)))))
116              (when (isinop result '%integrate)
117                ;; Laplace has not found a result but returns a definit
118                ;; integral. This integral can contain internal integration
119                ;; variables. Replace such a result with the noun form.
120                (setq result (list '(%laplace) fun var parm)))
121              ;; Check if we have a result, when not call $specint.
122              (check-call-to-$specint result fun parm)))))))
123
124;;; Check if laplace has found a result, when not try $specint.
125
126(defun check-call-to-$specint (result fun parm)
127  (cond
128    ((or (isinop result '%laplace)
129         (isinop result '%limit)   ; Try $specint for incomplete results
130         (isinop result '%at))     ; which contain %limit or %at too.
131     ;; laplace returns a noun form or a result which contains %limit or %at.
132     ;; We pass the function to $specint to look for more results.
133     (let (res)
134       ;; laplace assumes the parameter s to be positive and does a
135       ;; declaration before an integration is done. Therefore we declare
136       ;; the parameter of the Laplace transform to be positive before
137       ;; we call $specint too.
138       (with-new-context (context)
139         (progn
140           (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
141           (setq res ($specint (mul fun (power '$%e (mul -1 var parm))) var))))
142       (if (or (isinop res '%specint)  ; Both symobls are possible, that is
143               (isinop res '$specint)) ; not consistent! Check it! 02/2009
144           ;; $specint has not found a result.
145           result
146           ;; $specint has found a result
147           res)))
148       (t result)))
149
150(defun laplus (fun parm)
151  (simplus (cons '(mplus) (mapcar #'(lambda (e) (laplace e parm)) (cdr fun))) 1 t))
152
153(defun laptimes (fun parm)
154       ;;;EXPECTS A LIST (PERHAPS EMPTY) OF FUNCTIONS MULTIPLIED TOGETHER WITHOUT THE MTIMES
155       ;;;SEES IF IT CAN APPLY THE FIRST AS A TRANSFORMATION ON THE REST OF THE FUNCTIONS
156  (cond ((null fun) (list '(mexpt) parm -1.))
157	((null (cdr fun)) (laplace (car fun) parm))
158	((freeof var (car fun))
159	 (simptimes (list '(mtimes)
160			  (car fun)
161			  (laptimes (cdr fun) parm))
162		    1
163		    t))
164	((eq (car fun) var)
165	 (simptimes (list '(mtimes) -1 (sdiff (laptimes (cdr fun) parm) parm))
166		    1
167		    t))
168	(t
169	 (let ((op (caaar fun)))
170	   (cond ((eq op 'mexpt)
171		  (lapexpt (car fun) (cdr fun) parm))
172		 ((eq op 'mplus)
173		  (laplus ($multthru (fixuprest (cdr fun)) (car fun)) parm))
174		 ((eq op '%sin)
175		  (lapsin (car fun) (cdr fun) nil parm))
176		 ((eq op '%cos)
177		  (lapsin (car fun) (cdr fun) t parm))
178		 ((eq op '%sinh)
179		  (lapsinh (car fun) (cdr fun) nil parm))
180		 ((eq op '%cosh)
181		  (lapsinh (car fun) (cdr fun) t parm))
182		 ((eq op '$delta)
183		  (lapdelta (car fun) (cdr fun) parm))
184		 (t
185		  (lapshift (car fun) (cdr fun) parm)))))))
186
187(defun lapexpt (fun rest parm)
188       ;;;HANDLES %E**(A*T+B)*REST(T), %E**(A*T**2+B*T+C),
189       ;;; 1/SQRT(A*T+B), OR T**K*REST(T)
190  (prog (ab base-of-fun power result)
191     (setq base-of-fun (cadr fun) power (caddr fun))
192     (cond
193       ((and
194	 (freeof var base-of-fun)
195	 (setq
196	  ab
197	  (isquadraticp
198	   (cond ((eq base-of-fun '$%e) power)
199		 (t (simptimes (list '(mtimes)
200				     power
201				     (list '(%log)
202					   base-of-fun))
203			       1
204			       nil)))
205	   var)))
206	(cond ((equal (car ab) 0) (go %e-case-lin))
207	      ((null rest) (go %e-case-quad))
208	      (t (go noluck))))
209       ((and (eq base-of-fun var) (freeof var power))
210	(go var-case))
211       ((and (alike1 '((rat) -1 2) power) (null rest)
212	     (setq ab (islinear base-of-fun var)))
213	(setq result (div* (cdr ab) (car ab)))
214	(return (simptimes
215		 (list '(mtimes)
216		       (list '(mexpt)
217			     (div* '$%pi
218				   (list '(mtimes)
219					 (car ab)
220					 parm))
221			     '((rat) 1 2))
222		       (exponentiate (list '(mtimes) result parm))
223		       (list '(mplus)
224			     1
225			     (list '(mtimes)
226				   -1
227				   (list '(%erf)
228					 (list '(mexpt)
229					       (list '(mtimes)
230						     result
231						     parm)
232					       '((rat) 1 2)))
233				   ))) 1 nil)))
234       (t (go noluck)))
235     %e-case-lin
236     (setq result
237      (cond
238	(rest (sratsimp ($at (laptimes rest parm)
239			     (list '(mequal)
240				   parm
241				   (list '(mplus)
242					 parm
243					 (afixsign (cadr ab)
244						   nil))))))
245	(t (list '(mexpt)
246		 (list '(mplus)
247		       parm
248		       (afixsign (cadr ab) nil))
249		 -1))))
250     (return (simptimes (list '(mtimes)
251			      (exponentiate (caddr ab))
252			      result)
253			1
254			nil))
255     %e-case-quad
256     (setq result (afixsign (car ab) nil))
257     (setq
258      result
259      (list
260       '(mtimes)
261       (div* (list '(mexpt)
262		   (div* '$%pi result)
263		   '((rat) 1 2))
264	     2)
265       (exponentiate (div* (list '(mexpt) parm 2)
266			   (list '(mtimes) 4 result)))
267       (list '(mplus)
268	     1
269	     (list '(mtimes)
270		   -1
271		   (list '(%erf)
272			 (div* parm
273			       (list '(mtimes)
274				     2
275				     (list '(mexpt)
276					   result
277					   '((rat) 1 2)))))
278		   ))))
279     (and (null (equal (cadr ab) 0))
280	  (setq result
281		(maxima-substitute (list '(mplus)
282					 parm
283					 (list '(mtimes)
284					       -1
285					       (cadr ab)))
286				   parm
287				   result)))
288     (return (simptimes  (list '(mtimes)
289			       (exponentiate (caddr ab))
290			       result) 1 nil))
291     var-case
292     (cond ((or (null rest) (freeof var (fixuprest rest)))
293	    (go var-easy-case)))
294     (cond ((posint power)
295	    (return (afixsign (apply '$diff
296				     (list (laptimes rest parm)
297					   parm
298					   power))
299			      (even power))))
300	   ((negint power)
301	    (return (mydefint (hackit power rest parm)
302			      (createname parm (- power))
303			      parm parm)))
304	   (t (go noluck)))
305     var-easy-case
306     (setq power
307	   (simplus (list '(mplus) 1 power) 1 t))
308     (or (eq (asksign power) '$positive) (go noluck))
309     (setq result (list (list '(%gamma) power)
310			(list '(mexpt)
311			      parm
312			      (afixsign power nil))))
313     (and rest (setq result (nconc result rest)))
314     (return (simptimes (cons '(mtimes) result) 1 nil))
315     noluck
316     (return
317       (cond
318	 ((and (posint power)
319	       (member (caar base-of-fun)
320		     '(mplus %sin %cos %sinh %cosh) :test #'eq))
321	  (laptimes (cons base-of-fun
322			  (cons (cond ((= power 2) base-of-fun)
323				      (t (list '(mexpt)
324					       base-of-fun
325					       (1- power))))
326				rest)) parm))
327	 (t (lapshift fun rest parm))))))
328
329;;;INTEGRAL FROM A TO INFINITY OF F(X)
330(defun mydefint (f x a parm)
331  (let ((tryint (and (not ($unknown f))
332                     ;; $defint should not throw a Maxima error,
333                     ;; therefore we set the flags errcatch and $errormsg.
334                     ;; errset catches the error and returns nil
335                     (with-new-context (context)
336                       (progn
337                         (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
338                         (meval `(($assume) ,@(list (list '(mgreaterp) x 0))))
339                         (meval `(($assume) ,@(list (list '(mgreaterp) a 0))))
340                         (let ((errcatch t) ($errormsg nil))
341                           (errset ($defint f x a '$inf))))))))
342    (if tryint
343	(car tryint)
344	(list '(%integrate) f x a '$inf))))
345
346 ;;;CREATES UNIQUE NAMES FOR VARIABLE OF INTEGRATION
347(defun createname (head tail)
348  (intern (format nil "~S~S" head tail)))
349
350;;;REDUCES LAPLACE(F(T)/T**N,T,S) CASE TO LAPLACE(F(T)/T**(N-1),T,S) CASE
351(defun hackit (exponent rest parm)
352  (cond ((equal exponent -1)
353	 (let ((parm (createname parm 1)))
354	   (laptimes rest parm)))
355	(t (mydefint (hackit (1+ exponent) rest parm)
356		     (createname parm (- -1 exponent))
357		     (createname parm (- exponent)) parm))))
358
359(defun afixsign (funct signswitch)
360       ;;;MULTIPLIES FUNCT BY -1 IF SIGNSWITCH IS NIL
361  (cond (signswitch funct)
362	(t (simptimes (list '(mtimes) -1 funct) 1 t))))
363
364(defun lapshift (fun rest parm)
365  (cond ((atom fun) (merror "LAPSHIFT: expected a cons, not ~M" fun))
366	((or (member 'laplace (car fun) :test #'eq) (null rest))
367	 (lapdefint (cond (rest (simptimes (cons '(mtimes)
368						 (cons fun rest)) 1 t))
369			  (t fun)) parm))
370	(t (laptimes (append rest
371			     (ncons (cons (append (car fun)
372						  '(laplace))
373					  (cdr fun)))) parm))))
374
375;;;COMPUTES %E**(W*B*%I)*F(S-W*A*%I) WHERE W=-1 IF SIGN IS T ELSE W=1
376(defun mostpart (f parm sign a b)
377  (let ((substinfun ($at f
378			 (list '(mequal)
379			       parm
380			       (list '(mplus) parm (afixsign (list '(mtimes) a '$%i) sign))))))
381    (if (zerop1 b)
382	substinfun
383	(list '(mtimes)
384	      (exponentiate (afixsign (list '(mtimes) b '$%i) (null sign)))
385	      substinfun))))
386
387 ;;;IF WHICHSIGN IS NIL THEN SIN TRANSFORM ELSE COS TRANSFORM
388(defun compose (fun parm whichsign a b)
389  (let ((result (list '((rat) 1 2)
390		      (list '(mplus)
391			    (mostpart fun parm t a b)
392			    (afixsign (mostpart fun parm nil a b)
393				      whichsign)))))
394    (sratsimp (simptimes (cons '(mtimes)
395			       (if whichsign
396				   result
397				   (cons '$%i result)))
398			 1 nil))))
399
400 ;;;FUN IS OF THE FORM SIN(A*T+B)*REST(T) OR COS
401(defun lapsin (fun rest trigswitch parm)
402  (let ((ab (islinear (cadr fun) var)))
403    (cond (ab
404	   (cond (rest
405		  (compose (laptimes rest parm)
406			   parm
407			   trigswitch
408			   (car ab)
409			   (cdr ab)))
410		 (t
411		  (simptimes
412		   (list '(mtimes)
413		    (cond ((zerop1 (cdr ab))
414			   (if trigswitch parm (car ab)))
415			  (t
416			   (cond (trigswitch
417				  (list '(mplus)
418					(list '(mtimes)
419					      parm
420					      (list '(%cos) (cdr ab)))
421					(list '(mtimes)
422					      -1
423					      (car ab)
424					      (list '(%sin) (cdr ab)))))
425				 (t
426				  (list '(mplus)
427					(list '(mtimes)
428					      parm
429					      (list '(%sin) (cdr ab)))
430					(list '(mtimes)
431					      (car ab)
432					      (list '(%cos) (cdr ab))))))))
433		    (list '(mexpt)
434			  (list '(mplus)
435				(list '(mexpt) parm 2)
436				(list '(mexpt) (car ab) 2))
437			  -1))
438		   1 nil))))
439	  (t
440	   (lapshift fun rest parm)))))
441
442 ;;;FUN IS OF THE FORM SINH(A*T+B)*REST(T) OR IS COSH
443(defun lapsinh (fun rest switch parm)
444  (cond ((islinear (cadr fun) var)
445	 (sratsimp
446	  (laplus
447	   (simplus
448	    (list '(mplus)
449		  (nconc (list '(mtimes)
450			       (list '(mexpt)
451				     '$%e
452				     (cadr fun))
453			       '((rat) 1 2))
454			 rest)
455		  (afixsign (nconc (list '(mtimes)
456					 (list '(mexpt)
457					       '$%e
458					       (afixsign (cadr fun)
459							 nil))
460					 '((rat) 1 2))
461				   rest)
462			    switch))
463	    1
464	    nil) parm)))
465	(t (lapshift fun rest parm))))
466
467 ;;;FUN IS OF THE FORM LOG(A*T)
468(defun laplog (fun parm)
469  (let ((ab (islinear (cadr fun) var)))
470    (cond ((and ab (zerop1 (cdr ab)))
471	   (simptimes (list '(mtimes)
472			    (list '(mplus)
473				  (subfunmake '$psi '(0) (ncons 1))
474				  (list '(%log) (car ab))
475				  (list '(mtimes) -1 (list '(%log) parm)))
476			    (list '(mexpt) parm -1))
477		      1 nil))
478	  (t
479	   (lapdefint fun parm)))))
480
481(defun raiseup (fbase exponent)
482  (if (equal exponent 1)
483      fbase
484      (list '(mexpt) fbase exponent)))
485
486;;TAKES TRANSFORM OF DELTA(A*T+B)*F(T)
487(defun lapdelta (fun rest parm)
488  (let ((ab (islinear (cadr fun) var))
489	(sign nil)
490	(recipa nil))
491    (cond (ab
492	   (setq recipa (power (car ab) -1) ab (div (cdr ab) (car ab)))
493	   (setq sign (asksign ab) recipa (simplifya (list '(mabs) recipa) nil))
494	   (simplifya (cond ((eq sign '$positive)
495			     0)
496			    ((eq sign '$zero)
497			     (list '(mtimes)
498				   (maxima-substitute 0 var (fixuprest rest))
499				   recipa))
500			    (t
501			     (list '(mtimes)
502				   (maxima-substitute (neg ab) var (fixuprest rest))
503				   (list '(mexpt) '$%e (cons '(mtimes) (cons parm (ncons ab))))
504				   recipa)))
505		      nil))
506	  (t
507	   (lapshift fun rest parm)))))
508
509(defun laperf (fun parm)
510  (let ((ab (islinear (cadr fun) var)))
511    (cond ((and ab (equal (cdr ab) 0))
512	   (simptimes (list '(mtimes)
513			    (div* (exponentiate (div* (list '(mexpt) parm 2)
514						      (list '(mtimes)
515							    4
516							    (list '(mexpt) (car ab) 2))))
517				  parm)
518			    (list '(mplus)
519				  1
520				  (list '(mtimes)
521					-1
522					(list '(%erf) (div* parm (list '(mtimes) 2 (car ab)))))))
523		      1
524		      nil))
525	  (t
526	   (lapdefint fun parm)))))
527
528(defun lapdefint (fun parm)
529  (prog (tryint mult)
530     (and ($unknown fun)(go skip))
531     (setq mult (simptimes (list '(mtimes) (exponentiate
532					    (list '(mtimes) -1 var parm)) fun) 1 nil))
533     (with-new-context (context)
534       (progn
535         (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
536         (setq tryint
537               ;; $defint should not throw a Maxima error.
538               ;; therefore we set the flags errcatch and errormsg.
539               ;; errset catches an error and returns nil.
540               (let ((errcatch t) ($errormsg nil))
541                 (errset ($defint mult var 0 '$inf))))))
542     (and tryint (not (eq (and (listp (car tryint))
543			       (caaar tryint))
544			  '%integrate))
545	  (return (car tryint)))
546     skip (return (list '(%laplace) fun var parm))))
547
548
549(defun lapdiff (fun parm)
550;;;FUN IS OF THE FORM DIFF(F(T),T,N) WHERE N IS A POSITIVE INTEGER
551  (prog (difflist degree frontend resultlist newdlist order arg2)
552     (setq newdlist (setq difflist (copy-tree (cddr fun))))
553     (setq arg2 (list '(mequal) var 0))
554     a    (cond ((null difflist)
555		 (return (cons '(%derivative)
556			       (cons (list '(%laplace)
557					   (cadr fun)
558					   var
559					   parm)
560				     newdlist))))
561		((eq (car difflist) var)
562		 (setq degree (cadr difflist)
563		       difflist (cddr difflist))
564		 (go out)))
565     (setq difflist (cdr (setq frontend (cdr difflist))))
566     (go a)
567     out  (cond ((null (posint degree))
568		 (return (list '(%laplace) fun var parm))))
569     (cond (frontend (rplacd frontend difflist))
570	   (t (setq newdlist difflist)))
571     (cond (newdlist (setq fun (cons '(%derivative)
572				     (cons (cadr fun)
573					   newdlist))))
574	   (t (setq fun (cadr fun))))
575     (setq order 0)
576     loop (decf degree)
577     (setq resultlist
578	   (cons (list '(mtimes)
579		       (raiseup parm degree)
580		       ($at ($diff fun var order) arg2))
581		 resultlist))
582     (incf order)
583     (and (> degree 0) (go loop))
584     (setq resultlist (cond ((cdr resultlist)
585			     (cons '(mplus)
586				   resultlist))
587			    (t (car resultlist))))
588     (return (simplus (list '(mplus)
589			    (list '(mtimes)
590				  (raiseup parm order)
591				  (laplace fun parm))
592			    (list '(mtimes)
593				  -1
594				  resultlist))
595		      1 nil))))
596
597 ;;;FUN IS OF THE FORM INTEGRATE(F(X)*G(T)*H(T-X),X,0,T)
598(defun lapint (fun parm dvar)
599  (prog (newfun parm-list f var-list var-parm-list)
600     (and dvar (go convolution))
601     (setq dvar (cadr (setq newfun (cdr fun))))
602     (and (cddr newfun)
603	  (zerop1 (caddr newfun))
604	  (eq (cadddr newfun) var)
605	  (go convolutiontest))
606     notcon
607     (setq newfun (cdr fun))
608     (cond ((cddr newfun)
609	    (cond ((and (freeof var (caddr newfun))
610			(freeof var (cadddr newfun)))
611		   (return (list '(%integrate)
612				 (laplace (car newfun) parm dvar)
613				 dvar
614				 (caddr newfun)
615				 (cadddr newfun))))
616		  (t (go giveup))))
617	   (t (return (list '(%integrate)
618			    (laplace (car newfun) parm dvar)
619			    dvar))))
620     giveup
621     (return (list '(%laplace) fun var parm))
622     convolutiontest
623     (setq newfun ($factor (car newfun)))
624     (cond ((eq (caar newfun) 'mtimes)
625	    (setq f (cadr newfun) newfun (cddr newfun)))
626	   (t (setq f newfun newfun nil)))
627     gothrulist
628     (cond ((freeof dvar f)
629	    (setq parm-list (cons f parm-list)))
630	   ((freeof var f) (setq var-list (cons f var-list)))
631	   ((freeof dvar
632		    (sratsimp (maxima-substitute (list '(mplus)
633						       var
634						       dvar)
635						 var
636						 f)))
637	    (setq var-parm-list (cons f var-parm-list)))
638	   (t (go notcon)))
639     (cond (newfun (setq f (car newfun) newfun (cdr newfun))
640		   (go gothrulist)))
641     (and
642      parm-list
643      (return
644	(laplace
645	 (cons
646	  '(mtimes)
647	  (nconc parm-list
648		 (ncons (list '(%integrate)
649			      (cons '(mtimes)
650				    (append var-list
651					    var-parm-list))
652			      dvar
653			      0
654			      var)))) parm dvar)))
655     convolution
656     (return
657       (simptimes
658	(list
659	 '(mtimes)
660	 (laplace ($expand (maxima-substitute var
661					      dvar
662					      (fixuprest var-list))) parm dvar)
663	 (laplace
664	  ($expand (maxima-substitute 0 dvar (fixuprest var-parm-list))) parm dvar))
665	1
666	t))))
667
668(declare-top (special varlist ratform ils ilt))
669
670(defmfun $ilt (exp ils ilt)
671 ;;;EXP IS F(S)/G(S) WHERE F AND G ARE POLYNOMIALS IN S AND DEGR(F) < DEGR(G)
672  (let (varlist ($savefactors t) checkfactors $ratfac $keepfloat)
673		;;; MAKES ILS THE MAIN VARIABLE
674    (setq varlist (list ils))
675    (newvar exp)
676    (orderpointer varlist)
677    (setq var (caadr (ratrep* ils)))
678    (cond ((and (null (atom exp))
679		(eq (caar exp) 'mequal))
680	   (list '(mequal)
681		 ($ilt (cadr exp) ils ilt)
682		 ($ilt (caddr exp) ils ilt)))
683	  ((zerop1 exp) 0)
684	  ((freeof ils exp)
685	   (list '(%ilt) exp ils ilt))
686	  (t (ilt0 exp)))))
687
688(defun maxima-rationalp (le v)
689  (cond ((null le))
690	((and (null (atom (car le))) (null (freeof v (car le))))
691	 nil)
692	(t (maxima-rationalp (cdr le) v))))
693
694 ;;;THIS FUNCTION DOES THE PARTIAL FRACTION DECOMPOSITION
695(defun ilt0 (exp)
696  (prog (wholepart frpart num denom y content real factor
697	 apart bpart parnumer ratarg ratform)
698     (and (mplusp exp)
699	  (return (simplus  (cons '(mplus)
700				  (mapcar #'(lambda (f) ($ilt f ils ilt)) (cdr exp))) 1 t)))
701     (and (null (atom exp))
702	  (eq (caar exp) '%laplace)
703	  (eq (cadddr exp) ils)
704	  (return (cond ((eq (caddr exp) ilt) (cadr exp))
705			(t (subst ilt
706				  (caddr exp)
707				  (cadr exp))))))
708     (setq ratarg (ratrep* exp))
709     (or (maxima-rationalp varlist ils)
710	 (return (list '(%ilt) exp ils ilt)))
711     (setq ratform (car ratarg))
712     (setq denom (ratdenominator (cdr ratarg)))
713     (setq frpart (pdivide (ratnumerator (cdr ratarg)) denom))
714     (setq wholepart (car frpart))
715     (setq frpart (ratqu (cadr frpart) denom))
716     (cond ((not (zerop1 (car wholepart)))
717	    (return (list '(%ilt) exp ils ilt)))
718	   ((zerop1 (car frpart)) (return 0)))
719     (setq num (car frpart) denom (cdr frpart))
720     (setq y (oldcontent denom))
721     (setq content (car y))
722     (setq real (cadr y))
723     (setq factor (pfactor real))
724     loop (cond ((null (cddr factor))
725		 (setq apart real
726		       bpart 1
727		       y '((0 . 1) 1 . 1))
728		 (go skip)))
729     (setq apart (pexpt (car factor) (cadr factor)))
730     (setq bpart (car (ratqu real apart)))
731     (setq y (bprog apart bpart))
732     skip (setq frpart
733		(cdr (ratdivide (ratti (ratnumerator num)
734				       (cdr y)
735				       t)
736				(ratti (ratdenominator num)
737				       (ratti content apart t)
738				       t))))
739     (setq
740      parnumer
741      (cons (ilt1 (ratqu (ratnumerator frpart)
742			 (ratti (ratdenominator frpart)
743				(ratti (ratdenominator num)
744				       content
745				       t)
746				t))
747		  (car factor)
748		  (cadr factor))
749	    parnumer))
750     (setq factor (cddr factor))
751     (cond ((null factor)
752	    (return (simplus (cons '(mplus) parnumer) 1 t))))
753     (setq num (cdr (ratdivide (ratti num (car y) t)
754			       (ratti content bpart t))))
755     (setq real bpart)
756     (go loop)))
757
758(declare-top (special q z))
759
760(defun ilt1 (p q k)
761  (let (z)
762    (cond ((onep1 k) (ilt3 p ))
763	  (t (setq z (bprog q (pderivative q var)))
764	     (ilt2 p k)))))
765
766
767 ;;;INVERTS P(S)/Q(S)**K WHERE Q(S)  IS IRREDUCIBLE
768 ;;;DOESN'T CALL ILT3 IF Q(S) IS LINEAR
769(defun ilt2 (p k)
770  (prog (y a b)
771     (and (onep1 k) (return (ilt3 p)))
772     (decf k)
773     (setq a (ratti p (car z) t))
774     (setq b (ratti p (cdr z) t))
775     (setq y (pexpt q k))
776     (cond
777       ((or (null (equal (pdegree q var) 1))
778	    (> (pdegree (car p) var) 0))
779	(return
780	  (simplus
781	   (list
782	    '(mplus)
783	    (ilt2
784	     (cdr (ratdivide (ratplus a (ratqu (ratderivative b var) k)) y))
785	     k)
786	    ($multthru (simptimes (list '(mtimes)
787					ilt
788					(power k -1)
789					(ilt2 (cdr (ratdivide b y)) k))
790				  1
791				  t)))
792	   1
793	   t))))
794     (setq a (disrep (polcoef q 1))
795	   b (disrep (polcoef q 0)))
796     (return
797       (simptimes (list '(mtimes)
798			(disrep p)
799			(raiseup ilt k)
800			(simpexpt (list '(mexpt)
801					'$%e
802					(list '(mtimes)
803					      -1
804					      ilt
805					      b
806					      (list '(mexpt) a -1)))
807				  1
808				  nil)
809			(list '(mexpt)
810			      a
811			      (- -1 k))
812			(list '(mexpt)
813			      (factorial k)
814			      -1))
815		  1
816		  nil))))
817
818(declare-top(notype k))
819
820;;(DEFUN COEF MACRO (POL) (SUBST (CADR POL) (QUOTE DEG)
821;;  '(DISREP (RATQU (POLCOEF (CAR P) DEG) (CDR P)))))
822
823(defmacro coef (pol)
824  `(disrep (ratqu (polcoef (car p) ,pol) (cdr p))))
825
826(defun lapsum (&rest args)
827  (cons '(mplus) args))
828
829(defun lapprod (&rest args)
830  (cons '(mtimes) args))
831
832(defun expo (&rest args)
833  (cons '(mexpt) args))
834
835;;;INVERTS P(S)/Q(S) WHERE Q(S) IS IRREDUCIBLE
836(defun ilt3 (p)
837  (prog (discrim sign a c d e b1 b0 r term1 term2 degr)
838     (setq e (disrep (polcoef q 0))
839	   d (disrep (polcoef q 1))
840	   degr (pdegree q var))
841     (and (equal degr 1)
842	  (return
843	    (simptimes (lapprod
844			(disrep p)
845			(expo d -1)
846			(expo '$%e (lapprod -1 ilt e (expo d -1))))
847		       1
848		       nil)))
849     (setq c (disrep (polcoef q 2)))
850     (and (equal degr 2) (go quadratic))
851     (and (equal degr 3) (zerop1 c) (zerop1 d)
852	  (go cubic))
853     (return (list '(%ilt) (div* (disrep p)(disrep q)) ils ilt))
854     cubic (setq  a (disrep (polcoef q 3))
855		  r (simpnrt (div* e a) 3))
856     (setq d (div* (disrep p)(lapprod a (lapsum
857					 (expo ils 3)(expo '%r 3)))))
858     (return (ilt0 (maxima-substitute r '%r ($partfrac d ils))))
859     quadratic (setq b0 (coef 0) b1 (coef 1))
860
861     (setq discrim
862	   (simplus (lapsum
863		     (lapprod 4 e c)
864		     (lapprod -1 d d))
865		    1
866		    nil))
867     (setq sign (cond ((free discrim '$%i) (asksign discrim)) (t '$positive))
868	   term1 '(%cos)
869	   term2 '(%sin))
870     (setq degr (expo '$%e (lapprod ilt d (power c -1) '((rat) -1 2))))
871     (cond ((eq sign '$zero)
872	    (return (simptimes (lapprod degr (lapsum (div* b1 c)
873						     (lapprod
874						      (div* (lapsum (lapprod 2 b0 c) (lapprod -1 b1 d))
875							    (lapprod 2 c c)) ilt))) 1 nil))
876	    )		   ((eq sign '$negative)
877	    (setq term1 '(%cosh)
878		  term2 '(%sinh)
879		  discrim (simptimes (lapprod -1 discrim) 1 t))))
880     (setq discrim (simpnrt discrim 2))
881     (setq sign
882      (simptimes
883       (lapprod
884	(lapsum
885	 (lapprod 2 b0 c)
886	 (lapprod -1 b1 d))
887	(expo discrim -1))
888       1
889       nil))
890     (setq c (power c -1))
891     (setq discrim (simptimes (lapprod
892			       discrim
893			       ilt
894			       '((rat) 1 2)
895			       c)
896			      1
897			      t))
898     (return
899       (simptimes
900	(lapprod c degr
901	 (lapsum
902	  (lapprod b1 (list term1 discrim))
903	  (lapprod sign (list term2 discrim))))
904	1
905	nil))))
906
907(declare-top (unspecial ils ilt *nounl* q ratform var
908			varlist z))
909