share/affine/polyd.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 (defvar $last_answer nil)
  11
  12 (defun which-centrals-to-add (&aux *previously-checked-pairs* ($dot_simplifications))
  13   "Checks you aren't adding a trivial product of centrals etc.
  14 It affects the dot simplifications which become the dot simplifications
  15 modulo ones-to-add checked thru the highest degree of ones to add."
  16                                         ;  (unwind-protect
  17   (progn
  18     (loop for v in (cdr $centrals_so_far)
  19            until (eql (length ones-to-add) 3)
  20            do
  21            (cond (($zerop ($dotsimp v))(setq v 0))
  22                  (t (format t "~%Supposedly checking overlaps to deg ~a" ($nc_degree v))
  23                     (setq $rank_function (hilbert-modulo
  24                                           (append  ones-to-add (list v))))
  25                     (setq *previously-checked-pairs* nil)
  26                     ($check_overlaps ($nc_degree  v) t nil)))
  27            when (not ($zerop ($dotsimp v)))
  28            collecting v into ones-to-add
  29            and do ($add_relation_to_dot_simplifications v)
  30            finally (return ones-to-add))))
  31 ;    nil)
  32  ; (setq $dot_simplifications dot-simps)))
  33
  34 (defun find-number-of-centrals
  35        (n &optional(from-degree 2) (to-degree 10) (reset-centrals-so-far nil))
  36   "Finds at least  n centrals if there are n thru to-degree
  37 It also checks overlaps and so affects the dot_simplifications.  It stores
  38 them in $centrals_so_far.  If later asked it won't check below the centrals so
  39 far degree"
  40   (cond (reset-centrals-so-far (setq $centrals_so_far '((mlist)))))
  41   (cond ((> 1 (length $centrals_so_far))
  42          (setq from-degree (1+ (apply 'max from-degree
  43                                      (mapcar '$nc_degree (cdr $centrals_so_far)))))))
  44   (loop for i from (1+ from-degree) to (1+ to-degree)
  45         do ($check_overlaps i t nil (1- i ))
  46 ;       (show i )(format t "~~%*****************")
  47 ;        (displa $centrals_so_far)
  48         (loop for j from 1 to i
  49               when
  50               (and $rank_function
  51                    (not (eql ($length ($mono $current_variables  j))
  52                         (funcall $rank_function j))))
  53               do
  54               (format t "~%Using Current Variables of..")(displa $current_variables)
  55               (format t "   and rank function ~A." $rank_function)
  56               (throw 'check_a_case
  57                       (format nil "In spite of checking overlaps to degree ~A the monomial dimensions were ~A and the rank functions was ~A in degree ~A.  We conclude not finite global dimension." i ($length ($mono $current_variables  j))
  58                               (funcall $rank_function j) j)))
  59         (setq $centrals_so_far
  60               ($append $centrals_so_far
  61                        ($fast_central_elements
  62                          $current_variables (1- i))))
  63         (format t "~%Central elements thru degree ~A" (1- i))
  64         (displa $centrals_so_far)
  65         when (>= ($length $centrals_so_far) n)
  66         do
  67         (return 'done)))
  68
  69 (defun monomials-and-rank-function-agree (n &aux mono rank)
  70   (loop for i from 1 to n
  71         do
  72         (setq mono ($length ($mono $current_variables i))
  73               rank (funcall $rank_function i))
  74         (format t "~%In degree ~A the monomial dimension is ~A and the rank function ~A."
  75                 i mono rank)
  76         (cond ((not (eql mono rank))(return nil)))
  77         collecting mono into tem
  78         summing mono into total
  79         finally (push 1 tem) (push "total" tem) (push (incf total) tem)
  80         (format t "~%The total is ~A." total)(return  tem)))
  81
  82
  83 (eval-when
  84     (:load-toplevel)
  85   (cond ((not (member 'hilbert *all-rank-functions* :test #'eq))
  86          (push 'hilbert *all-rank-functions*))))
  87
  88 (defun $check_a_case
  89        (relations  &optional (to-degree 10) (rank-function '$global_dimension_3)
  90         (variables nil)
  91         &aux tem ones-to-add dims answer free-dot-simps full-dot-simps)
  92   "look for centrals thru degree to-degree then mod them out"
  93   (setq $centrals_so_far '((mlist)))
  94   (cond (variables (setq $current_variables variables))
  95         (t (setq $current_variables ($list_variables relations "X" "Y" "Z"))
  96            (loop for v in (cdr $current_variables)
  97                  when (or ($scalarp v) (search "%" (string v) :test #'char-equal))
  98                  do (setq $current_variables (delete v $current_variables :test #'equal)))))
  99   (show $current_variables)
 100   (funcall '$set_up_dot_simplifications  relations )
 101   (setq full-dot-simps $dot_simplifications)
 102   (setq answer
 103         (catch 'check_a_case
 104           (setq free-dot-simps $dot_simplifications)
 105           (loop for i from 3 until (or (eql i 6) (eql (length ones-to-add) 3))
 106                 ;; 6 ;;enough?
 107                 do
 108                 (setq $dot_simplifications full-dot-simps)
 109                 (setq $rank_function rank-function)
 110                 (setq *previously-checked-pairs* nil)
 111                 ;;you've changed them in ones-to-add
 112                 (find-number-of-centrals i 2 (1+ to-degree))
 113                 ;;gets centrals thru to-degree.
 114                 ;;may throw if not finite global dim.
 115                 (setq full-dot-simps $dot_simplifications)
 116                 ;;should be free ones to deg i
 117                 (setq ones-to-add (which-centrals-to-add))
 118                 ;; dotsimps from which.. stay
 119                 finally
 120                 (cond ((>= (length ones-to-add) 3)
 121                        (setq to-degree ($nc_degree (car (last ones-to-add))))))
 122                 (setq $rank_function
 123                       (hilbert-modulo ones-to-add))
 124                 (cond (ones-to-add
 125                        (cond((setq dims
 126                                    (monomials-and-rank-function-agree
 127                                      (floor (1+ (* 2.2 ($nc_degree (car (last ones-to-add))))))))
 128                              (return nil))))
 129                       (t (return (setq answer
 130                                        (list '(mlist) relations
 131                                              (format nil "There are no centrals up to degree ~A"
 132                                                      to-degree)
 133                                              ($ratsimp $dot_simplifications))))))
 134                 ;;hopefully we don't get here and tessing worked.
 135                 (format t "~%~%HOW DID WE GET HERE?~%")
 136                 (setq *previously-checked-pairs* nil)
 137                 (funcall '$set_up_dot_simplifications
 138                          (append relations ones-to-add ))
 139                 ($check_overlaps 12 t nil))))
 140   (cond (answer answer)
 141         (t (cond ((null dims)
 142                   (setq dims
 143                         (with-no-query
 144                           (cdr ($monomial_dimensions 15))))))
 145            (setq $last_answer
 146                  (list '(mlist)
 147                        "Relations:" relations
 148                        "Centrals added:" (cons '(mlist) ones-to-add)
 149                        "Monomial dimensions:" (cons '(mlist) dims)
 150                        (format nil "All Centrals to Degree ~A:" to-degree)
 151                        $centrals_so_far
 152                        "Dot Simplifications:"
 153                        ($ratsimp $dot_simplifications) ))))
 154   (format t "~%Conclusion:~%The algebra is a finite module of rank ~D over a polynomial ring.
 155  The polynomial ring is generated by elements in degrees ~D .  The actual generators are "  (second (nth 6 $last_answer)) (mapcar '$nc_degree (cdr $centrals_so_far)))
 156   (displa $centrals_so_far)
 157   (format t "~%The dimension of the spaces of monomials in degrees 0 through ~D are~%"
 158           (1-  (length (setq tem (cdddr (delete 0 (nth 6 $last_answer) :test #'equal))))))
 159   (format t "  ~A "  tem)
 160   'qed)
 161
 162 (defvar $last_answer nil)
 163 (defun te (&aux tem)
 164  (format t "~%Conclusion:~%The algebra is a finite module of rank ~D over a polynomial with 3 generators.
 165  The polynomial ring is generated by 3 elements in degrees ~D .  The generators are "  (second (nth 6 $last_answer)) (mapcar '$nc_degree (cdr $centrals_so_far)))
 166 ;  (displa $centrals_so_far)
 167   (format t "~%The dimension of the spaces of monomials in degrees 0 through ~D are~%"
 168         (1-  (length (setq tem (cdddr (delete 0 (nth 6 $last_answer) :test #'equal))))))
 169  (format t "  ~A "  tem) 'qed)
 170
 171 (defun $check_programs (list-relations)
 172   (loop for v in (cdr list-relations)
 173         do (reset-vgp)
 174         collecting ($check_a_case v) into tem
 175         finally (return (cons '(mlist) tem))))
 176
 177
 178 (defun new-convert-relation-to-dot-simp (relat  &optional (ordering $order_function)
 179                                        (then-simplify nil) &aux
 180                                      cof worst  )
 181           (setq relat ($dotsimp relat))
 182           (setq relat (num (cdr relat)))
 183           (setq worst (get (car relat) 'disrep)
 184                 cof (third  relat))
 185           (setq relat  (nred relat cof))
 186           (setq relat (n-  worst relat))     ;;; should take cdddr etc. not subtractt
 187           (setq relat (header-poly relat))
 188           (format t "~%Adding  ")
 189           (displa (cons '(mlist simp) (list worst relat)))
 190           (format t "    to dot_simplifications")
 191           (setq $dot_simplifications (append $dot_simplifications (list worst relat)))
 192           ($sort_dot_simplifications ordering)
 193           (cond (then-simplify ($simplify_dot_simplifications)))
 194           $dot_simplifications)
 195
 196 ;;This is only for $new_fast_dotsimp  tried to avoid calling itself and
 197 ;;added an unwind protect so that the dotsimps don't get screwed up if
 198 ;;you abort part way thru.
 199 (defun $simplify_dot_simplifications (&optional (from-degree 0) &aux
 200                                       relat resimplify dot-simps leng)
 201   (check-arg $new_fast_dotsimp (eq $new_fast_dotsimp t) "Use the old $simplify_dot_simplifications in polynomials.lisp for other dotsimps")
 202    ($sort_dot_simplifications)
 203      (format t "~%starting to resimplify dot simplifications..")
 204      (cond ((< (setq leng ($length $dot_simplifications)) 10) (displa $dot_simplifications))
 205            (t (format t "~%There are ~A of them." (truncate leng 2))))
 206
 207    (setq dot-simps $dot_simplifications)
 208    (loop named sue for v on (cdr $dot_simplifications) by #'cddr
 209          for i from 1 by 2
 210          do
 211          (unwind-protect
 212            (progn
 213              (setq $dot_simplifications (append (firstn i $dot_simplifications) (cddr v)))
 214              (cond ((and ;(not ($zerop (second v)))
 215                          (or ($must_replacep (car v))($must_replacep (second v))))
 216                     (setq relat (cond ($new_fast_dotsimp
 217                                        (header-poly (n- (car v) (second v))))
 218                                       (t (vsub* (car v) (second v)))))
 219                     (setq relat ($dotsimp relat))
 220                     (cond (($zerop relat))
 221                           (t (new-convert-relation-to-dot-simp relat )))
 222                     (setq dot-simps $dot_simplifications)
 223                     (setq resimplify t)
 224                     (return-from sue 'start-over))))
 225            (setq $dot_simplifications dot-simps)))
 226    (cond (resimplify ($simplify_dot_simplifications))
 227          (t (format t "~%They were OK")))
 228    $dot_simplifications)
 229
 230 (defun $monomial_and_degree_lessp (x y)
 231   (let ((x-deg ($nc_degree x)) (y-deg ($nc_degree y)))
 232     (cond ((> x-deg y-deg) t)
 233           ((eq x-deg y-deg)($monomial_alphalessp x y))
 234           (t nil))))
 235
 236 (defun mono-dimension (n) ($length ($mono $current_variables n)))
 237
 238 (defun hilbert (n &rest l)
 239   (cond ((< n 0) 0)
 240         ((eql n 0) 1)
 241         ((null l)($global_dimension_3 n))
 242         (t (- (apply 'hilbert n (cdr l))
 243               (apply 'hilbert (- n (car l)) (cdr l))))))
 244
 245 ;(defun $global_dimension_3 (n)
 246 ;  (case ($length $current_variables)
 247 ;    (2 (rank-dimension-three n))
 248 ;    (3 (polynomial-ring-1-1-1 n))))
 249
 250
 251 (defun our-dimensions (n)
 252  "for 2 3 and 4 variables we have unique sequences so we plug them in here."
 253   (case (length (cdr $current_variables))
 254         (0 (fsignal "no variables"))
 255         (1 1)
 256         (2 (dimension-from-sequence '(1 2 2 1) '(1 2 1) n 2))
 257         (3 (dimension-from-sequence '(1 3 3 1) '(1 1 1) n 3))
 258         (4  (dimension-from-sequence '(1 4 6 4 1) '(1 1 1 1) n 4))
 259         (t (fsignal "not handled yet"))))
 260
 261 (defun $global_dimension_3 (n) (our-dimensions n))
 262
 263 (defun dimension-from-sequence
 264        (list-powers list-degree-maps degree number-variables &aux ( answ 0))
 265   "the sequences start from the right tail."
 266   (check-arg list-powers (eql (car list-powers) 1) "first should be one")
 267   (cond ((< degree 0) 0)
 268         ((zerop degree) 1)
 269         ((eql 1 degree) number-variables)
 270         ;;(< degree degree-relations) (expt  number-variables degree))
 271         (t
 272          (loop for deg-map in  list-degree-maps
 273                for u in (cdr list-powers)
 274                for i from 0
 275                with j = degree
 276                while (>= j 0)
 277                do (setq j (- j deg-map))        ;(show (list degree j answ))
 278                (setq answ (+ answ (* (expt -1 i) u  (dimension-from-sequence list-powers list-degree-maps j
 279                                                                           number-variables ))))
 280 ;              (show answ)(show u)
 281                finally (return (abs answ))))))
 282
 283 (defun $sum_hilbert (&rest l &aux tem)
 284   (loop for i below 100  until (zerop (setq tem (apply #'hilbert i 3 l)))
 285         summing tem))
 286
 287 ;;for creating a one argument function that can be funcalled
 288 (defun $hilbert_modulo (l)
 289   (check-arg l $listp "macsyma list")
 290   (funcall 'hilbert-modulo (cdr l)))
 291
 292 (defun hilbert-modulo (l)
 293   (cond ((or (null l) (numberp (car l)))
 294          `(lambda (n) (apply 'hilbert n ',(copy-list l))))
 295         (t (hilbert-modulo  (mapcar #'$nc_degree l)))))
 296
 297 (defun $tes (relat)
 298   (new-convert-relation-to-dot-simp relat $order_function t))
 299
 300 (defun $add_relation_to_dot_simplifications (relat )
 301   (new-convert-relation-to-dot-simp relat $order_function t))
 302
 303 (defun $rat_the_dot_simplifications (&aux res v)
 304  "Does the new-rat for New_fast_dotsimp"
 305   (loop for w on (cdr $dot_simplifications) by #'cddr
 306         do (push (car w) res)
 307         (setq v (cadr w))
 308         (push (if ($ratp v) (header-poly (new-rat (ratdisrep v))) (header-poly  (new-rat v)))
 309                   res)
 310         )
 311   (setq $dot_simplifications (cons '(mlist) (nreverse res))))
 312
 313 ;
 314 ;
 315 ;
 316 ;
 317 ;(defun x$-macro-read (&optional  subchar
 318 ;                     (stream *standard-input*) arg &aux ( meval-flag t))
 319 ;
 320 ;  (declare (ignore subchar)) subchar
 321 ;  (setf (fill-pointer *my-stream* ) 0)
 322 ;  (cond ((eq #\$ (tyipeek t stream))(send stream :tyi)
 323 ;        (setq meval-flag nil)))
 324 ;  (with-output-to-string (st *my-stream*)
 325 ;    (let (char)
 326 ;      (loop while (not (eql char #\$))
 327 ;           do
 328 ;           (setq char (send stream :tyi))
 329 ;           (send st :tyo char)
 330 ;           finally (cond ((not (eql  char #\$))
 331 ;                          (zwei:barf "There was no matching $" ))))))
 332 ;  (cond (meval-flag
 333 ;        (list 'meval* (list 'quote   (parse-string *my-stream*))))
 334 ;       (t  (list 'quote   (parse-string *my-stream*)))))
 335 ;
 336 ;(si:SET-SYNTAX-/#-MACRO-CHAR #\$ 'x$-macro-read)
 337
 338 (defun $fast_central_elements_given_commutator (variables deg comut   used_var
 339                                                &aux  f unknowns eqns tem1 parameters answer )
 340   (setq variables (cons '(mlist) (sort (cdr variables) $order_function)))
 341   (setq $commutators nil $centralizers nil)
 342   (let* ((monoms ($mono variables deg))
 343 ;       (monoms-higher ($mono variables (1+ deg)))
 344           monoms-higher)
 345     (setq f ($general_sum monoms $aaaa))
 346     (setq unknowns ($list_variables f "aa" "par"))
 347 ;    (setq f prev-comut)
 348     (setq tem1 `((mlist) ,comut))                    ;;; (list '(mlist) ($com f v))
 349 ;         (setq $commutators (append `((mlist) ,v ,(second tem1)) (cdr $commutators)))
 350 ;         (setq eqns ($extract_linear_equations tem1 monoms-higher))
 351     (setq monoms-higher ($list_nc_monomials tem1))
 352     (setq eqns ($extract_linear_equations tem1 monoms-higher ))
 353     (show unknowns)
 354      (break t)
 355     (setq answer ($fast_linsolve eqns unknowns))
 356     (setq f ($ratsimp ($sublis answer f)))
 357     (setq $centralizers (append `((mlist) xy ,f) (cdr $commutators)))
 358     (displa f)
 359     (loop for v in (nthcdr ($length used_var) (cdr variables))
 360           do
 361           (setq unknowns ($list_variables f "aa" "par"))
 362           (setq parameters (loop for vv in (cdr unknowns)
 363                                  when (search "par" (string vv) :test #'char-equal)
 364                                  collecting vv))
 365           (setq tem1 (cdr $aaaa))
 366           (loop for vv in parameters
 367                 do (loop while tem1
 368                          when (not (member (car tem1) unknowns :test #'eq))
 369                          do (setq f (subst (car tem1) vv f))
 370                          (setq unknowns (subst (car tem1) vv unknowns))
 371                          (format t "~%Replacing ~A by ~A in f." vv (car tem1))
 372                          (setq tem1 (cdr tem1)) (return 'done)
 373                          do (setq tem1 (cdr tem1))))
 374
 375           (setq tem1 (list '(mlist) ($com f v)))
 376           (setq $commutators (append `((mlist) ,v ,(second tem1)) (cdr $commutators)))
 377 ;         (setq eqns ($extract_linear_equations tem1 monoms-higher))
 378           (setq eqns ($extract_linear_equations tem1 ($list_nc_monomials tem1)))
 379           (show unknowns)
 380           (setq answer ($fast_linsolve eqns unknowns))
 381           (setq f ($ratsimp ($sublis answer f)))
 382           (setq $centralizers (append `((mlist) ,v ,f) (cdr $commutators)))
 383           (displa f)
 384           finally (return ($separate_parameters f)))))
 385
 386
 387
 388 (defun $read_Lisp_string (a-string)
 389   (eval (read-from-string (string-trim "&" a-string))))
 390
 391
 392 (defun $try (u &aux tem centrals)
 393 ;  (funcall '$set_up_dot_simplifications (append $relations (list u)))
 394   ($add_relation_to_dot_simplifications u)
 395   (setq centrals nil)
 396   (setq *previously-checked-pairs* nil)
 397   (setq $rank_function nil)
 398   ($check_overlaps 10 t nil)
 399   (loop for i from 1 below 8
 400         when (setq tem (cdr  ($fast_central_elements $current_variables i)))
 401         do (loop for v in tem do ($add_relation_to_dot_simplifications v))
 402         (setq centrals (append centrals tem))
 403         ($check_overlaps 10 t nil)
 404         finally (return (cons '(mlist) centrals))))
 405
 406 (defun list-variables (expr &aux *all-vars*)
 407     (declare (special *all-vars*))
 408   (list-variables1 expr)
 409   *all-vars*)
 410
 411 (defun list-variables1 (expr)
 412   (declare (special *all-vars*))
 413    (cond   ((atom expr) nil)
 414            ((affine-polynomialp expr)(pushnew (p-var expr) *all-vars*)
 415             (loop for (deg cof) on (cdr expr) by #'cddr
 416                   do (list-variables1 cof)))
 417            ((rational-functionp expr)(list-variables1 (car expr))
 418             (list-variables1 (cdr expr)))
 419            (t (loop for v in expr do (list-variables1 v)))))
 420
 421 (defun reset-rat (tree &aux tem)
 422   (cond ((atom tree)
 423          (cond ((and (symbolp tree)
 424                      (setq tem (get tree 'disrep)))
 425                 (add-newvar tem))
 426                (t tree)))
 427         (t (setf (car tree) (reset-rat (car tree)))
 428            (setf (cdr tree) (reset-rat (cdr tree))))))
 429
 430 (defun npfactor (p  )
 431   (let ((varl (list-variables p)))
 432     (let ((genvar (reverse (sort varl #'pointergp))))
 433       (pfactor p))))
 434
 435 (defun nplcm (x y &aux (genvar *genvar*))
 436   (plcm x y))
 437
 438 (defun npgcd (f g)
 439   (let ((varl (list-variables (list f g))))
 440     (let ((genvar (reverse (sort varl 'pointergp))))
 441       (pgcd f g))))