;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package "MAXIMA")
(macsyma-module pade)

(declare-top(special tlist tpf))

(load-macsyma-macros rzmac mhayat ratmac)

(defmacro red (p) `(cdr ,p))

(defmacro p0? (p) `(null ,p))

(defmacro e0? (expon) `(zerop (car ,expon)))

(defmacro epos? (expon) `(signp g (car ,expon)))

(defmacro eneg? (expon) `(signp l (car ,expon)))

(defmacro num (r) `(car ,r))

(defmacro denom (r) `(cdr ,r))

(defmacro eneg (exp) `(cons (minus (car ,exp)) (cdr ,exp)))

(defmacro pade-lexp (poly) `(cond ((p0? ,poly) (ezero))
				  (t (le ,poly))))

(defun eshift (poly expon)
       (mapcar #'(lambda (tm) (term (e+ (e tm) expon)
				    (c tm)))
	poly))

(defmacro psmake (p)
	  `(make-ps (gvar-o tpf) (poly-data tpf) (reverse ,p)))

(defmfun $pade (taylor-form n-bound d-bound)
       (cond ((not (and ($ratp taylor-form) (memq 'trunc (car taylor-form))))
	      (merror "First argument to Pade must be a taylor series")))
       (let (((nil nil varlist genvar tlist) (car taylor-form)))
	    (cons '(mlist)
		  (pade (cdr taylor-form)
			(cons (car (tay-order n-bound))
			      (car (tay-order d-bound)))
			(orig-trunc (car (mrat-tlist taylor-form)))))))
	     
(defun pade (tpf bounds trunc)
  (let ((ord (le (terms tpf))) (poly (terms tpf)))
       (setq trunc (e+ trunc (egcd trunc (psexpon-gcd poly))))
       (unless (e0? ord)
	       (setq poly (eshift poly (eneg ord)))
	       (setq trunc (e- trunc ord))
	       (cond ((epos? ord) (rplaca bounds
					  (e- (car bounds) ord)))
		     (t (rplacd bounds (e+ (cdr bounds) ord)))))
       (let ((r1 (cons (list (term trunc (rcone))) nil))
	     (r2 (pade-monize (reverse poly) (list (term (ezero) (rcone))))))
	    (do ((ans (cond ((or (e> (le (last poly)) (num bounds))
				 (eneg? (denom bounds))) nil)
			    (t (list r2)))
		      (cond ((or (e> (le (num r2)) (num bounds))
				 (e> (le (denom r2)) (denom bounds)))
			     ans)
			    ((epos? (le (last (num r2))))
			     ans)			;not a true approximate
			    (t (cons r2 ans)))))
		((or (e0? (pade-lexp (num r2)))
		     (e>= (pade-lexp (denom r2)) (denom bounds)))
		 (mapcar #'(lambda (r) ($ratsimp
					(m// (srdis (psmake (eshift (num r)
								    ord)))
					     (srdis (psmake (denom r))))))
			 ans))
	      (setq r1 (pade1 r1 r2))
	      (exch r1 r2)))))

(defun pade1 (r1 r2)
       (do ((quoterm) (lcoef (lc (num r2))) (expon (le (num r2)))
	    (p2 (red (num r2))) (q2 (denom r2)) (p1 (num r1)) (q1 (denom r1)))
	   ((or (p0? p1) (e> expon (le p1)))
		 (pade-monize p1 q1))
	   (setq quoterm (term (e- (le p1) expon)
			    (rcminus (rcquo (lc p1) lcoef))))
	   (setq p1 (upoly+ (red p1)
			 (term* quoterm p2)))
	   (setq q1 (upoly+ q1 (term* quoterm q2)))))

(defun term* (term p)
       (let ((p1 nil) (coef (c term)) (expon (e term)) coef1)
	     (mapc #'(lambda (term1)
			     (setq coef1 (rctimes coef (c term1)))
			     (if (not (rczerop coef1))
				 (push (term (e+ expon (e term1))
					     coef1)
				       p1)))
		   p)
	     (nreverse p1)))

(defun pade-monize (num den)				;monicize numerator
       (let ((mult (term (ezero) (rcinv (lc num)))))
	    (cons (term* mult num) (term* mult den))))
	     
(defun upoly+ (p1 p2)
	(do ((p) (c))
	    (nil)
	    (cond ((p0? p1) (return (nreconc p p2)))
		  ((p0? p2) (return (nreconc p p1)))
		  ((e> (le p2) (le p1)) (push (pop p2) p))
		  ((e> (le p1) (le p2)) (push (pop p1) p))
		  (t
		   (setq c (rcplus (lc p1) (lc p2)))
		   (if (not (rczerop c))
		       (push (term  (le p1) c) p))
		   (setq p1 (red p1) p2 (red p2))))))

#-NIL
(declare-top(unspecial tlist tpf))