share/affine/polybas.lisp

   1 ;;; -*- mode: lisp; package: cl-maxima; syntax: common-lisp -*-
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;                                                                    ;;;;;
   4 ;;;     Copyright (c) 1984 by William Schelter,University of Texas     ;;;;;
   5 ;;;     All rights reserved                                            ;;;;;
   6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   7
   8 (in-package :maxima)
   9
  10 (defvar *varlist* nil)
  11 (defvar *genvar* nil)
  12
  13 (defun affine-polynomialp (x)
  14   (cond  ((numberp x))
  15          ((atom x) nil)
  16          ((and (atom (car x)) (symbolp (car x)) (get (car x) 'disrep)))
  17          (t nil)))
  18
  19 (Defun header-poly (expr)
  20    (cond ((numberp expr) expr)
  21          ((and (numberp (car expr))
  22                (numberp (cdr expr)))
  23           (cond ((zerop (car expr))0)
  24                 (t (list '(rat simp) (car expr) (cdr expr)))))
  25          ((affine-polynomialp expr)(cons (list 'mrat 'simp *varlist* *genvar*) (cons expr 1)))
  26          ((rational-functionp expr)(cons (list 'mrat 'simp *varlist* *genvar* ) expr ))
  27          (t expr)))
  28
  29 (defun $zerop ( n &aux type-of-n tem )
  30   (cond ((numberp n) (zerop n))
  31         ((atom n) nil)
  32         ((affine-polynomialp n) nil)
  33         ((rational-functionp n)(rzerop   n))
  34         (($bfloatp n) (eql (car n) 0))
  35         (t (setf type-of-n (caar n))
  36            (cond ((member type-of-n '(mrat rat) :test #'eq)
  37                   (equal (cdr n) (rzero)))
  38                  (t (and (numberp (setq tem ($ratsimp n)))(zerop tem)))))))
  39
  40 (defun rational-functionp (x)
  41   (cond ((numberp x))
  42         ((atom x) nil)
  43         (t (and (affine-polynomialp (car x))
  44                 (affine-polynomialp (cdr x))))))
  45
  46 (defun psublis (a-list denom poly &key degree vars-to-sub)
  47   "does a general sublis : a-list of form (list (cons old-var repl-poly)....)
  48    denom is a poly "
  49   (let ((tem (cond (vars-to-sub vars-to-sub)
  50                    (t (loop for v in a-list collecting (car v)))))
  51         deg)
  52     (cond ((affine-polynomialp poly)
  53            (cond (degree (setq deg degree))
  54                  (t (setq deg (poly-degree poly tem))))
  55            (psublis1 a-list denom poly
  56                      deg  tem))
  57           ((rational-functionp poly) (rsublis a-list denom poly :degree degree :vars-to-sub
  58                                               vars-to-sub :reduce t))
  59           ((atom poly) poly)
  60           (t (loop for v in poly collecting (psublis a-list denom v :degree degree
  61                                                    :vars-to-sub vars-to-sub))))))
  62
  63 ;;should take into account when the variables don't need replacing.
  64 (defun psublis1 (a-list denom poly degree varl)
  65    (cond ((atom poly) (ptimes poly (pexpt denom degree)))
  66         ((member (p-var poly) varl :test #'eq)
  67          (loop for (deg cof) on (cdr poly) by #'cddr
  68                with answer = 0
  69                do (setq answer
  70                         (pplus answer
  71                                (ptimes
  72                                         (psublis1 a-list denom cof
  73                                                       (- degree deg) varl)
  74                                        (pexpt (cdr (assoc (p-var poly) a-list :test #'equal))
  75                                               deg))))
  76
  77                finally (return answer)))
  78         ((> (valget (p-var poly))
  79              (loop for v in varl minimize (valget v)))
  80          (loop for (deg cof) on  (cdr poly ) by #'cddr
  81                with answer = 0
  82                with mon = (list (p-var poly ) 1 1)
  83                do (setq answer (pplus answer
  84                                       (ptimes
  85                                        (pexpt mon deg)
  86                                        (psublis1 a-list denom cof degree varl))))
  87
  88                finally (return answer)))
  89         (t (ptimes (pexpt denom degree) poly))))
  90
  91 (defun zero-sublis (poly &rest list-vars)
  92   (pcoeff poly 1 list-vars))
  93
  94 (defun pcomplexity  (poly)
  95   (cond ((atom poly) 0)
  96         (t (loop for (deg cof) on (cdr poly) by #'cddr
  97                  summing (+ deg (pcomplexity cof))))))
  98
  99 (defun $numerator (expr)
 100   (cond ((atom expr) expr)
 101         ((mbagp expr) (cons (car expr) (mapcar '$numerator (cdr expr))))
 102         (($ratp expr) (cons (car expr) (cons (car (cdr expr)) 1)))
 103         (($totaldisrep (header-poly (cons (num (new-rat expr)) 1))))))
 104
 105 (defmacro allow-rest-argument (new-funct binary-op  answer-for-null-argument rest-arg)
 106   `(cond  ((null ,rest-arg) ,answer-for-null-argument)
 107           (t (case (length ,rest-arg)
 108                (,2 (,binary-op (first ,rest-arg)(second ,rest-arg)))
 109                (,1 (car ,rest-arg))
 110                (otherwise (apply ',new-funct (,binary-op (first ,rest-arg)
 111                                                          (second ,rest-arg)) (cddr ,rest-arg)))))))
 112
 113 (defmacro polyop (x y identity-x identity-y number-op poly-op rat-op &optional rat-switch )
 114   (cond (rat-switch (setq rat-switch '(t))))
 115   `(cond
 116      ((and (numberp ,x)(numberp ,y))(,number-op ,x ,y))
 117      ((eq ,x ,identity-x) ,y)
 118      ((eq ,y ,identity-y) ,x)
 119      (t
 120       (let ((xx (poly-type ,x)) answer
 121             (yy (poly-type ,y)))
 122
 123 ;       (cond ((or (eq xx '$rat)(eq yy '$rat))(with-vgp(check-vgp-correspond))));;remove later
 124         (cond ((null xx)(setq ,x (cdr ($new_rat ,x)) xx ':rational-function))
 125               ((eq xx ':$rat)
 126                (setq ,x (cdr ,x) xx ':rational-function))
 127               ((eq xx ':rat ) (setq ,x (cons (second ,x) (third ,x))
 128                                    xx ':rational-function)))
 129
 130         (cond
 131           ((null yy)(setq ,y (cdr ($new_rat ,y)) yy ':rational-function))
 132           ((eq yy ':$rat)(setq ,y (cdr ,y) yy ':rational-function))
 133           ((eq yy ':rat ) (setq ,y (cons (second ,y) (third ,y))
 134                                yy ':rational-function)))
 135         (cond ((and (eq xx ':rational-function) (eql (denom ,x) 1))
 136                (setq ,x (car ,x) xx :polynomial)))
 137         (cond ((and (eq yy ':rational-function) (eql (denom ,y) 1))
 138                (setq ,y (car ,y) yy :polynomial)))
 139         (setq answer
 140               (case xx
 141                 (:number   (case yy
 142                              (:number (,number-op ,x ,y))
 143                              (:polynomial (,poly-op ,x ,y))
 144                              (:rational-function
 145                               (,rat-op (cons ,x 1) ,y ,@ rat-switch))))
 146                 (:polynomial
 147                  (case yy
 148                    (:number (,poly-op ,x ,y))
 149                    (:polynomial (,poly-op ,x ,y))
 150                    (:rational-function (,rat-op (cons ,x 1) ,y,@ rat-switch))
 151                    (otherwise (merror "unknown type for polyop "))))
 152                 (:rational-function
 153                  (case yy
 154                    (:number (,rat-op ,x  (cons ,y 1) ,@ rat-switch))
 155                    (:polynomial (,rat-op ,x (cons ,y 1)  ,@ rat-switch))
 156                    (:rational-function (,rat-op ,x ,y ,@ rat-switch))))
 157                 (otherwise (merror "unknown arg"))))
 158         (cond ((affine-polynomialp answer) answer)
 159               ((rational-functionp answer)
 160                (cond ((eq 1 (cdr answer))(car answer))
 161                      (t answer)))
 162               (t answer))))))
 163
 164 (defun n* (&rest l)
 165   (case (length l)
 166     (0 1)
 167     (1 (car l))
 168     (2 (n1* (car l) (second l)))
 169     (t (n1* (car l) (apply 'n* (cdr l))))))
 170
 171 (defun n+ (x y)
 172   (polyop x y  0 0 + pplus ratplus ))
 173
 174 (defun n1* (x y)
 175   (polyop x y 1 1 * ptimes rattimes t))
 176
 177 (defun n- (x y)
 178   (polyop x y nil 0 - pdifference ratdifference))
 179
 180 (defun nred (x y &aux answer)
 181   (setq answer (polyop x y nil 1  ratreduce ratreduce ratquotient))
 182   (setq answer (rationalize-denom-zeta3 answer))
 183   (cond ((numberp answer) answer)
 184         ((eql (denom answer) 1) (car answer))
 185         (t answer)))
 186
 187 (defun new-disrep (expr)
 188   (cond ((atom expr) expr)
 189         (t
 190          (let ((type (poly-type expr)))
 191            (case type
 192              (:number expr)
 193              (:polynomial ($ratdisrep (header-poly expr)))
 194              (:rational-function
 195               (cond (($numberp expr)
 196                      (cond ((zerop (num expr)) 0)
 197                            (t
 198                      (list '(rat simp) (num expr) (denom expr)))))
 199                     (t
 200                      ($ratdisrep (header-poly expr)))))
 201              (otherwise  (cond ((mbagp expr)(cons (car expr) (mapcar 'new-disrep expr)))
 202                                (($ratp expr)($ratdisrep expr))
 203                                (t expr))))))))
 204  (defun poly-degree (poly varl)
 205   (cond ((atom poly) 0)
 206         ((member (p-var poly) varl :test #'eq)
 207          (loop for (deg cof) on (cdr poly) by #'cddr
 208                maximize (+ deg (poly-degree cof varl))))
 209         (t
 210          (loop for (deg cof) on (cdr poly) by #'cddr
 211                maximize (poly-degree cof varl)))))
 212
 213 (defun shl (l)
 214   (cond ($display2d (mapcar 'sh l))
 215        (t (loop for v in l
 216                 for i from 0
 217                 initially (format t "~%[")
 218                 when (numberp v) do (princ v)
 219                 else
 220                 do  (sh (header-poly v))
 221                 when (< i (1- (length l))) do(format t ",~%")
 222                 finally (format t "]")))))
 223
 224 (defun sh (expr &optional (stream *standard-output*) &aux (disp 'string-grind) answ)
 225   (declare (special *fake-rat*))
 226   (cond ( $display2d (setq disp 'displa)))
 227   (let ((*standard-output* stream))
 228     (setq answ(cond ((numberp expr) (funcall disp expr))
 229                     ((affine-polynomialp expr)(funcall disp (cons  *fake-rat* (cons expr 1))))
 230                     ((rational-functionp expr)(funcall disp (cons *fake-rat* expr)))
 231                     (t (funcall disp expr))))
 232     (cond (answ (format t "~A" answ)))))
 233
 234 (defmacro setq-num-den (num den expr &aux form .expr.)
 235   (cond ((symbolp expr) (setq .expr. expr))
 236         (t (setq .expr. '.expr.)))
 237   (setq form
 238         `(cond ((affine-polynomialp  , .expr.) (setq ,num , .expr. , den  1))
 239                ((rational-functionp , .expr.) (setq ,num (num , .expr.) ,den (denom , .expr.)))
 240                (t (fsignal "expr is supposed to be poly or rational-functionp"))))
 241   (cond ((symbolp expr)
 242          form)
 243         (t `(let ((.expr. ,expr))
 244               ,form))))
 245
 246 (defun splice-in (after-nth item a-list )
 247   (cond ((eql after-nth -1)(cons item a-list))
 248         (t (nconc (subseq a-list 0 (1+ after-nth))
 249                   (cons item (cdr (nthcdr after-nth a-list)))))))
 250
 251
 252 (defun nsplice-in (after-nth item a-list &aux tem)
 253   (cond ((eql after-nth -1)(cons item a-list))
 254         (t
 255 ;        (rplacd (setq tem (nthcdr after-nth a-list))
 256 ;                (cons item (cdr tem)))
 257          (setf (cdr (setq tem (nthcdr after-nth a-list)))
 258                  (cons item (cdr tem)))
 259
 260          a-list)))
 261
 262 (defun nc-equal (p1 p2)
 263   (or (and  (atom p1) (eql p1 p2))
 264       (and (listp p1) (listp p2) (equal (cdr p1 ) (cdr p2)))))
 265
 266 (defvar *check-order* nil)
 267
 268 (defun find-in-ordered-list (x a-list &optional  (order-function (function alphalessp))
 269                              (order-equal #'nc-equal)
 270                              &aux (offset 0)
 271                              (leng 0) tem (after 0)
 272                              (mid 0) (ub 0) (lbound 0) already-there)
 273   (declare (fixnum leng mid ub lbound after))
 274   (setq leng (length a-list))
 275
 276   (cond ((funcall order-function x (car a-list))(setq after -1))
 277         ((funcall order-equal  x (car a-list))(setq after 0)(setq already-there a-list))
 278         ((funcall order-equal  x (car(setq tem  (last a-list))))(setq after (1- leng))
 279          (setq already-there tem))
 280         ((not (funcall order-function x (car tem)))(setq after (1- leng)))
 281
 282         (t (setq offset 0 lbound 0 ub (1- leng))
 283            (loop while (> (- ub lbound) 1)
 284                  do
 285                  (setq mid (truncate (+ ub lbound) 2))
 286                  (cond ((funcall order-function x (car(setq tem
 287                                                         (nthcdr mid a-list))))(setq ub mid))
 288                        ((funcall order-equal  (car tem) x)(setq after (+ offset mid)
 289                                            already-there tem)(return 'done))
 290                        (t (setq offset (+ mid offset)
 291                                 a-list tem
 292                                 ub (- ub mid))))
 293                  finally(setq after offset))))
 294   (and *check-order*
 295        (cond ((setq tem (member x a-list :test order-equal))
 296               (iassert (eql (not already-there) (not t)))
 297               (iassert (eql after (- leng (length tem )))))))
 298   (values after already-there a-list))
 299
 300 (defun disrep-list (expr)
 301   "Just substitutes the any disrep property for a symbol for the symbol"
 302   (cond ((atom expr)(cond ((symbolp expr)(get expr 'disrep))
 303                           (t expr)))
 304                 (t (cons (disrep-list (car expr)) (disrep-list (cdr expr))))))
 305
 306 ;;Perhaps we should have our own prep1 function
 307 ;;It just goes through a general form converting it to cre, and
 308 ;;when it comes to a symbol it can't find in genpairs it
 309 ;;calls newsym which is much like our add-newvar . We could
 310 ;;have it look things up in a hash table rather than genpairs
 311 ;;The checking for being on tellratlist is done by the ratsetup2
 312 ;;One has to be careful with the with-vgp form since nested ones
 313 ;;will cause havoc.  Also if inside one you must refer to varlist
 314 ;;not *varlist* yet outside you may want *varlist*.  So for example
 315 ;;in
 316
 317
 318 (defun check-repeats (varl)
 319   (loop for v on varl
 320         when (member (car v) (cdr v) :test 'nc-equal)
 321           do (show (car v) ) (fsignal "repeat in varlist")))
 322
 323 (defun check-order (varl)
 324   (declare (special $order_function))
 325   (loop for (v w) on varl
 326         while w
 327           when (not (funcall $order_function  v w))
 328             do (fsignal "bad order ~A and ~A "  v w )))
 329
 330
 331 (defun generate-new-variable-gensym (variable  &aux tem)
 332        (setq tem (gensym))
 333                        (putprop tem variable 'disrep)
 334                        (push (cons variable (rget tem)) genpairs )
 335                     tem   )
 336
 337 (defmacro with-vg (&rest body)
 338   `(let ((varlist *varlist*)
 339          (genvar *genvar*))
 340      ,@body
 341      (setq *varlist* varlist)
 342      (setq *genvar* genvar)))
 343
 344 (defmacro with-vgp (&rest body)
 345   `(let ((varlist *varlist*)
 346         (genvar *genvar*)
 347         (genpairs *genpairs*))
 348         (prog1 (progn ,@ body)
 349         (setq *varlist* varlist *genpairs* genpairs)
 350         (setq *genvar* genvar))))
 351
 352 (defun zl-union (&rest l)
 353   (case (length l)
 354         (1 (loop for v on (car l) unless (member (car v) (cdr v))
 355                    collect (car v)))
 356         (2 (loop for v in (car l) with sec = (second l)
 357                   unless (member v sec)
 358                   collect v into result
 359                   finally (return (nconc result (zl-union sec)))))
 360         (t (zl-union (car l) (apply #'zl-union (cdr l))))))