src/rat3b.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 rat3b)
  14
  15 ;;      THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 2.
  16 ;;      IT INCLUDES RATIONAL FUNCTIONS ONLY.
  17
  18 (defun ralgp (r) (or (palgp (car r)) (palgp (cdr r))))
  19
  20 (defun palgp (poly)
  21   (cond ((pcoefp poly) nil)
  22         ((alg poly) t)
  23         (t (do ((p (cdr poly) (cddr p)))
  24                ((null p))
  25              (and (palgp (cadr p)) (return t))))))
  26
  27
  28 (defun ratdx (e *x*)
  29   (declare (special *x*))
  30   (prog (varlist flag v* genvar *a a trunclist)
  31      (declare (special v* *a flag trunclist))
  32      (and (member 'trunc (car e) :test #'eq) (setq trunclist (cadddr (cdar e))))
  33      (cond ((not (eq (caar e) (quote mrat))) (setq e (ratf e))))
  34      (setq varlist (caddar e))
  35      (setq genvar (car (cdddar e)))
  36      ;; Next cond could be flushed if genvar would shrink with varlist
  37      (cond ((> (length genvar) (length varlist))
  38             ;; Presumably this produces a copy of GENVAR which has the
  39             ;; same length as VARLIST.  Why not rplacd?
  40             (setq genvar (mapcar #'(lambda (a b) (declare (ignore a)) b)
  41                                  varlist genvar))))
  42      (setq *x* (fullratsimp *x*))
  43      (newvar *x*)
  44      (setq a (mapcan #'(lambda (z)
  45                          (prog (ff)
  46                             (newvar
  47                              (setq ff (fullratsimp (sdiff z *x*))))
  48                             (orderpointer varlist)
  49                             (return (list z ff)))) varlist))
  50      (setq *a (cons nil a))
  51      (mapc #'(lambda(z b)
  52                (cond ((null (old-get *a z))(putprop b (rzero) 'diff))
  53                      ((and(putprop b(cdr (ratf (old-get *a z))) 'diff)
  54                           (alike1 z *x*))
  55                       (setq v*  b))
  56                      (t (setq flag t)))) varlist genvar)
  57
  58      ;;; causing lisp error - [ 2010843 ] diff of Taylor poly
  59      ;;(cond ((and (signp n (cdr (old-get trunclist v*)))
  60      ;;          (car (old-get trunclist v*))) (return 0)))
  61
  62      (and trunclist
  63           (return (cons (list 'mrat 'simp varlist genvar trunclist 'trunc)
  64                         (cond (flag (psdp (cdr e)))
  65                               (t (psderivative (cdr e) v*))))))
  66      (return (cons (list 'mrat 'simp varlist genvar)
  67                    (cond (flag (ratdx1 (cadr e) (cddr e)))
  68                          (t (ratderivative (cdr e) v*)))))))
  69
  70 (defun ratdx1 (u v)
  71   (ratquotient (ratdif (rattimes (cons v 1) (ratdp u) t)
  72                        (rattimes (cons u 1) (ratdp v) t))
  73                (cons (pexpt v 2) 1)))
  74
  75 (defun ratdp (p)
  76   (cond ((pcoefp p) (rzero))
  77         ((rzerop (get (car p) 'diff))
  78          (ratdp1 (cons (list (car p) 'foo 1) 1) (cdr p)))
  79         (t (ratdp2 (cons (list (car p) 'foo 1) 1)
  80                    (get (car p) 'diff)
  81                    (cdr p)))))
  82
  83 (defun ratdp1 (x v)
  84   (cond ((null v) (rzero))
  85         ((equal (car v) 0) (ratdp (cadr v)))
  86         (t (ratplus (rattimes (subst (car v) 'foo x) (ratdp (cadr v)) t)
  87                     (ratdp1 x (cddr v))))))
  88
  89 (defun ratdp2 (x dx v)
  90   (cond ((null v) (rzero))
  91         ((equal (car v) 0) (ratdp (cadr v)))
  92         ((equal (car v) 1)
  93          (ratplus (ratdp2 x dx (cddr v))
  94                   (ratplus (rattimes dx (cons (cadr v) 1) t)
  95                            (rattimes (subst 1 'foo x)
  96                                      (ratdp (cadr v)) t))))
  97         (t (ratplus (ratdp2 x dx (cddr v))
  98                     (ratplus (rattimes dx
  99                                        (rattimes (subst (1- (car v))
 100                                                         'foo
 101                                                         x)
 102                                                  (cons (ptimes (car v)
 103                                                                (cadr v))
 104                                                        1)
 105                                                  t)
 106                                        t)
 107                              (rattimes (ratdp (cadr v))
 108                                        (subst (car v) (quote foo) x)
 109                                        t))))))
 110
 111 (defun ratderivative (rat  var)
 112   (let ((num (car rat))
 113         (denom (cdr rat)))
 114     (cond ((equal 1 denom) (cons (pderivative num var) 1))
 115           (t (setq denom (pgcdcofacts denom (pderivative denom var)))
 116              (setq num (ratreduce (pdifference (ptimes (cadr denom)
 117                                                        (pderivative num var))
 118                                                (ptimes num (caddr denom)))
 119                                         ;RATREDUCE ONLY NEEDS TO BE DONE WITH CONTENT OF GCD WRT VAR.
 120                                   (car denom)))
 121              (cond ((pzerop (car num)) num)
 122                    (t (rplacd num (ptimes (cdr num)
 123                                           (pexpt (cadr denom) 2)))))))))
 124
 125 (defun ratdif (x y)
 126   (ratplus x (ratminus y)))
 127
 128 (defun ratfact (x fn)
 129   (cond ((and $keepfloat (or (pfloatp (car x)) (pfloatp (cdr x)))
 130               (setq fn 'floatfact) nil))
 131         ((not (equal (cdr x) 1))
 132          (nconc (funcall fn (car x)) (fixmult (funcall fn (cdr x)) -1)))
 133         (t (funcall fn (car x)))))
 134
 135 (defun floatfact (p)
 136   (destructuring-let (((cont primp) (ptermcont p)))
 137     (setq cont (monom->facl cont))
 138     (cond ((equal primp 1) cont)
 139           (t (append cont (list primp 1))))))
 140
 141 (defun ratinvert (y)
 142   (ratalgdenom
 143    (cond ((pzerop (car y)) (rat-error "`quotient' by `zero'"))
 144          ((and modulus (pcoefp (car y)))
 145           (cons (pctimes (crecip (car y)) (cdr y)) 1))
 146          ((and $keepfloat (floatp (car y)))
 147           (cons (pctimes (/ (car y)) (cdr y)) 1))
 148          ((pminusp (car y)) (cons (pminus (cdr y)) (pminus (car y))))
 149          (t (cons (cdr y) (car y))))))
 150
 151 (defun ratminus (x)
 152   (cons (pminus (car x)) (cdr x)))
 153
 154 (defun ratalgdenom (x)
 155   (cond ((not $ratalgdenom) x)
 156         ((pcoefp (cdr x)) x)
 157         ((and (alg (cdr x))
 158               (ignore-rat-err
 159                 (rattimes (cons (car x) 1)
 160                           (rainv (cdr x)) t))))
 161         (t x)))
 162
 163 (defun ratreduce (x y &aux b)
 164   (cond ((pzerop y) (rat-error "`quotient' by `zero'"))
 165         ((pzerop x) (rzero))
 166         ((equal y 1) (cons x 1))
 167         ((and $keepfloat (pcoefp y) (or $float (floatp y) (pfloatp x)))
 168          (cons (pctimes (quotient 1.0 y) x) 1))
 169         (t (setq b (pgcdcofacts x y))
 170            (setq b (ratalgdenom (rplacd (cdr b) (caddr b))))
 171            (cond ((and modulus (pcoefp (cdr b)))
 172                   (cons (pctimes (crecip (cdr b)) (car b)) 1))
 173                  ((pminusp (cdr b))
 174                   (cons (pminus (car b)) (pminus (cdr b))))
 175                  (t b)))))
 176
 177 (defun ptimes* (p q)
 178   (cond ($ratwtlvl (wtptimes p q 0))
 179         (t (ptimes p q))))
 180
 181 (defun rattimes (x y gcdsw)
 182   (cond ($ratfac (facrtimes x y gcdsw))
 183         ((and $algebraic gcdsw (ralgp x) (ralgp y))
 184          (let ((w  (rattimes x y nil)))
 185            (ratreduce (car w) (cdr w))))
 186         ((equal 1 (cdr x))
 187          (cond ((equal 1 (cdr y)) (cons (ptimes* (car x) (car y)) 1))
 188                (t (cond (gcdsw (rattimes (ratreduce (car x) (cdr y))
 189                                          (cons (car y) 1) nil))
 190                         (t (cons (ptimes* (car x) (car y)) (cdr y)))))))
 191         ((equal 1 (cdr y)) (rattimes y x gcdsw))
 192         (t (cond (gcdsw (rattimes (ratreduce (car x) (cdr y))
 193                                   (ratreduce (car y) (cdr x)) nil))
 194                  (t (cons (ptimes* (car x) (car y))
 195                           (ptimes* (cdr x) (cdr y))))))))
 196
 197 (defun ratexpt (x n)
 198   (cond ((equal n 0) '(1 . 1))
 199         ((equal n 1) x)
 200         ((minusp n) (ratinvert (ratexpt x (- n))))
 201         ($ratwtlvl (ratreduce (wtpexpt (car x) n) (wtpexpt (cdr x) n)))
 202         ($algebraic (ratreduce (pexpt (car x) n) (pexpt (cdr x) n)))
 203         (t (cons (pexpt (car x) n) (pexpt (cdr x) n)))))
 204
 205 (defun ratplus (x y &aux q n)
 206   (cond ($ratfac (facrplus x y))
 207         ($ratwtlvl
 208          (ratreduce
 209           (pplus (wtptimes (car x) (cdr y) 0)
 210                  (wtptimes (car y) (cdr x) 0))
 211           (wtptimes (cdr x) (cdr y) 0)))
 212         ((and $algebraic (ralgp x) (ralgp y))
 213          (ratreduce
 214           (pplus (ptimeschk (car x) (cdr y))
 215                  (ptimeschk (car y) (cdr x)))
 216           (ptimeschk (cdr x) (cdr y))))
 217         ((equal 1 (cdr x))
 218          (cond ((equal 0 (car x)) y)
 219                ((equal 1 (cdr y)) (cons (pplus (car x) (car y)) 1))
 220                (t (cons (pplus (ptimes (car x) (cdr y)) (car y)) (cdr y)))))
 221         ((equal 1 (cdr y))
 222          (cond ((equal 0 (car y)) x)
 223                (t (cons (pplus (ptimes (car y) (cdr x)) (car x)) (cdr x)))))
 224         (t (setq q (pgcdcofacts (cdr x) (cdr y)))
 225            (setq n (pplus (ptimes (car x)(caddr q))
 226                           (ptimes (car y)(cadr q))))
 227            (if (cadddr q)               ; denom factor from algebraic gcd
 228                (setq n (ptimes n (cadddr q))))
 229            (ratreduce n
 230                       (ptimes (car q)
 231                               (ptimes (cadr q) (caddr q)))))))
 232
 233 (defun ratquotient (x y)
 234   (rattimes x (ratinvert y) t))
 235
 236 ;;      THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 2.
 237 ;;      IT INCLUDES RATIONAL FUNCTIONS ONLY.