src/rat3b.lisp

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