share/affine/dim-3.lisp

   1 ;;; -*- Mode:Lisp; Package:CL-MAXIMA; Syntax:COMMON-LISP; Base:10 -*-
   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 point-of-tensors (special-tensor standard-tensor
  11                                &key (nc-monomials ($list_nc_monomials standard-tensor)))
  12  (let ((  st-tensor (st-rat standard-tensor))
  13        ( sp-tensor (st-rat special-tensor))
  14        (nc-monoms  (st-rat nc-monomials)))
  15   (setq nc-monoms  (sort nc-monoms 'alphalessp))
  16    (loop for v in nc-monoms
  17          collecting (pcoeff st-tensor v) into variable-names
  18          collecting (pcoeff sp-tensor v) into variable-values
  19          finally (return
  20                    (loop for vv in variable-names
  21                          for ww in variable-values
  22                          when (and (consp vv) (affine-polynomialp vv)
  23                                    (eql (length vv) 3))
  24                          collecting (cons (car vv) ww))))))
  25
  26 (defun gen-matrix-rows (mat )
  27   (cond ((matrix-p mat)(cdr mat))
  28         (($matrixp mat)  (mapcar 'cdr (cdr mat)))
  29         (t (merror "not a matrix"))))
  30
  31 (defun point-of-algebras (special-tensor standard-tensor q standard-q &aux tensor-subs)
  32   (setq tensor-subs (point-of-tensors special-tensor standard-tensor))
  33   (append tensor-subs
  34           (loop for v in (mapcar #'st-rat (apply #'append (gen-matrix-rows q)))
  35              for w in (mapcar #'st-rat (apply #'append (gen-matrix-rows standard-q)))
  36         collecting (cons (car w) v))))
  37
  38 (defun $dimensions_for_replacements ($repl-list &optional (degree 7))
  39  (let (($dot_simplifications
  40          (cons '(mlist) (loop for v in (cdr $repl-list)
  41                               collecting v collecting 0)))
  42        ($current_variables ($list_variables $repl-list)))
  43    (loop for i below degree do
  44         (format t "~%The dimension in deg ~A is ~A and hilbert is ~A"
  45                 ($length ($mono $current_variables i))
  46                 i
  47                 (hilbert i)))))
  48
  49 (defun $normalizing_element_p (elmt &aux tem answ eqns solns  oth)
  50   (loop for v in (cdr $current_variables)
  51         do (setq tem (sub*  (ncmul* v elmt)
  52                             (ncmul*  elmt (setq oth ($general_sum $current_variables $aaaa )))))
  53         (setq answ ($totaldisrep ($numerator ($dotsimp tem))))
  54         (setq eqns ($extract_linear_equations  (list '(mlist) answ)))
  55         (setq solns ($fast_linsolve eqns (subseq $aaaa 0 (length $current_variables))))
  56         collecting ($sublis solns oth) into final
  57         do (mshow (cons '(mlist)final))
  58         finally (return (cons '(mlist) final))))
  59
  60 (defmacro while (test &body body)
  61        `(loop
  62            (when (null ,test) (return))
  63            ,@body))
  64
  65 (defun collect-atoms (tree &optional (pred #'(lambda (x) t)) &aux result)
  66    (declare (special result))
  67    (collect-atoms1 tree pred) result)
  68
  69 (defun collect-atoms1 (tree pred)
  70    (declare (special result))
  71   (cond ((atom (car tree))(cond ((funcall pred (car tree))
  72                                  (pushnew (car tree) result :test 'equal))))
  73         (t (collect-atoms1 (car tree) pred)))
  74   (cond ((atom (cdr tree))(cond ((funcall pred (cdr tree))
  75                                  (pushnew (cdr tree) result :test 'equal))))
  76         (t (collect-atoms1 (cdr tree) pred))))
  77
  78
  79 (defun $find_coefficients_for_linear_combination( items in-terms-of monoms &optional (cofs $aaaa)
  80                  &aux (gen-sum 0) rest-monoms final-answ answ tem dif)
  81    (setq monoms (nreverse (st-rat  monoms)))
  82   (while (setq in-terms-of (cdr in-terms-of))
  83     (setq cofs (cdr cofs))
  84     (setq gen-sum (n+ gen-sum (n* (car cofs) (car in-terms-of)))))
  85   (while (setq items (cdr items))
  86     (setq dif (n- (car items) gen-sum))
  87     (setq answ nil)
  88     (setq rest-monoms  monoms)
  89     (while (and (setq tem (car rest-monoms))(setq rest-monoms (cdr rest-monoms)))
  90               (setq tem (gencoeff dif tem))
  91              (push (new-disrep tem) answ))
  92     (push (cons '(mlist) answ) final-answ))
  93   (cons '(mlist) final-answ))
  94
  95 (defun $generate_matrix (fun i j)
  96   (cons '($matrix)
  97         (loop for ii from 1 to i
  98            collecting (cons'(mlist)
  99                            (loop for jj from 1 to j
 100                                   collecting (meval* `((,fun) ,ii ,jj)))))))
 101
 102
 103 (defun gen-matrix-row-plus (row1 row2 &key row1-multiple  row2-multiple)
 104        "form row1+row2 where the multiple you may specify multiples."
 105        (cond (row1-multiple (setq  row1-multiple (st-rat  row1-multiple))))
 106        (cond (row1-multiple (setq  row2-multiple (st-rat  row2-multiple))))
 107        (loop for v in (cdr row1) for w in (cdr row2)
 108              when row1-multiple do (setq v (n* row1-multiple v))
 109              when row2-multiple do (setq w (n* row2-multiple w))
 110              collecting (n+ v w) into result
 111              finally (return (cons '(mlist) result))))
 112
 113
 114 (defun $interchange_matrix_columns (matrix i j &optional copy)
 115   (cond (copy (setq matrix
 116                     (cons (car matrix)
 117                           (loop for u in (cdr matrix)
 118                                 with max =  (1+ (max i j))
 119                                 collecting (nconc (subseq u 0 max) (nthcdr max u)))))))
 120
 121   (loop for u in (cdr matrix)
 122         do (swapf (nth i u) (nth j u)))
 123   matrix)
 124
 125 (defun $interchange_matrix_rows (matrix i j &optional copy)
 126   "interchange MATRIX rows I and J.  If COPY then copy then don't do it destructively"
 127   (cond (copy (setq matrix (nconc (subseq matrix 0 (setq copy (1+ (max i j))))
 128                                    (nthcdr copy matrix )))))
 129   (swapf (nth i matrix) (nth j matrix))
 130   matrix)
 131
 132 (defmacro matrix-entry (matrix i j)
 133   `(nth ,j (nth ,i ,matrix)))
 134
 135 (defun $multiply_row (list factor)
 136   (cons (car list) (loop for v in (cdr list) collecting (n* factor v))))
 137
 138 ;;still must remove the factors along the diagonal.
 139 (defun $decompose_symmetric_matrix (matrix &aux (ident ($ident ($length matrix))) fact
 140                                     ident-pivot-row pivot-row orig )
 141   "Write MATRIX as a product of P.DIAGONAL.TRANSPOSE(P).  This returns value P."
 142   (setq orig matrix)
 143   (setq matrix (cons (car matrix) (loop for v in (cdr matrix) collecting (copy-list v))))
 144   (loop for ii from 1 below ($length matrix)
 145         do
 146    (loop for i from ii to ($length matrix)
 147          when (not (pzerop (matrix-entry matrix i ii)))
 148          do
 149          (cond ((eq i ii) (return))
 150                (t (setf (nth ii matrix) (gen-matrix-row-plus (nth i matrix) (nth ii matrix)))
 151                   (setf (nth ii ident) (gen-matrix-row-plus (nth i ident) (nth ii ident)))
 152                   (loop for i1 from ii to ($length matrix)
 153                         do (setf (matrix-entry matrix i1 ii)
 154                                  (n+ (matrix-entry matrix i1 ii) (matrix-entry matrix i1 i))))
 155                   (return)))
 156          finally (fsignal "This implementation assumes the matrix is nonsingular"))
 157    (setq pivot-row ($multiply_row (nth ii matrix) (setq fact (nred 1 (matrix-entry matrix ii ii)))))
 158    (setf ident-pivot-row ($multiply_row (nth ii ident) fact))
 159    (loop for i from (1+ ii) to ($length matrix)
 160          do (setf (nth i matrix) (gen-matrix-row-plus (nth i matrix) pivot-row
 161                                     :row2-multiple (setq fact (n* -1 (matrix-entry matrix i ii)))))
 162          (setf (nth i ident)(gen-matrix-row-plus (nth i ident) ident-pivot-row
 163                                     :row2-multiple fact))))
 164 ;should be diag; (let(( answ ($new_disrep ident)))    (show (ncmul* answ orig ($transpose answ))))
 165  ($new_disrep ident))
 166
 167 (defun $new_disrep (item)
 168   (cond ((atom item) (new-disrep item))
 169         ((and (listp (car item))(mbagp item)) (cons (car item) (mapcar '$new_disrep (cdr item))))
 170         (t (new-disrep item))))
 171
 172 (defun $symbol_matrix(symbol rows columns)
 173   "Generate the matrix of symbols starting with SYMBOL with ROWS rows
 174   and COLUMNS columns."
 175   (cons '($matrix simp)
 176         (loop for i from 1 to rows
 177               collecting (loop for j from 1 to columns
 178                                collecting
 179                                ($concat symbol i j) into tem
 180                                finally (return (cons '(mlist simp) tem))))))
 181
 182 (defun replace-symbols (form containing-string &key prefix replacements &aux subs vari)
 183   (setq vari ($list_variables form containing-string))
 184   (setq subs (cond (prefix (loop for v in (cdr vari)
 185                                  collecting (cons v ($concat prefix v))))
 186                    (t (check-arg replacements (and (listp replacements)
 187                                                    (>= ($length replacements)
 188                                                        ($length vari))) "long enough macsyma-list")
 189                       (cond (prefix (setf (cdr vari) (sort (cdr vari) 'alphalessp))))
 190                       (pairlis (cdr vari) (subseq (cdr replacements) 0 ($length vari))))))
 191   (sublis subs form))
 192
 193 (defun check-cases (roots-of-one &aux relats)
 194   (loop for v in (cdr roots-of-one)
 195         do (setq relats  (subst v '$b $relationsb))
 196       ($set_up_dot_simplifications relats)
 197       (loop for i in '(2 5 10 3 4)
 198             do ($check_overlaps (1+ i) t nil)
 199             (push ($fast_central_elements $current_variables i) $centrals_so_far)
 200             (displa (cons '(mlist) $centrals_so_far)))))
 201
 202
 203 (defun check-centrals (in-degree &key relations set_up &aux)
 204   (cond (set_up     ($set_up_dot_simplifications relations)))
 205       (loop for i in in-degree
 206             do ($check_overlaps (1+ i) t nil)
 207             (push ($fast_central_elements $current_variables i) $centrals_so_far)
 208             (displa (cons '(mlist) $centrals_so_far))))
 209
 210
 211
 212
 213 (defun $check_transpose_condition_2_1 (mm q &optional (variables $current_variables) &aux he)
 214     (setq he (ncmul*  variables mm))
 215     (sub* he ($transpose (ncmul* q mm ($transpose variables)))))
 216
 217 (defun $m_from_f (f &aux n)
 218   (loop for v in (cdr f) when (not (numberp v)) do (setq n ($nc_degree v)) (return))
 219   (let ((monoms ($mono $current_variables (1- n))))
 220     (cons '($matrix)
 221           (loop for fi in (cdr f)
 222                 collecting
 223                 (cons '(mlist)
 224                       (loop for v in (cdr $current_variables)
 225                             collecting
 226                             (loop for w in (cdr monoms)
 227                                   with answ = 0 with mon
 228                                   do (setq answ (add* answ (mul*  ($ratcoef fi (ncmul* w v))
 229                                                                   w)))
 230
 231                                   finally (return answ))))))))
 232
 233 (defun $diagonal (lis)
 234   (cons '($matrix simp) (loop for v in (cdr lis)
 235                               for i from 0
 236                               with tem
 237                               do (setq tem (make-list (length (cdr lis)) :initial-element 0))
 238                               (setf (nth i tem) v)
 239                               collecting (cons '(mlist) tem))))
 240
 241 ;;checks that the relations of specified type, satisfy the xm=qf rule
 242 ;;I did it for all diagonal types and they were ok.
 243 ;(defun $dim3_relations_check(type number-variables &aux rels subs-and-relations mat qq)
 244 ;  (declare (special $relations_2 $relations_3 $rtx $rtx3))
 245 ;  (setq subs-and-relations
 246 ;       (case number-variables
 247 ;         (2 (setq rels $relations_2)
 248 ;            (setq mat $rtx)
 249 ;            (setq vars '((mlist) $%alpha $%beta ))
 250 ;            (case type
 251 ;              ($a (nth 1 rels))
 252 ;              ($e (nth 2 rels))
 253 ;              ($h (nth 3 rels))
 254 ;              ($s1(nth 4 rels))
 255 ;              ($s2(nth 5 rels))
 256 ;              ($s2p(nth 6 rels))
 257 ;              (t (fsignal "type one of a,b,h,s1,s2,or s2p"))))
 258 ;         (3 (setq rels $relations_3)
 259 ;            (setq vars '((mlist) $%alpha $%beta $%gamma))
 260 ;            (setq mat $rtx3)
 261 ;            (case type
 262 ;              ($a (nth 1 rels))
 263 ;              ($b (nth 2 rels))
 264 ;              ($e (nth 3 rels))
 265 ;              ($h(nth 4 rels))
 266 ;              ($s1(nth 5 rels))
 267 ;              ($s1p(nth 6 rels))
 268 ;              ($s1pp(nth 7 rels))
 269 ;              ($s2(nth 8 rels))
 270 ;              (t (fsignal "type one of a,b,e,h,s1,s1p,s1pp,  or s2p"))))))
 271 ; (setq rels  ($totaldisrep  ($sub_list `((mlist) ((mequal)
 272 ;                                                 ,vars ,(second subs-and-relations)))
 273 ;                                      ($list_matrix_entries
 274 ;                                        (ncmul*
 275 ;                                          (third subs-and-relations) mat)))))
 276 ; (setq qq ($diagonal (second subs-and-relations)))
 277 ;  ($ratsimp ($check_transpose_condition_2_1 ($m_from_f rels) qq)))
 278
 279 (defun $dim3_relations(type number-variables &aux rels subs-and-relations mat)
 280   (declare (special $relations_2 $relations_3 $rtx $rtx3))
 281   (setq subs-and-relations
 282         (case number-variables
 283           (2 (setq rels $relations_2)
 284              (setq mat $rtx)
 285              (setq vars '((mlist) $%alpha $%beta ))
 286              (case type
 287                ($a (nth 1 rels))
 288                ($e (nth 2 rels))
 289                ($h (nth 3 rels))
 290                ($s1(nth 4 rels))
 291                ($s2(nth 5 rels))
 292                ($s2p(nth 6 rels))
 293                (t (fsignal "type one of a,b,h,s1,s2,or s2p"))))
 294           (3 (setq rels $relations_3)
 295              (setq vars '((mlist) $%alpha $%beta $%gamma))
 296              (setq mat $rtx3)
 297              (case type
 298                ($a (nth 1 rels))
 299                ($b (nth 2 rels))
 300                ($e (nth 3 rels))
 301                ($h(nth 4 rels))
 302                ($s1(nth 5 rels))
 303                ($s1p(nth 6 rels))
 304                ($s1pp(nth 7 rels))
 305                ($s2(nth 8 rels))
 306                (t (fsignal "type one of a,b,e,h,s1,s1p,s1pp,  or s2p"))))))
 307   ($totaldisrep ($sub_list `((mlist) ((mequal)
 308                                                   ,vars ,(second subs-and-relations)))
 309                                        ($list_matrix_entries
 310                                          (ncmul*
 311                                            (third subs-and-relations) mat)))))
 312
 313 (defun $maybe_ldata_solve (eqns &key (inequality 1) yes &aux ld)
 314  (declare (special $answer $last_equations ))
 315  (setq $last_equations eqns )
 316  (format t "The value of Last_equations is :..")(displa eqns)
 317  (cond ((or yes (y-or-n-p "solve the system"))
 318          (setq ld (make-ldata :eqns (st-rat eqns)))
 319          (setq $answer (simplify-ldata ld :open-g (st-rat inequality)))
 320          (setq $answer (delete-redundant-ldata $answer :ignore-ldata-inequalities t))
 321          (format t "~%**The final answer is (stored in $answer) :")
 322           ($ldata_disrep $answer))
 323         (t eqns)))
 324
 325
 326
 327 (defun $normalizing_conditions (deg-or-element &optional (auto-data) &aux  tem ld eqns )
 328     (declare (special $answer gen $current_variables  $bbbb))
 329   "Find the conditions for an ELEMENT or if element is a number for the
 330   general nc-polynomial in that degree to be normalized.  You may
 331   supply the optional equations for the automorphispm group of the
 332   algebra to assist in the computation.  The current variables and
 333   dotsimps are used."
 334   (cond ((numberp deg-or-element)
 335          (setq gen ($general_sum ($mono $current_variables deg-or-element) $bbbb)))
 336         (t (setq gen deg-or-element)))
 337   (setq tem (cons '(mlist) (loop for v in (cdr $current_variables)
 338                                  collecting (sub* (ncmul* v gen)
 339                                                   (ncmul* gen ($scalar_sum
 340                                                                ($concat '$c (string-grind v))
 341                                                                 $current_variables))))))
 342   (mshow tem)
 343   (setq tem ($dotsimp tem))
 344   (mshow tem)
 345   (setq eqns ($extract_linear_equations ($totaldisrep ($numerator tem))))
 346   (cond (auto-data (loop for v in (cdr auto-data )
 347                          collecting ($maybe_ldata_solve ($append v eqns) :yes t ) into all
 348                          finally (return (apply '$append all))))
 349         (t
 350          ($maybe_ldata_solve eqns))))
 351
 352
 353 (defun $ldata_disrep (lis)
 354   (cond ((ldatap lis)
 355          (cons '(mlist) (loop for v in (ldata-eqns lis)
 356                collecting (new-disrep v))))
 357         ((ldatap (car lis))
 358          (cons '(mlist) (mapcar '$ldata_disrep lis)))
 359         (t (fsignal "ldata_disrep only ldata or lists of ldata"))))
 360
 361 (defun $linear_automorphism_group (&optional
 362                                     (relations ($relations_from_dot_simps))
 363                                     (set_up_simps nil)
 364                                     (variables $current_variables)
 365                                             &aux mat subs subs-for-sublis eqns
 366                                             image deg)
 367   (mshow relations)
 368   (or variables (seta $current_variables ($list_nc_variables relations)))
 369   (mshow $current_variables)
 370   (cond (set_up_simps
 371          (setq deg (apply 'max (mapcar '$nc_degree (cdr relations))))
 372          ($set_up_dot_simplifications relations deg)))
 373   (setq subs (loop for v in (cdr variables)
 374                     collecting ($scalar_sum
 375                                 ($concat '$c (string-grind v))
 376                                 variables)))
 377   (mshow (cons '(mlist)subs))
 378   (setq mat  ($coefmatrix (cons '(mlist) subs) variables))
 379   (setq subs-for-sublis
 380         (cons '(mlist)
 381               (loop for v in (cdr variables)
 382                  for w in subs
 383                  collecting `((mequal) ,v ,w))))
 384   (mshow subs-for-sublis)
 385   (setq image ($dotsimp ($expand ($sublis subs-for-sublis
 386                        relations))))
 387   (setq eqns ($extract_linear_equations ($totaldisrep ($numerator image))))
 388   ($maybe_ldata_solve eqns :inequality ($determinant mat)))
 389
 390 (defremember $dotsimp_remember (elmnt &optional ($dot_simplifications $dot_simplifications))
 391   ($dotsimp elmnt))
 392
 393 (defun $algebra_trace (element  basis &aux cof tem rbasis  den num)
 394   (setq rbasis (st-rat basis))
 395   (setq vars (loop for v in rbasis
 396                    when (listp v) collecting (car v)
 397                   and  do (assert (get (car v) 'disrep))))
 398   (loop for v in (cdr basis)
 399         for rv in rbasis
 400         do (show v rv)
 401         with tra = 0
 402         do (show v) (setq tem  ($dotsimp_remember (ncmul* element v)))
 403         (setq den (function-denominator tem))
 404         (setq num (function-numerator tem))
 405         (setq cof  (pcoeff num rv vars))
 406         when (not (pzerop cof)) do (mshow num cof)
 407         do
 408         (setq tra (n+ tra (nred cof den)))
 409         finally (return (new-disrep tra))))
 410
 411 (defun poly-scalar-p (poly)
 412   (and (affine-polynomialp poly) (or (atom poly) ($scalarp (get (car poly) 'disrep)))))
 413
 414 (defun set-up-trace-subs (basis)
 415     (check-arg basis $listp "maxima list")
 416   (loop for v in  (cdr basis)
 417         collecting (add-newvar v) into old
 418         collecting ($algebra_trace v basis) into new
 419         finally (return (list old new))))
 420
 421 (defun $algebra_trace_matrix (basis &aux elt subs ar result simp )
 422   (setq subs  (set-up-trace-subs basis))
 423   (setq tr (make-pairing-function (cdr basis) (second subs) :test 'nc-equal))
 424   (check-arg basis $listp "maxima list") (setq basis (cdr basis))
 425   (setq ar (make-array (list (length basis) (length  basis))))
 426   (loop for i from 0
 427         for v in basis
 428         do
 429         (loop for j from i below (length  basis)
 430               do (setq elt (ncmul* (nth i basis) (nth j basis)))
 431               (show elt)
 432               (setq simp (new-rat ($dotsimp_remember elt)))
 433               (setq result (new-disrep (apply-linear-function tr  simp)))
 434               (setf (aref ar i j) result)
 435               (setf (aref ar j i) result)))
 436   (values (maxima-matrix-from-array ar) subs tr))
 437
 438 (defun $trace_general_element (elt basis basis-trace &aux result simp)
 439   (setq tr (make-pairing-function (cdr basis) (cdr basis-trace) :test 'nc-equal))
 440   (setq simp (new-rat ($dotsimp_remember elt)))
 441   (setq result (new-disrep (apply-linear-function tr  simp))))
 442
 443
 444 (defun maxima-matrix-from-array  (ar &aux (dims  (array-dimensions ar)))
 445   (cons '($matrix) (loop for i below (car dims)
 446                          collecting
 447                          (cons '(mlist) (loop for j below (second dims)
 448                                              collecting (aref ar i j))))))
 449
 450 (defun make-pairing-function (domain range &key (test 'alike))
 451   #'(lambda (u) (loop for v in domain
 452                       for w in range
 453                       when (funcall test u v)
 454                       do (return w)
 455                       finally (error "u was not in the domain ~A" domain))))
 456
 457
 458 (defun apply-linear-function (fn expr &aux (den 1) result num-answ)
 459   "fn should be a function which is defined on the monomials and this extends it to all"
 460    (cond ((affine-polynomialp expr) (setq den 1))
 461          ((rational-functionp expr)   (setq den (denom expr))(setq expr (num expr)))
 462          (t (fsignal "expr must be a polynomial or rational function")))
 463    (show den)
 464    (setq num-answ
 465          (cond ((atom expr)(nred (n* expr (funcall fn 1)) den))
 466                (t (cond (($scalarp (get (car expr) 'disrep))(show expr)
 467                          (n* expr (funcall fn 1)))
 468                         (t (assert (poly-scalar-p (p-cof expr)))
 469                            (setq result(n* (p-cof expr)
 470                                            (funcall fn (get (car expr) 'disrep))))
 471                            (cond ((cdddr expr)
 472                                   (assert (and (eql (fourth expr) 0)
 473                                                (null (nthcdr 5 expr))))
 474                                   (n+ (apply-linear-function fn (fifth expr))
 475                                       result))
 476                                  (t result)))))))
 477    (nred num-answ den))
 478
 479 (defun dual-basis-element (number-of-element-in-basis trace-matrix &optional (basis $basis)
 480                                                       &aux result cond0 actual-conditions
 481                                                       vari)
 482   (declare( special $basis))
 483   "if trace-matrix is the <ui,uj> for an inner product, find the elt u, such that
 484    <u,ui>=1 and <u,uj>=0 for i not = j, where i= number-of-element-in-basis."
 485   (setq cond0 ($list_matrix_entries
 486                 (ncmul* trace-matrix (setq vari (subseq $aaaa 0 (length basis))))))
 487   (setq actual-conditions   (cons '(mlist) (loop for v in (cdr cond0 )
 488                                                 for i from 1
 489                                                 when (eql i number-of-element-in-basis)
 490                                                 collecting (sub* v 1)
 491                                                 else collecting v)))
 492   (setq result ($fast_linsolve actual-conditions vari))
 493   ($sublis result ($general_sum basis $aaaa)))
 494
 495 (defun $linear_variables (eqns &optional constants &aux tem vari var-eqns)
 496   (cond (constants (setq constants (list-variables (st-rat constants)))))
 497   (setq eqns (st-rat eqns))
 498   (setq vari (list-variables eqns))
 499   (setq var-eqns (mapcar 'list-variables eqns))
 500   (loop for v in vari
 501         with mon = (list nil 1 1)
 502         do
 503         (loop for w in var-eqns
 504               for eqn in eqns
 505               when (member v w)
 506               do (setf (car mon) v) (setq tem (pcoeff eqn mon))
 507               (cond ((> (pdegree eqn v) 1)(setq vari (delete v vari)))
 508                     (t
 509                      (cond ((numberp tem))
 510                            ((and constants (every #'(lambda (var)(member var constants))
 511                                                   (list-variables tem))))
 512                            (t
 513                             (setq vari (delete v vari))))))))
 514   (values (cons '(mlist) (loop for v in vari collecting (get v 'disrep)))
 515           vari))
 516
 517 (defun $fast_linsolve_inconsistent_equations ( &optional (disrep t))
 518   (let ((eqns (find-extra-conditions (pv-the-sparse-matrix $poly_vector))))
 519     (cond (disrep (cons '(mlist) (mapcar 'new-disrep eqns)))
 520           (t eqns))))
 521
 522 (defun find-extra-conditions (sp-mat)
 523   (loop for i below (array-total-size (sp-column-used-in-row sp-mat))
 524      when( and  (null (aref (sp-column-used-in-row sp-mat) i))
 525                 (not (new-zerop (aref (sp-constants-column sp-mat) i))))
 526      collecting (aref (sp-constants-column sp-mat) i)))
 527
 528 (defun $my_phi(form)
 529   (setq form ($ratsimp form))
 530   (rota1 form))
 531
 532 (defun rota1 (form)
 533    (cond ((atom form) form)
 534          ((eql (caar form) 'mnctimes)
 535           (cons (car form) (cons (car (last form)) (butlast (cdr form) 1))))
 536          (t (cons (car form) (loop for v in (cdr form ) collect (rota1 v))))))
 537
 538 (defun $central_rels( lis)
 539  (cons '(mlist) (loop for v in (cdr lis)
 540          append
 541          (loop for w in (cdr ($append lis  $current_variables))
 542             collect (new-disrep (n- (n. v w) (n. w v)))))))
 543
 544 ;;Make lis1 commute with lis2
 545
 546 (defun $centralize_rels (lis1 lis2)
 547   (cons '(mlist)
 548   (loop for v in (cdr lis1)
 549          append
 550          (loop for w in (cdr lis2)
 551                 collect
 552                 (new-disrep (n- (n. v w) (n. w v)))))))
 553
 554 (defun $set_up (rels &optional n)
 555   (setq *previously-checked-pairs* nil)
 556   ($set_up_dot_simplifications rels)
 557     (setq $bil (cons '(mlist) *all-dotsimp-denoms*))
 558     (setq  *all-dotsimp-denoms* nil)
 559     ( format t "~%Beginning check overlaps")
 560      (prog1 ($check_overlaps n t nil)
 561  (setq $denoms (cons '(mlist) *all-dotsimp-denoms* ))))
 562
 563 (defun $half_quantum_rels (rel3 m)
 564   (let* ($dot_simplifications
 565          (vec ($sort ($list_variables ($list_nc_monomials rel3))))
 566          (matt ($generate_matrix '(lambda (i j) ($concat '$u i j)) m m))
 567          (joe ($sub_list `((mlist)
 568                            ((mequal)
 569                             ,vec ,($list_matrix_entries
 570                                    (ncmul* matt ($transpose vec)))))
 571                          rel3)))
 572     (format t "Using variables:" )(displa vec)
 573     (setq joe (simplifya joe nil))
 574     (displa joe)
 575     (let ((ctls ($centralize_rels ($list_matrix_entries matt) vec))
 576           dotted mon)
 577           (displa ctls)
 578       ($set_up ($append ctls rel3) 4)
 579       (setq dotted ($dotsimp joe))
 580       (setq mon ($mono vec 2))
 581       ($set_up (sub* mon ($firstn ($length mon) $aaaa)) 1)
 582       ($print $dot_simplifications)
 583       (setq dotted ($dotsimp dotted))
 584       (apply  '$append (cdr ($separate_parameters dotted "aa"))))))
 585
 586
 587
 588
 589
 590
 591 (defvar *nvars* 6)
 592 (defun list-ordered-pairs-with-repeat-count(i j n nrepeats)
 593   ;; return the list of 2*n integers all less than i
 594   ;; such that there are no more than NREPEATS.
 595   (let ((*nvars* i))
 596     (list-ordered-pairs-with-repeat-count1 i j n nrepeats)))
 597
 598 (defun repeat-ok1 (i j lis nrepeats &aux (n 0))
 599
 600   ;; return t if (append (list i j) lis) has no more than nrepeats
 601   (cond ((eql i j) (setq n 1)))
 602   (cond ((member i lis) (setq n (+ n 1))))
 603   (cond ((member j lis) (setq n (+ n 1))))
 604   (cond ((> n nrepeats) nil)
 605         (t (loop for x on lis when (member (car x) (cdr x))
 606                   do (setq n (1+ n)))
 607            (cond ((> n nrepeats) nil)
 608                  (t t)))))
 609
 610 (defun list-ordered-pairs-with-repeat-count1 (i j number-pairs number-repeats)
 611   (cond ((<= number-pairs 0) (list nil))
 612         (t
 613          (loop for ii to i
 614                 when ( < ii *nvars*)
 615                 append
 616                 (loop for jj below (if (eql ii i) j *nvars* )
 617                        append
 618                        (loop for u in (list-ordered-pairs-with-repeat-count1
 619                                         ii jj (- number-pairs 1) number-repeats)
 620                               when (repeat-ok1 ii jj u number-repeats)
 621                                       collect (append (list ii jj) u)))))))
 622
 623
 624
 625 (defun extract_cof1 (f vars)
 626   (cond ((atom f) (list f))
 627         ((unless (affine-polynomialp f)
 628             (progn  (setq f (numerator f)) nil)))
 629         ((member (car f) vars)
 630          (loop for (deg cof) on (cdr f) by #'cddr
 631                 nconc (extract_cof1 cof vars)))
 632         (t (list f))))
 633
 634 (defun $extract_linear_from_commutative (form vars)
 635   (assert  (listp form))
 636   (setq form (st-rat form))
 637   (setq vars (list-variables (st-rat vars)))
 638   (cons '(mlist)
 639         (mapcar #'new-disrep
 640                 (loop for f in form
 641                    appending
 642                      (extract_cof1 f vars)))))