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 ((*nopoint t) *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))) 969