share/affine/sub-proj.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 (defun rational-subst (poly var denom )
  11   "does a substitution like var --> var/denom  where var is a genvar
  12 and denom is a poly.  It calculates the numerator after bringing a common
  13 denominator of denom^(pdegree poly var)"
  14    (let ((deg  (pdegree poly var)))
  15      (setq poly (copy-tree poly))
  16      (sub-rat2 poly var denom deg)))
  17
  18 (cond ((not (fboundp 'p-var)) (load "msm:ratmac.fasl")))
  19
  20 (defmacro with-main-variable (mvar &body body)
  21   `(let ((,mvar (make-symbol "zzzzzz")))
  22      (setf (symbol-value ,mvar) 40000)
  23      (setf (get ,mvar 'disrep) 'zzzzzzz)
  24      ,@ body))
  25
  26 (defmacro with-main  (main form-to-reorder resubstitute &body body)
  27   (cond ((null resubstitute)
  28          `(let ((newvar (make-symbol "zzzzzz"))
  29                 (oldmain ,main))
  30             (setf (symbol-value newvar) 40000)
  31             (setq ,main newvar)
  32             (setf (get newvar 'disrep) 'zzzzzzz)
  33             (setf ,form-to-reorder (psublis (list (cons oldmain (list newvar 1 1))) 1
  34                                              ,form-to-reorder))
  35             (progn ,@ body)))
  36         (t
  37          `(let ((newvar (make-symbol "zzzzzz"))
  38                 (oldmain ,main))
  39             (setf (symbol-value newvar) 40000)
  40             (setq ,main newvar)
  41             (setf (get newvar 'disrep) 'zzzzzzz)
  42             (setf ,form-to-reorder
  43                   (psublis (list  (cons oldmain (list newvar 1 1))) 1
  44                             ,form-to-reorder))
  45             (psublis  (list  (cons newvar (list oldmain 1 1))) 1
  46                       (progn ,@ body))))))
  47
  48
  49 (defmacro sum-over-polynomial (deg cof poly summand)
  50   `(do ((-tail- (cdr ,poly) (cddr -tail-))
  51         (-answer- 0)(,deg)(,cof))
  52        ((null -tail-) -answer-)
  53     (setq ,cof (cadr -tail-))
  54     (setq ,deg (car -tail-))
  55     (setq -answer- (pplus -answer- ,summand))))
  56 (defmacro replace-polynomial-coefficients (deg cof poly replacement)
  57   `(do ((-tail- (cdr ,poly) (cddr -tail-))
  58        (,deg)(,cof))
  59       ((null -tail-) ,poly)
  60     (setq ,cof (cadr -tail-))
  61     (setq ,deg (car -tail-))
  62     (setf (cadr -tail-) ,replacement)))
  63
  64 (defun $degree (form var)
  65   (setq form (st-rat form))
  66   (pdegree  form (car (st-rat var))))
  67
  68 (defun sub-rat2 (poly var denom deg)
  69   (cond ((atom poly) (setq poly (ptimes poly (pexpt denom deg))))
  70         ((eq (car poly) var)
  71          (cond ((pointergp var (car denom))
  72                 (replace-polynomial-coefficients
  73                   degree cof poly
  74                   (ptimes cof
  75                           (pexpt denom (- deg degree)))))
  76                (t
  77                 (setq poly
  78                 (sum-over-polynomial degree cof poly
  79                   (ptimes (list (car poly) degree 1)
  80                           (ptimes cof (pexpt denom (- deg degree)))))))))
  81         ((pointergp (car poly) var)
  82          (cond ((pointergp (car poly) (car denom))
  83                 (print 'a)
  84                 (replace-polynomial-coefficients
  85                   degree cof poly
  86                   (sub-rat2 cof var denom deg)))
  87                (t (print 'b)
  88                 (setq poly
  89                 (sum-over-polynomial
  90                   degree cof poly
  91                   (ptimes (list (car poly) degree 1)
  92                           (sub-rat2 cof var denom deg)))))))
  93         (t(setq poly (ptimes poly (pexpt denom deg)))))
  94   poly)
  95
  96 (defun $sub1 (poly var denom &optional deg)
  97   deg
  98   (setq poly (st-rat poly))
  99   (setq var (car (st-rat var)))
 100   (setq denom (st-rat denom))
 101   (header-poly (rational-subst poly var denom )))
 102
 103
 104 (defun poly-subst (poly var repl)
 105   "This substitutes for var (a genvar) the repl (a poly) staying in cre form"
 106 ;;  (cond   (not (and (symbolp var)   ;error check
 107   (poly-subst1 (copy-tree poly) var repl))
 108
 109
 110 (defun poly-subst1 (poly var repl)
 111         (cond ((atom poly) poly)
 112               ((eq (car poly) var)
 113                (do ((tail (cdr poly)(cddr tail))
 114                     (answer 0))
 115                    ((null tail) (setq poly answer))
 116                  (setq answer (pplus
 117                                (ptimes (pexpt repl(car tail))
 118                                        (cadr tail))
 119                                answer))))
 120               ((pointergp (car poly) (car repl))
 121                (do ((tail (cdr poly)(cddr tail)))
 122                    ((null tail))
 123                  (setf (cadr tail)
 124                        (poly-subst (cadr tail) var repl))))
 125               (t (cond
 126                   ((eql 0 (pdegree poly var)) nil)
 127                   (t
 128                    (do ((tail (cdr poly)(cddr tail))
 129                         (answer 0))
 130                        ((null tail) (setq poly answer))
 131                      (setq answer
 132                            (pplus answer
 133                                   (ptimes
 134                                    (list (car poly) (car tail) 1)
 135
 136                                    (poly-subst (cadr tail) var repl)))))))))
 137 poly)
 138
 139 (defun $polysub (poly var repl)
 140   (setq poly (st-rat poly))
 141   (setq var (car (st-rat var)))
 142   (setq repl (st-rat repl))
 143   (header-poly (poly-subst poly var repl )))
 144
 145 (defun rsublis (a-list denom poly &key degree vars-to-sub reduce &aux answ)
 146   "does a general sublis : a-list of form (list (cons old-var repl-poly)....)
 147    denom is a poly "
 148   (let ((tem (cond (vars-to-sub vars-to-sub)
 149                    (t (loop for v in a-list collecting (car v)))))
 150         deg)
 151     (cond ((affine-polynomialp poly)
 152            (cond (degree (setq deg degree))
 153                  (t(setq deg (poly-degree
 154                                poly
 155                                tem ))))
 156            (setq answ (psublis1 a-list denom poly
 157                                 deg
 158                                 tem))
 159            (cond (reduce (ratreduce  answ  (pexpt denom deg)))
 160                  (t (cons answ  (pexpt denom deg)))))
 161           ((rational-functionp poly)
 162            (ratquotient (rsublis a-list denom (num poly) :degree degree :vars-to-sub
 163                                  vars-to-sub )
 164                         (rsublis a-list denom (denom poly) :degree degree :vars-to-sub
 165                                  vars-to-sub)))
 166           (t (merror "bad type for poly : should be ratl fn or poly")))))
 167
 168 (defun $psublis (a-list denom poly)
 169   "use psublis([y=x^2,v=u^3],denom,poly)"
 170   (header-poly(psublis  (loop for (u v  repl)  in (cdr a-list) by  #'cdr
 171                               do (check-arg u (eq (car u) 'mequal) "Type a=repl")
 172                               collecting
 173                               (cons (p-var (st-rat v))
 174                                     (st-rat repl)))
 175                         (st-rat denom) (st-rat poly))))
 176
 177 (defun $coll_linear (expr &aux answ)
 178   (cond ((mbagp expr)
 179          (setq answ (loop for v in (cdr expr) collecting (coll-linear (st-rat v)))))
 180         (t (setq answ (coll-linear expr))))
 181   (setq answ (apply 'append answ))
 182   (loop for w in answ
 183         collecting (get w 'disrep) into tem
 184         finally (return (cons '(mlist) tem))))
 185
 186 (defun psubst ( repl var poly)
 187   (psublis (list (cons var repl)) 1 poly))
 188 (defun  pdiscriminant (poly var &aux main old-var answ)  ;;main variable
 189    ;;change to main variable if necessary
 190   (cond ((eq var :main)(setq var (p-var poly)))
 191         ((not (eql var (p-var poly)))
 192          (setq main (gensym))
 193          (setf (symbol-value main) 40000000)
 194          (setf (get main 'disrep) 'main-var)
 195          (show main)(setq old-var var)
 196          (setq poly (psubst (list main 1 1) var poly))
 197          (print 'hi)
 198          (setq var main)))
 199   (show poly)
 200    (mshow poly)
 201    (setq answ   (resultant poly (pderivative poly var) ))
 202    (cond (main
 203           (setq answ (psublis (list (cons main (list old-var 1 1))) 1 answ)))
 204          (t answ)))
 205
 206 (defun presultant (p1 p2 var)
 207   (cond ((and (consp p1) (consp p2)(eq (p-var p1) (p-var p2)) (eq var (p-var p1)))
 208          (resultant p1 p2 ))
 209         (t
 210          (let ((args (list p1 p2)))
 211            (with-main var  args t
 212              (resultant (car args) (second args)))))))
 213
 214
 215
 216 ;;incorrect needs
 217 (defun gen-vrem (b divisor &aux tem)
 218   (cond
 219     ((atom divisor) 0)
 220     ((atom b)b )
 221     ((pointergp (p-var divisor) (p-var b))
 222      b)
 223     ((eq (p-var divisor) (p-var b))
 224          (second (vdivide b divisor)))
 225     ((pointergp (p-var b) (p-var divisor))
 226      (loop for (deg cof) on (cdr b) by #'cddr
 227            do (setq tem (gen-vrem  cof divisor))
 228            when (not  (pzerop  tem))
 229            collecting deg into lis
 230            collecting tem into lis
 231            finally (return (cond ((eql (car lis) 0) (second lis))
 232                                  (t (cons (p-var b) lis))))))))
 233
 234 ;;incorrect doesn't take into account the
 235 ;;denoms eg z^3+ z*p1(x) +p2(x) must find the c for p1 and for p2 and then mult
 236 ;;up appropriateley.
 237 ;;won't work for x^2*y+x+1 where
 238 (defun gen-vrem (b divisor &aux tem)
 239   (cond
 240     ((atom divisor) 0)
 241     ((atom b)b )
 242     ((pointergp (p-var divisor) (p-var b))
 243      b)
 244     ((eq (p-var divisor) (p-var b))
 245          (second (vdivide b divisor)))
 246     ((pointergp (p-var b) (p-var divisor))
 247      (loop for (deg cof) on (cdr b) by #'cddr
 248            do (setq tem (gen-vrem  cof divisor))
 249            when (not  (pzerop  tem))
 250            collecting deg into lis
 251            collecting tem into lis
 252            finally (return (cond ((eql (car lis) 0)
 253                                   (second lis))
 254                                  (t (cons (p-var b) lis))))))))
 255
 256 ;;for resubstitute = nil then we would not want to
 257 ;(defun pdeg (pol var)
 258 ;  (with-main var pol t
 259 ;     (print    (second pol))(print var) pol))
 260 ;
 261 ;(defun pdeg (pol var)
 262 ;  (with-main var pol nil
 263 ;     (print    (second pol)) pol))
 264 ;
 265
 266
 267
 268 (defun highest-power-dividing (f divisor &aux quot )
 269 ;  (declare (values power final-quotient))
 270   (loop for i from 1
 271         with prev-quot = f
 272         do
 273         (setq quot (testdivide prev-quot divisor))
 274         when (null quot)
 275         do (return (values (1- i) prev-quot))
 276         else do (setq prev-quot quot)))
 277
 278 ;c-reqd*f=g*quot+remaind
 279 (defun gen-prem (f g var &aux remainder c-reqd)
 280   (cond ((< (pdegree f var) (pdegree g var))
 281          (setq remainder f)(setq c-reqd 1))
 282         (t
 283          (cond ((and (eq (p-var f) (p-var g))
 284                      (eq (p-var f) var))
 285                  (setq remainder  (vdivide f g))
 286                  (setq c-reqd (third remainder))
 287                  (setq remainder (second remainder)))
 288                (t (with-main-variable u
 289                     (setq f(psublis (list (cons var (list u 1 1))) 1 f))
 290                     (setq g(psublis (list (cons var (list u 1 1))) 1 g))
 291                     (setq remainder  (vdivide f g))
 292                     ;;this should not involve main variable..
 293                     (setq c-reqd (third remainder))
 294                     (setq remainder (second remainder))
 295                    (setq remainder (psublis (list (cons u
 296                                                         (list var 1 1)))
 297                                             1  remainder)))))))
 298   (values remainder c-reqd))
 299
 300 ;
 301 ;;;tests that (gen-prem  f g var) does same as numerator(remainder (f,g,var))
 302 ;;;it was 25 times faster .
 303 ;(defun $test (f g var &aux ans2 ans1 dis)
 304 ; (user:tim   (setq ans1 ($numerator ($remainder f g var))))
 305 ; (user:tim (setq f (st-rat f)) (setq g (st-rat g)))
 306 ;  (setq ans2 (gen-prem  f  g (add-newvar var)))
 307 ;  (setq dis (new-disrep ans2))
 308 ;;  (setq dis ($totaldisrep (header-poly ans2)))
 309 ;   (list '(mlist) ans1 dis ($ratsimp (sub* ans1 dis))))
 310
 311 ;;works now c*b=q*divisor +r
 312 (defun vdivide (b divisor &aux rnew rfactor leading-gcd deltaq)
 313   (let ((q 0)(c 1)(r  b))
 314     (cond
 315       ((atom divisor) (list b 0 divisor)) ;;should b/gcd(c(b),divisor) 0 a/same
 316       ((atom b) (list 0 b 1))
 317       ((pointergp (p-var divisor) (p-var b))
 318        (list 0 b 1))
 319       ((pointergp (p-var b) (p-var divisor))
 320        (list b 0 divisor))
 321       (t (loop until (or (atom r)
 322                          (not (eql (p-var r) (p-var divisor)))
 323                          (< (p-le r) (p-le divisor)))
 324                with mon = (list (p-var r) 1 1)
 325                do
 326                (setq leading-gcd (pgcd (p-lc r) (p-lc divisor)))
 327                (setq rfactor (pquotient (p-lc divisor) leading-gcd))
 328                ;;make r so that you can form r-divisor*q ==> lower degree
 329                (setq rnew (ptimes rfactor r))
 330                (setq c (ptimes c rfactor))
 331                ;;want rnew-divisor*deltaq to be lower degre in main variable
 332                ;;so need (p-lc deltaq) = (p-lc rnew)/(p-lc divisor)
 333                ;;but (p-lc rnew)=(p-lc divisor) (p-lc r) /leading-gcd
 334                ;;so (p-lc r)/leading-gcd will be the right (p-lc deltaq)
 335                ;;then multiply by (p-var r)^ m where m=(differenc of leading degrees)
 336                (setq deltaq (ptimes (pquotient (p-lc r) leading-gcd)
 337                                     (pexpt mon (- (p-le r) (p-le divisor)))))
 338 ;              (mshow  rnew divisor (p-lc rnew) (p-lc divisor) deltaq)
 339
 340 ;              (setq prev-rdeg (p-le r))
 341                (setq r (pdifference rnew
 342                                     (ptimes divisor
 343                                             deltaq)))
 344
 345 ;              (cond ((not (< (p-le r) prev-rdeg)) (break t)))
 346                (setq q (pplus (ptimes rfactor q) deltaq))
 347                finally (return (list q r c)))))))
 348
 349
 350
 351
 352 ;(defmacro mshow (&rest l)
 353 ;  (loop for v in l
 354 ;       collecting `(format t "~%The value of ~A is.. " ',v) into tem
 355 ;       collecting `(sh ,v) into tem
 356 ;       finally (return (cons 'progn tem))))
 357 ;
 358 ;(defun test ( b divisor)
 359 ;  (let ((answ (vdivide b divisor)))
 360 ;    (pdifference (ptimes (third answ) b)
 361 ;                (pplus (ptimes (first answ) a)
 362 ;                       (second answ)))))
 363
 364 (defun coll-linear  (f &aux *linear*)
 365         (declare (special *linear*))
 366         (coll-linear1 f)
 367         *linear*)
 368
 369 (defun constant-polyp (f)
 370   (cond ((atom f) t)
 371         (t (cond ((get (car f)'constant) t)
 372                  (t (do ((tail (cdr f) (cddr tail)))
 373                         ((null tail) t)
 374                       (cond ((null (constant-polyp (cadr f)))
 375                              (return nil)))))))))
 376
 377 ;;checked thisn
 378 (defun coll-linear1 (poly)
 379   (declare (special *non-linear*))
 380   (declare (special *linear*))
 381   ;;here we collect things that look linear but u*z+u  would collect u
 382   (cond ((atom poly) nil)
 383         ((get (car poly) 'constant)
 384          (loop for (deg cof) on (cdr poly) by #'cddr
 385                do (coll-linear1 cof)))
 386         (t (cond ((and (eql (p-le poly) 1)
 387                        (constant-functionp (p-lc poly)))
 388                   (pushnew (p-var poly) *linear*)))
 389            (cond ((eql 0 (nth (- (length poly)  2) poly))
 390                   (coll-linear1 (car (last poly)))))))
 391   ;;now must check these are really linear to remove the u collected above.
 392   (cond ((consp poly)
 393          (let ((ldeg 0))
 394          (cond ((get (p-var poly) 'constant)(setq ldeg 1)))
 395          (loop for u in *linear*
 396                do
 397          (loop for (deg cof) on  (cdr poly) by 'cddr while
 398                (> deg 0)
 399                when (> (pdegree cof u) ldeg)
 400                do  (setq *linear* (delete u *linear* :test #'equal)))))))
 401   )
 402
 403 ;(compare-functions
 404
 405 (defun gen-psublis (old new poly)
 406   (multiple-value-bind (subs denom)
 407       (subs-for-psublis old new)
 408     (psublis subs denom poly)))
 409
 410 (defun gen-rat-sublis (old new f &optional (switch t) &aux answ)
 411   (simple-rat-sublis (subs-for-simple-rat-sublis old new) f))
 412
 413 (defun subs-for-psublis (old new &optional &aux (lcd 1) a-denom)
 414   (cond (($listp old) (setq old (cdr old))))
 415   (cond (($listp new) (setq new (cdr new))))
 416   (loop for v in old when  (not (get v 'disrep))
 417      do (return (loop for v in old
 418                    when (not (get v 'disrep ))
 419                    do (cond ((symbolp v)
 420                              (setq v (add-newvar v)))
 421                             (t (merror "only for replacing symbols")))
 422                    collecting v into tem
 423                    finally (setq old tem))))
 424   (setq new (loop for v in new
 425                when (symbolp v)
 426                do (cond ((get v 'disrep)
 427                          (setq v (list v 1 1)))
 428                         (t (setq v (st-rat v))))
 429                else when (affine-polynomialp v) nconc nil
 430                else do (setq v (new-rat v)) (setq a-denom t)
 431                  (setq lcd (plcm lcd (denom v)))
 432                collecting v))
 433   (cond (a-denom
 434          (loop for v on new
 435             when (affine-polynomialp (car v)) do (setf (car v)  (gen-ptimes lcd (car v)))
 436             else
 437             do (setf (car v) (gen-ptimes (num (car v)) (pquotient lcd (denom (car v)))))))
 438         (t (setq new (subseq new 0 (length old)))))
 439   (loop for v in old
 440      for w on new
 441      do (setf (car w) (cons v (car w))))
 442   (values new lcd))
 443
 444 (defun simple-rat-sublis ( subs f &optional (switch t) &aux sub)
 445   (cond ((atom f) (cons f 1))
 446         ((and (numberp (car f))(numberp (cdr f))) f)
 447         ((affine-polynomialp f)
 448          (setq sub (cdr (assoc (p-var f) subs :test #'equal)))
 449          (cond ((null sub)(setq sub (cons  (list  (p-var f) 1 1) 1))))
 450          (loop for (deg cof) on (cdr f) by #'cddr
 451                with answ = (cons 0 1)
 452                do (setq answ (ratplus answ (rattimes (ratexpt sub deg)
 453                                                      (simple-rat-sublis subs cof) switch)))
 454                finally (return answ)))
 455         ((rational-functionp f)
 456          (ratquotient (simple-rat-sublis subs (num  f) switch)
 457                (simple-rat-sublis subs (denom  f) switch)))
 458         ((loop for v in f collecting (simple-rat-sublis subs v switch)))))
 459
 460
 461 (defun subs-for-simple-rat-sublis (varl subs)
 462   (loop for v in varl
 463         for w in subs
 464         when (not  (get v 'disrep))
 465         do (setq v (add-newvar v))
 466         do (cond ((affine-polynomialp w)
 467                   (setq w (cons w 1)))
 468                  ((rational-functionp w) nil)
 469                  ((get w 'disrep)(setq w (cons (list w 1 1)1)))
 470                  (t  (setq w (new-rat w))))
 471         collecting (cons v w)))
 472
 473
 474
 475
 476 (defun gen-rat-sublis (old new f &optional (switch t) &aux )
 477    (simple-rat-sublis (subs-for-simple-rat-sublis old new) f))
 478
 479 ;
 480 ;;;I tested it on a number of fns and it worked.
 481 ;(defun test (old new fn)
 482 ;  (setq ans1 (gen-psublis old new fn))
 483 ;  (setq ans1 (ratreduce (num ans1) (denom ans1)))
 484 ;  (setq ans2  (simple-rat-sublis (subs-for-simple-rat-sublis old new) fn))
 485 ;  (mshow ans2 ans1)
 486 ;  (equal ans1 ans2))
 487
 488 (defun collect-monomial-coefficients (poly variables)
 489   "Returns a list of all coefficients of all monomials in VARIABLES occurring in POLY"
 490   (cond ((null variables) (list poly))
 491         (($zerop poly) (list poly))
 492         (t (loop for i from 0 to (pdegree poly (car variables))
 493                  when (>= i  1)
 494                  appending (collect-monomial-coefficients
 495                              (pcoeff poly (list (car variables) i 1)) (cdr variables))
 496                  else
 497                  appending
 498                   (collect-monomial-coefficients
 499                     (pcoeff poly 1 (subseq variables 0 1)) (cdr variables))))))
 500
 501
 502 (defun collect-monomial-coefficients-and-monomials (poly variables &optional (monomial 1))
 503   "Returns a list of two element lists,each whose first element is the monomial and
 504    whose second is the coefficient.  See the associated function which collects and
 505    then verifies it has the whole sum"
 506
 507   (cond ((null variables) (list (list monomial poly)))
 508         (t (loop for i from 0 to (pdegree poly (car variables))
 509                  for mon = (list (car variables) i 1)
 510                  when (>= i  1)
 511                  appending (collect-monomial-coefficients-and-monomials
 512                              (pcoeff poly mon)
 513                              (cdr variables)
 514                              (ptimes monomial mon))
 515                  else
 516                  appending
 517                  (collect-monomial-coefficients-and-monomials
 518                    (pcoeff poly 1 (subseq variables 0 1)) (cdr variables) monomial)))))
 519
 520
 521 (defun collect-and-verify-coefficients (poly variables &aux answ)
 522   (setq answ (collect-monomial-coefficients-and-monomials poly variables))
 523   (loop for v in answ
 524         with answer = 0
 525         do (setq answer (pplus answer (ptimes (first v) (second v))))
 526         finally (iassert (equal answer poly))
 527         (return (mapcar 'cadr answ))))
 528
 529 (defun va (h)
 530   (list-variables  (st-rat h)))