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. 243