share/affine/polyc.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
  11 ;;(newvar expr) will add in a sorted manner all the new variables to the beginning
  12 ;;of the varlist.
  13
  14 (defvar *genpairs* nil)
  15
  16 (defun reset-vgp ()  (setq *genpairs* nil  *genvar*  nil  *varlist* nil)
  17        (setq *xxx*
  18              (let (*print-radix*)
  19                (loop for i from 1 to 30 collecting
  20                       (add-newvar (intern (format nil "$X~A" i))))))
  21        'done)
  22
  23
  24 (defun show-vgp ()
  25   (show *genpairs* genpairs *genvar* genvar *varlist* varlist))
  26
  27 (defvar vlist nil)
  28
  29 (defun put-tellrat-property (va the-gensym)
  30   (loop for v in tellratlist
  31         when (equal va (car v))
  32         do (putprop the-gensym (cdr v) 'tellrat)))
  33
  34
  35 (defun poly-ncmul1 (mon   poly  mon1 &aux monom new-monom tem gen-sym)
  36   (cond ((null mon)(setq mon 1))
  37         ((null mon1) (setq mon1 1)))
  38   (cond ((not (affine-polynomialp poly))(break t)))
  39
  40   (cond ((and (numberp mon)(numberp mon1))
  41          (n* poly (n* mon mon1)))
  42         ((or (numberp poly)($scalarp (setq monom (get (car poly) 'disrep))))
  43          (setq gen-sym (add-newvar (setq monom (ncmul* mon mon1))))
  44          (n* (list gen-sym 1 1) ;;(assolike monom *genpairs*);;might be long list to look in
  45              poly))
  46         (t (setq new-monom (ncmul* mon monom mon1))
  47            (cond ((contains-a-zero-replacement new-monom)
  48                   (cond ((and (eql (second poly) 1)(eql (fourth poly) 0))
  49                          (poly-ncmul1 mon (fifth poly) mon1))
  50                         ((eql (second poly) 1) 0)
  51                         (t (merror "There is a bad order in nc polynomial ~A" poly))))
  52                  (t
  53                   (setq gen-sym (add-newvar new-monom))
  54                   (cond ((and (eql (second poly) 1)(eql (fourth poly) 0))
  55                          (setq tem (poly-ncmul1 mon (fifth poly) mon1))
  56                          (cond
  57                            ((eql tem 0)(list gen-sym 1 (third poly)))
  58                            (t
  59                             (list gen-sym 1 (third poly) 0
  60                                    tem))))
  61                         ((eql (second poly) 1)
  62                          (list gen-sym 1 (third poly)))
  63                         (t (merror "There is a bad order in nc polynomial ~A" poly))))))))
  64
  65
  66 ;the above checks for terms in the radical and zero simplifications.
  67 ;(defun poly-ncmul1 (mon   poly  mon1 &aux monom new-monom gen-sym)
  68 ;  (cond ((null mon)(setq mon 1))
  69 ;       ((null mon1) (setq mon1 1)))
  70 ;  (cond ((not (affine-polynomialp poly))(break t)))
  71 ;
  72 ;  (cond ((and (numberp mon)(numberp mon1))
  73 ;        (n* poly (n* mon mon1)))
  74 ;        ((or (numberp poly)($scalarp (setq monom (get (car poly) 'disrep))))
  75 ;        (setq gen-sym (add-newvar (setq monom (ncmul* mon mon1))))
  76 ;        (n* (list gen-sym 1 1) ;;(assolike monom *genpairs*);;might be long list to look in
  77 ;            poly))
  78 ;       (t (setq new-monom (ncmul* mon monom mon1))
  79 ;
  80 ;          (setq gen-sym (add-newvar new-monom))
  81 ;          (cond ((and (eql (second poly) 1)(eql (fourth poly) 0))
  82 ;                 (list gen-sym 1 (third poly) 0 (poly-ncmul1 mon (fifth poly) mon1)))
  83 ;                ((eql (second poly) 1)
  84 ;                 (list gen-sym 1 (third poly)))
  85 ;                (t (merror "There is a bad order in nc polynomial ~A" poly))))))
  86 ;(defun new-rat-ncmul (a nc-rat-expr c)
  87 ;  (cons (poly-ncmul1 a (num nc-rat-expr) c) (denom nc-rat-expr)))
  88
  89 (defun poly-ncmul (&rest ll)
  90   "broken"
  91   (loop for v in ll when (not (affine-polynomialp v))do (break t))
  92   (cond ((null ll) 1)
  93         ((eql (length ll) 1) (car ll))
  94         ((numberp (car ll))(n* (car ll) (apply 'poly-ncmul (cdr ll))))
  95         ((affine-polynomialp (car ll))
  96               (apply
  97                                   'poly-ncmul (poly-ncmul1 1 (car ll) (second ll))
  98                                   (nthcdr 2 ll)))
  99         ((affine-polynomialp (second ll))
 100          (cond ((affine-polynomialp (third ll)) (merror "poly-ncmul can't have two 'affine-polynomialp objects in its arg yet"))
 101                (t (apply 'poly-ncmul
 102                                  (poly-ncmul1 (car ll)  (second ll) (third ll))
 103                                  (nthcdr 3 ll)))))))
 104
 105 (defun reset-gen () (setq genvar nil varlist nil genpairs nil))
 106
 107 (defun show-vg (&optional (varl *varlist* ) (genv *genvar*))
 108   (loop for v in varl
 109         for w in genv
 110         do (format t "~%~A ~A disrep ~A value ~A" w v (get w 'disrep) (symbol-value w))
 111         when (not (equal (get w 'disrep) v))do (format t" **** bad disrep ***")
 112         when (not (equal (get-genvar v) w)) do (format t" **** bad *genpair***~A" (assolike v *genpairs*)))
 113   (format t "~%~A and ~A are the lengths of genv and varl " (length genv) (length varl)))
 114
 115
 116 (defun last2 (a-list)
 117   (loop for v on a-list
 118         while (cddr v)
 119         finally (return v)))
 120
 121 (defun lastn (n a-list)
 122   (loop for v on a-list
 123         while (nthcdr n v)
 124         finally (return v)))
 125
 126
 127
 128 (defun choose (m n)
 129  (quotient (factorial m) (* (factorial n) (factorial (- m n)))))
 130
 131 (defun polynomial-ring-1-1-1  (deg)
 132         (choose (+ deg 2) deg))
 133
 134 (defun three-times-n (n) (* 3 n))
 135
 136 (defun list-equations1 (expr)
 137   "converts a single macsyma equation into a lisp list of cons's suitable for sublis"
 138   (cond ((and (listp (car expr))(eq (caar expr) 'mequal))(setq expr (cdr expr)))
 139         (t (merror "not a equation ~A " expr)))
 140   (cond ((atom (first expr))(list (cons (first expr) (second expr))))
 141         ((member (caar (first expr)) '($matrix mlist) :test #'eq)
 142          (loop for v in (cdr (first expr))
 143                for w in (cdr (second expr))
 144                appending (list-equations1 (list v w))))))
 145
 146 (defun list-equations-macsyma1 (expr)
 147   (cond ((and (listp (car expr))(eq (caar expr) 'mequal))
 148          (cond ((atom (second expr)) (list expr))
 149                ((not (member (caar (second expr)) '($matrix mlist) :test #'eq)) (list expr))
 150                (t (setq expr (cdr expr))
 151                   (cond ((member (caar (first expr)) '($matrix mlist) :test #'eq)
 152                          (cond ((or (atom (second expr))(not (eql (caar (first expr))
 153                                                             (caar (second expr)))))
 154                                 (loop for v in (cdr (first expr))
 155                                       appending (list-equations-macsyma1
 156                                                   (list '(mequal) v (second expr)))))
 157                                (t
 158                          (loop for v in (cdr (first expr))
 159                                for w in (cdr (second expr))
 160                                appending (list-equations-macsyma1
 161                                            (list '(mequal) v w))))))))))))
 162
 163
 164 (defun $list_factors(poly)
 165   (cond
 166     ((atom poly ) poly)
 167     ((mbagp poly) (cons (car poly) (mapcar #'$list_factors (cdr poly))))
 168     ((eq (caar poly) 'mtimes)
 169      (cons '(mlist) (cdr poly)))
 170     (t poly)))
 171
 172 (defun $list_equations (a-list)
 173   (check-arg a-list '$listp "macsyma list")
 174   (loop for v in (cdr a-list)
 175         appending (list-equations-macsyma1 v) into tem
 176         finally (return (cons '(mlist) tem))))
 177
 178 (defun $sub_list (eqns expr)
 179   (check-arg eqns '$listp "macsyma list")
 180   (setq eqns ($list_equations eqns))
 181   ($sublis eqns expr))
 182
 183 ;  (setq eqns (loop for v in (cdr eqns)
 184 ;       appending (list-equations1 v) ))
 185 ;  (new-sublis eqns expr))
 186
 187 (defun new-sublis (subs expr &aux rat-form answer)
 188   (cond (($ratp expr)(setq expr ($ratdisrep  expr)) (setq rat-form t)))
 189     (cond ((mbagp expr)(cons (car expr)(mapcar #'(lambda (x) (new-sublis subs x))
 190                                                (cdr expr))))
 191           (t (setq answer (sublis subs expr))(cond (rat-form ($new_rat answer))
 192                                                    (t (meval* expr))))))
 193
 194 (defun remove-from-*genpairs* (expr)
 195   (setq *genvar* (delete (get-genvar expr) *genvar* :test #'equal))
 196   (loop for v in *genpairs*
 197         for i from 0
 198         when (equal (car v) expr)
 199         do
 200         (cond ((eql i 0)(setq *genpairs* (cdr *genpairs*)))
 201               (t (setq *genpairs* (delete v *genpairs* :test #'equal)))))
 202   (setq *varlist* (delete expr *varlist* :test #'equal)))
 203
 204 ;
 205 ;(defmacro must-replacep (symbol-to-set rat-poly)
 206 ;  `(let ((.chang.))
 207 ;     (multiple-value (,symbol-to-set .chang.)
 208 ;       (rat-polysimp ,rat-poly))
 209 ;     .chang.))
 210
 211
 212 (defun pe (&aux tem)
 213     (format t "Enter a Macsyma expression:")
 214   (setq tem (mread-noprompt *standard-input*))
 215   (cdr ($new_rat (meval* tem))))
 216
 217 (defvar *fake-rat* '(mrat nil nil nil))
 218 ;
 219 ;(defun sh (expr)
 220 ;  (cond ((numberp expr)(displa expr))
 221 ;       ((affine-polynomialp expr)(displa (cons  *fake-rat* (cons expr 1))))
 222 ;       ((rational-functionp expr)(displa (cons *fake-rat* expr)))
 223 ;       (t (displa expr)))
 224 ;  expr)
 225
 226
 227 (defvar $homework_done '((mlist)))
 228 (defquote $homework ( a-list &aux answer work )
 229   (loop for v in (cdr a-list)
 230         when (symbolp v)
 231         do (setq $homework_done (append $homework_done (list  (setq work (meval* v)) )))
 232         else
 233         do (setq $homework_done  (append $homework_done (list  (setq work  v))))
 234         do (format t "~%Working on ...") (displa work)
 235         do  (setq $homework_done (append $homework_done (list  (setq answer(meval* work)))))
 236         (displa answer))
 237   $homework_done)
 238
 239 (defun $eij (n i j)
 240   ($setelmx 1 i j ($ident n)))
 241
 242 (defun $unipotent_invariant_vectors (representation n1 n2 &aux answer eqns)
 243   "representation should take two arguments a (size n1) and b (size n2) where a is
 244    in GLn1 and b is a macsyma list of length n2. It returns the general vector b left invariant  by the upper triangular (unipotent)  matrices."
 245   (loop for i from 1 below n1
 246      appending (cdr ($list_equations
 247                       (list '(mlist)
 248                             (list '(mequal)
 249                                   (mfuncall representation ($eij n1 i (1+ i))
 250                                             ($firstn n2 $aaaa))
 251                                   ($firstn n2 $aaaa))))) into tem
 252     finally (setq answer ($fast_linsolve (setq eqns (cons '(mlist) tem)) ($firstn n2 $aaaa))))
 253
 254   ($separate_parameters ($sublis answer ($firstn n2 $aaaa))))
 255 (defun $matrix_from_list (a-list &aux tem)
 256   (let ((n (round (setq tem (expt ($length a-list) .5)))))
 257     (cond ((eql (expt n 2) ($length a-list)) 'fine)
 258           (t (merror "The length of list is not a square")))
 259     (setq a-list (cdr a-list))
 260     (loop while a-list
 261           collecting (cons '(mlist) (subseq a-list 0 n)) into tem1
 262           do (setq a-list (nthcdr n a-list))
 263           finally (return  (cons '($matrix ) tem1)))))
 264
 265 (defun kronecker-row-product (a b)
 266   (loop for v in (cdr b)
 267         appending
 268         (loop for w in (cdr a)
 269                  collecting (mul* v w))
 270         into tem
 271         finally (return (cons '(mlist) tem))))
 272
 273 (defun $kronecker_product (m n)
 274   (loop for row in (cdr n)
 275         appending
 276         (loop for rowm in (cdr m)
 277               collecting (kronecker-row-product rowm row))
 278         into tem1
 279         finally (return (cons '($matrix) tem1))))
 280
 281 (defun $standard_rep (mat v)($list_matrix_entries(ncmul* mat v)))
 282
 283 (defun $unlist_matrix_entries (vect &optional size  )
 284   (cond (size nil)
 285         (t (setq size (round (expt  (length vect) .5)))))
 286   (setq vect (cdr vect))
 287   (loop for i below size
 288         collecting (cons '(mlist) (subseq vect 0 size)) into tem
 289         do (setq vect (nthcdr  size vect))
 290         finally (return (cons '($matrix) tem))))
 291
 292 (defun sub-list (a-list expr)
 293   (loop for v on a-list by 'cddr
 294         collecting `((mequal) ,(car v) ,(second v)) into tem
 295         finally (return ($sub_list (cons '(mlist) tem) expr))))
 296 (defmacro defrep (name argu &rest body &aux v mat)
 297    (setq mat (first argu) v (second argu))
 298   (setq body (sublis (list (cons v 'gen-vector )(cons mat 'gen-mat )) body))
 299   `(defun ,name (,mat ,v &aux  info gen-mat gen-vector image)
 300      (remprop ',name 'representation)
 301      (cond ((setq info (get ',name 'representation))
 302             (sub-list (list (first info) ,mat
 303                             (second info) ,v)
 304                       (third info)))
 305            (t    (setq gen-mat ($general_matrix (1- (length (second ,mat))) '$p))
 306                  (setq gen-vector (subseq $aaaa 0 (length ,v)))
 307                  (setq image (progn ,@body))
 308                  (putprop  '$tensor4 (list gen-mat gen-vector image) 'representation)
 309                  ($tensor4 mat v)))))
 310
 311 ;(defrep sta (ma va)
 312 ; (ncmul* ma va))
 313
 314 (defun $tensor4 (mat v &aux info gen-mat gen-vector image)
 315   (cond ((setq info (get '$tensor4 'representation))
 316
 317          (sub-list (list (first info) mat
 318                          (second info) v)
 319                    (third info)))
 320         (t       (setq gen-mat ($general_matrix (1- (length (second mat))) '$p))
 321                  (setq gen-vector (subseq $aaaa 0 (length v)))
 322
 323                  (setq image ($tensor_product_representation
 324                                '$tensor2
 325                                4
 326                                '$tensor2
 327                                4
 328                                gen-mat gen-vector))
 329                  (putprop  '$tensor4 (list gen-mat gen-vector image) 'representation)
 330                  ($tensor4 mat v))))
 331
 332
 333 (defun $tensor4 (mat v)($tensor_product_representation
 334                                '$tensor2
 335                                4
 336                                '$tensor2
 337                                4
 338                          mat v))
 339
 340 (defun $tensor2 (mat v)($tensor_product_representation
 341                          '$standard_rep
 342                          2
 343                          '$standard_rep
 344                          2
 345                          mat v))
 346
 347 (defun $tensor3 (mat v)($tensor_product_representation
 348
 349                          '$tensor2
 350                          4
 351                          '$standard_rep
 352                          2
 353                          mat v))
 354
 355
 356
 357 (defun $tensor2 (mat v)($tensor_product_representation
 358                          '$standard_rep
 359                          3
 360                          '$standard_rep
 361                          3
 362                          mat v))
 363
 364 (defun $tensor3 (mat v)($tensor_product_representation
 365
 366                          '$tensor2
 367                          9
 368                         '$standard_rep
 369                          3
 370
 371                          mat v))
 372
 373 (defun $tensor3 (mat v)($tensor_product_representation
 374                          '$standard_rep
 375                          3
 376                          '$tensor2
 377                          9
 378                          mat v))
 379
 380
 381
 382
 383 (defun $sum5_standard_rep(mat v)
 384   (let ((size ($length mat)))
 385     (loop for i below 5
 386           appending (cdr ($standard_rep mat ($firstn size v))) into tem
 387           do (setq v (cons '(mlist) (nthcdr size (cdr v))))
 388           finally (return (cons '(mlist) tem)))))
 389
 390 ;;The representation of ei^Vej^Vek as a list is done by varying the first
 391 ;;index fastest thus if dim V = 2 then get
 392 ;;[e1^Ve1^Ve1,
 393 ;  e2^Ve1^Ve1,
 394 ;  e1^Ve2^Ve1,
 395 ;  e2^Ve2^Ve1,
 396 ;  e1^Ve1^Ve2,
 397 ;  e2^Ve1^Ve2,
 398 ;  e1^Ve2^Ve2,
 399 ;  e2^Ve2^Ve2]
 400
 401
 402 ;;the first two are in the above representation and the second two are for reversed
 403 ;;indexing where the subscripts vary the fastest at the right.(as usual in zeta lisp)
 404
 405 (defun list-ind-to-tensor-ind (dimv tensor-power ind)
 406   (setq ind (1- ind))
 407   (loop for i below tensor-power
 408         collecting (1+ (mod ind dimv) )
 409         do (setq ind (quotient ind dimv))))
 410
 411
 412 (defun tensor-ind-to-list-ind (dimv &rest ll)
 413   (loop for m in ll
 414         for i from 0
 415         summing (* (expt  dimv i) (1- m)) into tem
 416         finally (return (1+ tem))))
 417
 418 (defun tensor-ind-to-list-ind (dimv &rest ll)
 419   (loop for m in ll
 420         for i to 0 downfrom (1- (length ll))
 421         summing (* (expt  dimv i) (1- m)) into tem
 422         finally (return (1+ tem))))
 423
 424 (defun list-ind-to-tensor-ind (dimv tensor-power ind)
 425   (setq ind (1- ind))
 426   (loop for i below tensor-power
 427         collecting (1+ (mod ind dimv) ) into tem
 428         do (setq ind (quotient ind dimv))
 429         finally (return (nreverse tem))))
 430
 431
 432 (defun $tensor_product_representation (rep1 d1 rep2 d2 p w &aux mm1 mm2 )
 433   (setq mm1 ($find_matrix_of_operator
 434               `(lambda (v) (mfuncall ',rep1 ',p v))
 435                                       ($standard_basis d1)))
 436   (setq mm2 ($find_matrix_of_operator
 437               `(lambda (v) (mfuncall ',rep2 ',p v))($standard_basis d2)))
 438   ($list_matrix_entries (ncmul* (mfuncall '$kronecker_product mm1 mm2) w)))
 439 (defvar $binary_quartic_monomials
 440         ($ratsimp (subst 'mtimes 'mnctimes ($commutative_dot_monomials '((mlist) $x  $y) 4))))
 441
 442 (defun $find_matrix_of_operator (operator basis )
 443   "operator acts on the space spanned by basis"
 444   (let (answer zero-b unknowns gen-a gen-sum gen-b)
 445     (setq gen-b ($firstn ($length basis) $bbbb))
 446     (setq gen-sum ($general_sum basis $bbbb))
 447     (setq answer ($fast_linsolve
 448                   ($list_equations
 449                    (list '(mlist)
 450                          (list '(mequal)
 451                                (mfuncall operator gen-sum)
 452                                ($general_sum basis (setq unknowns
 453                                                          (subseq $aaaa 0 (length basis)))))))
 454                   unknowns))
 455     (setq gen-a ($sublis answer unknowns))
 456     (setq zero-b (loop for w in (cdr gen-b)
 457                     collecting (cons w 0)))
 458
 459     (loop for v in (cdr gen-b)
 460        collecting (meval* (sublis zero-b (subst 1 v  gen-a)))
 461        into tem
 462        finally (return
 463                  ($transpose  (cons '($matrix simp) tem))))))
 464
 465 (defun $find_matrix_of_operator (operator basis &aux answer zero-b
 466                                  unknowns gen-a gen-sum gen-b rhs lhs eqns)
 467   "operator acts on the space spanned by basis"
 468   (setq gen-b ($firstn ($length basis) $bbbb))
 469   (setq gen-sum ($general_sum basis $bbbb))
 470   (setq lhs (mfuncall operator gen-sum))
 471   (setq rhs    ($general_sum basis (setq unknowns (subseq $aaaa 0 (length basis)))))
 472
 473   (setq eqns ($list_equations (list '(mlist) (list '(mequal) lhs rhs))))
 474   (setq answer ($fast_linsolve eqns unknowns))
 475   (setq gen-a ($sublis answer unknowns))
 476   (setq zero-b (loop for w in (cdr gen-b) collecting (cons w 0)))
 477   (loop for v in (cdr gen-b)
 478         collecting (meval* (sublis zero-b (subst 1 v  gen-a)))
 479         into tem
 480         finally (return
 481                   ($transpose  (cons '($matrix simp) tem)))))
 482
 483 (defun $restriction_representation (repr elt-gln vector sub-basis)
 484   (let ((repr-of-p ($find_matrix_of_operator `(lambda (x) (mfuncall ',repr ',elt-gln x))
 485                                              sub-basis)))
 486     (ncmul* repr-of-p vector)))
 487 (defun $standard_basis (n)
 488   (cons '(mlist) (cdr ($ident n))))
 489
 490 (defun $find_highest_weight_vectors (rep n1 n2 &aux uni-vectors answer)
 491
 492                  (setq uni-vectors     ($unipotent_invariant_vectors rep n1 n2))
 493         (setq answer       ($diagonalize_torus rep n1 n2 uni-vectors))
 494         (list  '(mlist) uni-vectors answer))
 495
 496 (defvar *some-primes* '(2 3 5 7 11 13 17 19 23 29 31 37))
 497
 498 (defun $diagonalize_torus (representation n1 n2 basis &aux diag tem mat cpoly answer eigen-values weight-vector)
 499   "here operator has args a and b where a is in torus and b is a sum of elements of basis"
 500   (setq diag ($ident n1))
 501   (loop for i from 1 to n1
 502      do
 503        (setq diag ($setelmx (nth (1- i) *some-primes*) i i diag)))
 504
 505   (setq tem `(lambda (x)
 506                (mfuncall ',representation ',diag x)))
 507   (setq mat  ($find_matrix_of_operator tem basis))
 508   (setq cpoly ($charpoly mat '$x))
 509   (setq answer ($solve cpoly '$x) )
 510   (setq eigen-values (loop for u in (cdr answer)
 511                         collecting (third u)))
 512   (loop for u in eigen-values
 513      collecting ($find_eigenvectors mat u) into part-basis
 514      collecting (setq weight-vector (exponent-vector u n1)) into weight-vectors
 515      do (show weight-vector)
 516      finally (return (list '(mlist) (cons '(mlist) part-basis)
 517                            (cons '(mlist) weight-vectors)))))
 518 (defun $find_eigenvectors (mat eigenvalue &aux answer eqns)
 519   (let* ((row-length ($length (second mat)))
 520         (gen-vector ($firstn row-length $aaaa)))
 521     (setq eqns ($list_equations
 522                  (list '(mlist)
 523                  (list '(mequal) (ncmul* (sub* mat
 524                                                (mul* eigenvalue
 525                                                      ($ident row-length)))
 526                                          gen-vector)
 527                        0))))
 528     (setq answer ($fast_linsolve eqns gen-vector))
 529     ($separate_parameters ($sublis answer gen-vector))))
 530
 531 (defun $rank_of_representation (rep n1 list-of-generators &aux mat)
 532   (setq mat ($general_matrix n1 $cccc))
 533   ($find_rank_of_list_of_polynomials
 534     (loop for v in (cdr list-of-generators)
 535           appending (cdr (mfuncall rep mat v))
 536           into tem
 537           finally (return (cons '(mlist) tem)))))
 538
 539 (defun $general_matrix (n a-list)
 540   (cond ((listp a-list)
 541   ($unlist_matrix_entries a-list n))
 542         (t (loop for i from 1 to n
 543                  appending
 544                  (loop for j from 1 to n
 545                        collecting (new-concat a-list i j))
 546                  into tem
 547                  finally (return ($unlist_matrix_entries (cons '(mlist ) tem) n ))))))
 548
 549 (defun exponent-vector (n &optional (length-vector 10) &aux wnum denom num wdnom)
 550   (cond ((numberp n)
 551   (loop for u in *some-primes*
 552         for j below length-vector
 553         collecting
 554         (loop for i from 1
 555               when (not (zerop (mod n (expt  u i))))
 556               do (return (1- i)))))
 557         ((equal (caar n) 'rat)(setq num (second n) denom (third n))
 558          (setq wnum (exponent-vector num length-vector))
 559          (setq wdnom (exponent-vector denom length-vector))
 560          (loop for v in wnum
 561                for w in wdnom
 562                collecting (- v w)))))
 563
 564
 565 (defun replace-parameters-by-aa (expr &aux unknowns  parameters tem1)
 566           (setq unknowns ($list_variables expr "aa" "par"))
 567           (setq parameters (loop for vv in (cdr unknowns)
 568                                  when (search "par" (string vv) :test #'char-equal)
 569                                  collecting vv))
 570           (setq tem1 (cdr $aaaa))
 571
 572       (loop for vv in parameters
 573             appending (loop while tem1
 574                             when (not (member (car tem1) unknowns :test #'eq))
 575                             collecting (cons vv (car tem1)) into subs1
 576                             and do
 577                             (let ((*print-level* 3))
 578                               (format t "~%Replacing ~A by ~A in ~A." vv (car tem1) expr))
 579                             (setq tem1 (cdr tem1)) (return subs1)
 580                             else
 581                             do (setq tem1 (cdr tem1)))
 582
 583             into subs
 584             finally
 585             (show subs)
 586             (setq expr (sublis subs expr))
 587             (setq unknowns (sublis subs unknowns)))
 588       (values expr unknowns))
 589
 590 (defun function-denominator (pol)
 591   (cond ((numberp pol) (denominator pol))
 592         ((affine-polynomialp pol) 1)
 593         ((rational-functionp pol) (denom pol))
 594         (t (denom (new-rat pol)))))
 595
 596 (defun function-numerator (pol)
 597   (cond ((rationalp pol) (numerator pol))
 598         ((affine-polynomialp pol) pol)
 599         ((rational-functionp pol) (car pol))
 600         (t (car (new-rat pol)))))
 601
 602 (defun poly-type (x )
 603   (cond ((numberp x) ':number)
 604         ((affine-polynomialp x) ':polynomial)
 605         ((rational-functionp x) ':rational-function)
 606         (($ratp x) ':$rat)
 607         ((and (listp x) (eq (caar x) 'rat)) ':rat)
 608         (t nil)))
 609
 610 (defun check-rat-variables(expr)
 611   (let ((varlist (varlist expr))
 612         (genvar (genvar expr)))
 613     (check-vgp-correspond)))
 614
 615 (defun check-vgp-correspond()
 616   (loop for v in varlist
 617         for w in genvar
 618         when (or (not (equal v (get w 'disrep)))
 619                  (not (eq w (get-genvar v))))
 620        do (merror "bad ~A and ~A" v w)))
 621
 622
 623 (defun $slow_gcdlist (&rest ll)
 624   (cond ((null ll) 1)
 625         (($listp (car ll)) (apply '$gcdlist (cdar ll)))
 626         (t
 627   (case (length ll)
 628     (1 (car ll))
 629     (2 ($gcd (car ll) (second ll)))
 630     (otherwise ($gcd (car ll) (apply '$gcdlist (cdr ll))))))))
 631
 632 (defun $slow_projective (l)
 633   ($ratsimp (simplifya (list '(mquotient) l ($slow_gcdlist l)) nil)))
 634
 635 ;
 636 ;(defmacro allow-rest-argument (new-funct binary-op  answer-for-null-argument rest-arg)
 637 ; `(cond  ((null ,rest-arg) ,answer-for-null-argument)
 638 ;        (t     (case (length ,rest-arg)
 639 ;         (,2 (,binary-op (first ,rest-arg)(second ,rest-arg)))
 640 ;         (,1 (car ,rest-arg))
 641 ;         (otherwise (apply ',new-funct (,binary-op (first ,rest-arg)
 642 ;                                          (second ,rest-arg)) (cddr ,rest-arg)))))))
 643
 644
 645
 646 ;(defmacro polyop (x y identity-x identity-y number-op poly-op rat-op &optional rat-switch )
 647 ;  (cond (rat-switch (setq rat-switch '(t))))
 648 ;  `(cond
 649 ;     ((and (numberp ,x)(numberp ,y))(,number-op ,x ,y))
 650 ;     ((eq ,x ,identity-x) ,y)
 651 ;     ((eq ,y ,identity-y) ,x)
 652 ;     (t
 653 ;      (let ((xx (poly-type ,x)) answer
 654 ;           (yy (poly-type ,y)))
 655 ;
 656 ;;      (cond ((or (eq xx '$rat)(eq yy '$rat))(with-vgp(check-vgp-correspond))));;remove later
 657 ;       (cond ((null xx)(setq ,x (cdr ($new_rat ,x)) xx ':rational-function))
 658 ;             ((eq xx ':$rat)
 659 ;              (setq ,x (cdr ,x) xx ':rational-function))
 660 ;             ((eq xx ':rat ) (setq ,x (cons (second ,x) (third ,x))
 661 ;                                  xx ':rational-function)))
 662 ;
 663 ;       (cond
 664 ;         ((null yy)(setq ,y (cdr ($new_rat ,y)) yy ':rational-function))
 665 ;         ((eq yy ':$rat)(setq ,y (cdr ,y) yy ':rational-function))
 666 ;         ((eq yy ':rat ) (setq ,y (cons (second ,y) (third ,y))
 667 ;                              yy ':rational-function)))
 668 ;       (cond ((and (eq xx ':rational-function) (eq (denom ,x) 1))
 669 ;              (setq ,x (car ,x) xx :polynomial)))
 670 ;       (cond ((and (eq yy ':rational-function) (eq (denom ,y) 1))
 671 ;              (setq ,y (car ,y) yy :polynomial)))
 672 ;       (setq answer
 673 ;             (case xx
 674 ;               (:number   (case yy
 675 ;                            (:number (,number-op ,x ,y))
 676 ;                            (:polynomial (,poly-op ,x ,y))
 677 ;                            (:rational-function
 678 ;                             (,rat-op (cons ,x 1) ,y ,@ rat-switch))))
 679 ;               (:polynomial
 680 ;                (case yy
 681 ;                  (:number (,poly-op ,x ,y))
 682 ;                  (:polynomial (,poly-op ,x ,y))
 683 ;                  (:rational-function (,rat-op (cons ,x 1) ,y,@ rat-switch))
 684 ;                  (otherwise (merror "unknown type for polyop "))))
 685 ;               (:rational-function
 686 ;                (case yy
 687 ;                  (:number (,rat-op ,x  (cons ,y 1) ,@ rat-switch))
 688 ;                  (:polynomial (,rat-op ,x (cons ,y 1)  ,@ rat-switch))
 689 ;                  (:rational-function (,rat-op ,x ,y ,@ rat-switch))))
 690 ;               (otherwise (merror "unknown arg"))))
 691 ;       (cond ((affine-polynomialp answer) answer)
 692 ;             ((rational-functionp answer)
 693 ;              (cond ((eq 1 (cdr answer))(car answer))
 694 ;                    (t answer)))
 695 ;             (t answer))))))
 696 ;
 697 ;(defun n* (&rest l)
 698 ;  (case (length l)
 699 ;    (0 1)
 700 ;    (1 (car l))
 701 ;    (2 (n1* (car l) (second l)))
 702 ;    (t (n1* (car l) (apply 'n* (cdr l))))))
 703 ;
 704 ;(defun n+ (x y) (polyop x y  0 0 + pplus ratplus ))
 705 ;(defun n1* (x y)(polyop x y 1 1 * ptimes rattimes t))
 706 ;(defun n- (x y)(polyop x y nil 0 - pdifference ratdifference))
 707 ;(defun nred (x y &aux answer)
 708 ;  (setq answer (polyop x y nil 1  ratreduce ratreduce ratquotient))
 709 ;  (setq answer (rationalize-denom-zeta3 answer))
 710 ;  (cond ((numberp answer) answer)
 711 ;       ((eq (denom answer) 1) (car answer))
 712 ;       (t answer)))
 713
 714
 715
 716 ;(defun new-disrep (expr)
 717 ;  (cond ((atom expr) expr)
 718 ;       (t
 719 ;        (let ((type (poly-type expr)))
 720 ;          (case type
 721 ;            (:number expr)
 722 ;            (:polynomial ($ratdisrep (header-poly expr)))
 723 ;            (:rational-function
 724 ;             (cond (($numberp expr)
 725 ;                    (cond ((zerop (num expr) 0))
 726 ;                          (t
 727 ;                    (list '(rat simp) (num expr) (denom expr)))))
 728 ;                   (t
 729 ;                    ($ratdisrep (header-poly expr)))))
 730 ;            (otherwise  (cond ((mbagp expr)(cons (car expr) (mapcar 'new-disrep expr)))
 731 ;                              (($ratp expr)($ratdisrep expr))
 732 ;                              (t expr))))))))
 733
 734
 735
 736 (defun sp-add* (&rest llist)
 737   (let ((varlist))
 738     (apply 'add*  llist)))
 739
 740 (defun sp-minus* (a)
 741   (let ((varlist))
 742     (simplifya (list '(mminus) a) nil)))
 743
 744 (defun sp-sub* (a b )
 745   (sp-add* a (sp-minus* b)))
 746
 747
 748 (defun sp-mul* (&rest a-list)
 749   (let ((varlist ))
 750     (cond ((< (length a-list) 3)(apply 'mul* a-list))
 751           (t (sp-mul* (car a-list) (apply 'sp-mul* (cdr a-list)))))))
 752
 753 (defun sp-div* (a b)
 754   (let ((varlist))
 755    (simplifya  (list '(mquotient) a b) nil)))
 756
 757
 758 ;;the following don't always seem to work
 759 (defun sp-add* (&rest l)(allow-rest-argument sp-add* n+ 0 l))
 760 (defun sp-mul* (&rest l)(allow-rest-argument sp-mul* n* 1  l))
 761 (defun sp-div* (a b) (nred a b))
 762 (defun sp-minus* ( b)(n- 0 b))
 763 (defun sp-sub* (a b)(n- a b))
 764
 765
 766 (defun $numberp (n)
 767   "we changed this to allow polynomial-type"
 768   (cond ((numberp n) t)
 769         ((atom  n) nil)
 770         ((and (numberp (car n))(numberp (cdr n))) t)
 771         ((atom (car n)) nil)
 772         ((eq (caar n) 'rat))
 773         ((eq (caar n) 'mrat)(and (numberp (cadr n)) (numberp (cddr n))))
 774         (($floatnump n))
 775         (($bfloatp n))))
 776
 777 (defun poly-scalarp (a &aux tem)
 778   (cond ((numberp a))
 779         ((atom a)($scalarp a))
 780         ((setq tem(affine-polynomialp a)) ($scalarp tem))
 781         ((and (setq tem(affine-polynomialp (car a))) (affine-polynomialp (cdr a)))
 782          ($scalarp tem))))
 783
 784
 785 (defun $coefficient_matrix (a-list variables)
 786   (check-arg a-list  '$listp nil)
 787   (check-arg variables  '$listp nil)
 788   (loop for w in (cdr a-list)
 789         collecting
 790         (loop for v in (cdr variables)
 791               collecting
 792               ($nc_coeff w v) into tem
 793               finally (return (cons '(mlist) tem)))
 794         into bil
 795         finally (return (cons '($matrix) bil))))
 796 ;
 797 ;
 798 ;
 799 ;
 800 ;
 801 ;
 802 ;
 803 ;
 804 ;
 805 ;
 806 ;(DEFMFUN new-RATREP* (X)
 807 ;  (LET ()                                      ; (GENPAIRS)
 808 ;;          (ORDERPOINTER VARLIST)  ;;creates enough new genvar and numbers them done
 809 ;    (RATSETUP1 VARLIST GENVAR)
 810 ;    ;which is now in new-ratf
 811 ;
 812 ;    (MAPC #'(LAMBDA (Y Z) (cond ((assolike y genpairs)  nil)
 813 ;                               (t (PUSH (CONS Y (RGET Z)) GENPAIRS))))
 814 ;         VARLIST GENVAR)
 815 ;
 816 ;    (RATSETUP2 VARLIST GENVAR)
 817 ;    (XCONS (PREP1 X)                           ; PREP1 changes VARLIST
 818 ;          (LIST* 'MRAT 'SIMP VARLIST GENVAR    ;    when $RATFAC is T.
 819 ;                 (IF (AND (NOT (ATOM X)) (MEMber 'IRREDUCIBLE (CDAR X) :test #'eq))
 820 ;                     '(IRREDUCIBLE))))))
 821 ;
 822 ;(DEFMFUN new-RATREP* (X)
 823 ;  (LET ()                                      ; (GENPAIRS)
 824 ;;          (ORDERPOINTER VARLIST)  ;;creates enough new genvar and numbers them done
 825 ;    (RATSETUP1 VARLIST GENVAR)
 826 ;    ;which is now in new-ratf
 827 ;;
 828 ;;    (MAPC #'(LAMBDA (Y Z) (cond ((assolike y genpairs)  nil)           ;;do this in new
 829 ;;                              (t (PUSH (CONS Y (RGET Z)) GENPAIRS))))
 830 ;;        VARLIST GENVAR)
 831 ;
 832 ;    (RATSETUP2 VARLIST GENVAR)
 833 ;    (XCONS (PREP1 X)                           ; PREP1 changes VARLIST
 834 ;          (LIST* 'MRAT 'SIMP VARLIST GENVAR    ;    when $RATFAC is T.
 835 ;                 (IF (AND (NOT (ATOM X)) (MEMber 'IRREDUCIBLE (CDAR X) :test #'eq))
 836 ;                     '(IRREDUCIBLE))))))
 837 ;
 838 ;(Defun New-RATF (L &aux tem)
 839 ;  (with-vgp
 840 ;    (PROG (U *WITHINRATF*)
 841 ;         (SETQ *WITHINRATF* T)
 842 ;         (WHEN (EQ '%% (CATCH 'RATF (new-NEWVAR L)))
 843 ;           (SETQ *WITHINRATF* NIL) (RETURN (SRF L)))
 844 ;         (loop for v in varlist
 845 ;               for i from 1
 846 ;               when  (setq tem (get-genvar v))
 847 ;               collecting tem into a-list
 848 ;               else
 849 ;               collecting
 850 ;               (prog1 (setq tem (gensym))
 851 ;                      (putprop tem v 'disrep)
 852 ;                      (show (list tem v 'disrep))
 853 ;                      (push (cons v (rget tem)) genpairs )
 854 ;                      )
 855 ;               into a-list
 856 ;               do (setf (symbol-value tem) i)
 857 ;               finally (setq genvar a-list))
 858 ;         (SETQ U (CATCH 'RATF (new-RATREP* L)))        ; for truncation routines
 859 ;         (RETURN (OR U (PROG2 (SETQ *WITHINRATF* NIL) (SRF L)))))))
 860 ;
 861 ;(DEFUN new-NEWVAR (L  )
 862 ;  (let (( vlist varlist))
 863 ;
 864 ;  (NEWVAR1 L)
 865 ;  (setq varlist (sortgreat vlist))
 866 ;  vlist))
 867  ; (SETQ VARLIST (NCONC (SORTGREAT VLIST) VARLIST)))
 868
 869
 870 (defun $new_rat (expr)
 871   (cond
 872     ((and ($ratp expr)(equal *genvar* (genvar expr))(equal *varlist* (varlist expr)))
 873      expr)
 874     (($ratp expr)(cond ((loop for v in (varlist expr)   ;;what about minimize varlist here.
 875                               for w in (genvar expr)    ;or in new-ratf
 876                               when
 877                               (not (eq w (get-genvar v)))
 878                               do
 879                               (return 'bad-order))
 880                         (new-ratf expr))
 881                        (t expr)))
 882     ((mbagp expr) (cons (car expr) (mapcar #'$new_rat (cdr expr))) ($new_rat expr))
 883     ((and (not (atom expr))(not (atom (car expr))))(cond ((equal (caar expr) 'rat) expr)
 884                                                          (t (new-ratf expr))))
 885
 886     (t (new-ratf expr))))
 887 ;(if (mbagp exp) (cons (car exp) (mapcar #'rat0 (cdr exp))) (ratf exp)))
 888
 889
 890 (defun $shift (expr &optional (variables $current_variables) &aux tem)
 891  (sublis (loop for v on (cdr variables)
 892         when (second v)
 893         do (setq tem  (second v))
 894         else do (setq tem (second variables))
 895
 896         collecting (cons (car v) tem)) expr))
 897 (defun $spur (expr function n &aux (tem expr) (answer 0))
 898   (loop for i below n
 899         do
 900         (setq answer (add* answer (setq tem (funcall function tem)))))
 901  ($ratsimp answer))
 902
 903
 904 (defun $basis_of_invariant_ring (n &optional (variables $current_variables) &aux
 905                                  monoms trace_monoms f new-f eqns solns tem ar sp pivots)
 906      (setq monoms ($mono variables n))
 907         (setq trace_monoms (loop for u in (cdr monoms)
 908                             collecting ($spur u '$shift (1- (length variables)) ) into tem
 909                             finally (return (cons '(mlist ) tem))))
 910         (break t)
 911         (setq f ($general_sum trace_monoms $aaaa))
 912         (setq new-f ($dotsimp f))
 913         (setq eqns ($extract_linear_equations (list '(mlist) new-f) monoms))
 914         (setq solns ($fast_linsolve eqns (subseq $aaaa 0 (length monoms))))
 915         (setq sp(send $poly_vector :the-sparse-matrix))
 916     (setq ar (send sp ':column-used-in-row))
 917    (setq pivots (loop for i below (length (the cl:array ar))
 918           when (setq tem (aref ar i))
 919           collecting tem into tem1
 920           finally (return  tem1)))
 921    (loop for i in pivots
 922                  collecting (nth (1+ i) trace_monoms)
 923                  into tem1
 924                  finally (return (cons '(mlist) tem1))))
 925
 926 (defvar $mmm)
 927 (defvar $conditions_on_q nil)
 928
 929 (defun $case_rep (mat vect triple  &aux  answer  subs ($expop 100)
 930                  ind-used vmmm prod mat.rtx relat-mat)
 931   (setq triple (cdr triple))
 932   (setq subs (list (cons '$%alpha (first triple))
 933                    (cons '$%beta (second triple))
 934                    (cons '$%gamma (third triple))))
 935
 936   (let ((mmm (sublis subs $mmm)))
 937     (setq mmm ($ratsimp mmm))
 938     (displa mmm)
 939     (loop for  v in (cdr $conditions_on_q)
 940           for i from 1
 941           when ($zerop ($ratsimp (sublis subs v)))
 942           collecting (nth i mmm) into tem
 943           and
 944           collecting i into ind
 945           finally (setq ind-used ind)(setq relat-mat tem))
 946     (format t "~%Dimension of representation is ~A." (length ind-used))
 947     (show ind-used)
 948     (setq vmmm
 949           (loop for i from 1 to (length ind-used)
 950                 collecting (mul* (nth i vect)(nth (1- i) relat-mat)) into tem
 951                 finally (return (meval* (cons '(mplus) tem)))))
 952     (setq mat.rtx (ncmul* mat '(($matrix simp)
 953                                 ((mlist simp) $x) ((mlist simp) $y) ((mlist simp) $z))))
 954
 955     (setq prod
 956           (ncmul*  ($transpose mat)
 957                    ($sub_list `((mlist simp) ((mequal simp) (($matrix simp)
 958                                                              ((mlist simp) $x)
 959                                                              ((mlist simp) $y)
 960                                                              ((mlist simp) $z))
 961                                               ,mat.rtx)) vmmm) mat))
 962     (setq answer (mfuncall '$cof prod))
 963     (loop for i in ind-used
 964           collecting (nth i answer) into tem
 965           finally (return (cons '(mlist) tem)))))
 966
 967 (defun $make_commutative (expr)
 968   (resimplify (subst 'mtimes 'mnctimes expr)))