src/pade.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module pade)
  14
  15 (declare-top (special tlist))
  16
  17 (load-macsyma-macros rzmac mhayat ratmac)
  18
  19 (defmacro red (p) `(cdr ,p))
  20
  21 (defmacro p0? (p) `(null ,p))
  22
  23 (defmacro e0? (expon) `(zerop (car ,expon)))
  24
  25 (defmacro epos? (expon) `(signp g (car ,expon)))
  26
  27 (defmacro eneg? (expon) `(signp l (car ,expon)))
  28
  29 (defmacro num (r) `(car ,r))
  30
  31 (defmacro denom (r) `(cdr ,r))
  32
  33 (defmacro eneg (exp) `(cons (- (car ,exp)) (cdr ,exp)))
  34
  35 (defmacro pade-lexp (poly) `(cond ((p0? ,poly) (ezero))
  36                              (t (le ,poly))))
  37
  38 (defun eshift (poly expon)
  39   (mapcar #'(lambda (tm) (term (e+ (e tm) expon)
  40                                (c tm)))
  41           poly))
  42
  43 (defmacro psmake (p tpf)
  44   `(make-ps (gvar-o ,tpf) (poly-data ,tpf) (reverse ,p)))
  45
  46 (defmfun $pade (taylor-form n-bound d-bound)
  47   (cond ((not (and ($ratp taylor-form) (member 'trunc (car taylor-form) :test #'eq)))
  48          (merror (intl:gettext "pade: first argument must be a Taylor series; found: ~M") taylor-form)))
  49   (destructuring-let (((nil nil varlist genvar tlist) (car taylor-form)))
  50     (if (psp (cdr taylor-form))
  51         (cons '(mlist)
  52               (pade (cdr taylor-form)
  53                     (cons (car (tay-order n-bound))
  54                           (car (tay-order d-bound)))
  55                     (orig-trunc (car (mrat-tlist taylor-form)))))
  56       (list '(mlist) (srdis (cdr taylor-form))))))
  57
  58 (defun pade (tpf bounds trunc)
  59   (let ((ord (le (terms tpf))) (poly (terms tpf)))
  60     (setq trunc (e+ trunc (egcd trunc (psexpon-gcd poly))))
  61     (unless (e0? ord)
  62       (setq poly (eshift poly (eneg ord)))
  63       (setq trunc (e- trunc ord))
  64       (cond ((epos? ord) (rplaca bounds
  65                                  (e- (car bounds) ord)))
  66             (t (rplacd bounds (e+ (cdr bounds) ord)))))
  67     (let ((r1 (cons (list (term trunc (rcone))) nil))
  68           (r2 (pade-monize (reverse poly) (list (term (ezero) (rcone))))))
  69       (do ((ans (cond ((or (e> (le (last poly)) (num bounds))
  70                            (eneg? (denom bounds))) nil)
  71                       (t (list r2)))
  72                 (cond ((or (e> (le (num r2)) (num bounds))
  73                            (e> (le (denom r2)) (denom bounds)))
  74                        ans)
  75                       ((epos? (le (last (num r2))))
  76                        ans)             ;not a true approximate
  77                       (t (cons r2 ans)))))
  78           ((or (e0? (pade-lexp (num r2)))
  79                (e>= (pade-lexp (denom r2)) (denom bounds)))
  80            (mapcar #'(lambda (r) ($ratsimp
  81                                   (m// (srdis (psmake (eshift (num r) ord) tpf))
  82                                        (srdis (psmake (denom r) tpf)))))
  83                    ans))
  84         (setq r1 (pade1 r1 r2))
  85         (rotatef r1 r2)))))
  86
  87 (defun pade1 (r1 r2)
  88   (do ((quoterm) (lcoef (lc (num r2))) (expon (le (num r2)))
  89        (p2 (red (num r2))) (q2 (denom r2)) (p1 (num r1)) (q1 (denom r1)))
  90       ((or (p0? p1) (e> expon (le p1)))
  91        (pade-monize p1 q1))
  92     (setq quoterm (term (e- (le p1) expon)
  93                         (rcminus (rcquo (lc p1) lcoef))))
  94     (setq p1 (upoly+ (red p1)
  95                      (term* quoterm p2)))
  96     (setq q1 (upoly+ q1 (term* quoterm q2)))))
  97
  98 (defun term* (term p)
  99   (let ((p1 nil) (coef (c term)) (expon (e term)) coef1)
 100     (mapc #'(lambda (term1)
 101               (setq coef1 (rctimes coef (c term1)))
 102               (if (not (rczerop coef1))
 103                   (push (term (e+ expon (e term1))
 104                               coef1)
 105                         p1)))
 106           p)
 107     (nreverse p1)))
 108
 109 (defun pade-monize (num den)            ;monicize numerator
 110   (let ((mult (term (ezero) (rcinv (lc num)))))
 111     (cons (term* mult num) (term* mult den))))
 112
 113 (defun upoly+ (p1 p2)
 114   (do ((p) (c))
 115       (nil)
 116     (cond ((p0? p1) (return (nreconc p p2)))
 117           ((p0? p2) (return (nreconc p p1)))
 118           ((e> (le p2) (le p1)) (push (pop p2) p))
 119           ((e> (le p1) (le p2)) (push (pop p1) p))
 120           (t
 121            (setq c (rcplus (lc p1) (lc p2)))
 122            (if (not (rczerop c))
 123                (push (term  (le p1) c) p))
 124            (setq p1 (red p1) p2 (red p2))))))
 125
 126 #-nil
 127 (declare-top(unspecial tlist))