share/affine/polyb.lisp

   1 ;;; -*- mode: lisp; package: cl-maxima; syntax: common-lisp -*-
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;                                                                    ;;;;;
   4 ;;;     Copyright (c) 1984 by William Schelter,University of Texas     ;;;;;
   5 ;;;     All rights reserved                                            ;;;;;
   6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   7
   8 (in-package :maxima)
   9
  10 (defun $circle_times (&rest l)
  11   (case  (length l)
  12     (0 1)
  13     (1 (car l))
  14     (2 (let (firs sec (u (car l)) (v (second l)))
  15          (setq firs (sub* (mul* (nth 1 u) (nth 1 v))
  16                           (mul* (nth 2 u) (nth 2 v))))
  17          (setq sec (add* (mul* (nth 1 u) (nth 2 v))
  18                          (mul* (nth 2 u) (nth 1 v))))
  19        (list '(mlist simp) firs sec)))
  20     (otherwise
  21      ($circle_times (car l) (apply '$circle_times (cdr l))))))
  22
  23 (defun sort-grouped-list (a-list group-size pred &aux if-necessary answer tem (prev (car a-list)))
  24   (loop for v in a-list
  25      for i from 0
  26      when (zerop (mod i group-size))
  27      do (cond ((funcall pred v prev)
  28                (setq if-necessary t)
  29                (return t)))
  30        (setq prev v))
  31   (cond (if-necessary
  32          (loop  with v = a-list
  33             while v
  34             collect v into all
  35             do (setf v (nthcdr group-size v))
  36             finally
  37               (setq all (sort all pred :key 'car ))
  38               (return
  39                 (loop for v in all
  40                    nconc
  41                      (loop for w in v for i below group-size
  42                         collect w)))))
  43         (t a-list)))
  44
  45 (defun $sort_dot_simplifications (&optional (pred $order_function) &aux rev-pred)
  46   (cond ((eq pred '$monomial_alphalessp)(setq rev-pred '$monomial_alphagreatp))
  47         (t (setq rev-pred #'(lambda (x y)
  48                               (and (not (funcall pred x y))(not (equal x y)))))))
  49   (rplacd $dot_simplifications
  50           (sort-grouped-list (cdr $dot_simplifications) 2 rev-pred)))
  51
  52 ;see polyd.lisp for current def
  53 ;(defun $simplify_dot_simplifications (&optional (from-degree 0) &aux temp1 repl temp2 temp3
  54 ;                                     relat
  55 ;                                     (dot-simps $dot_simplifications))
  56 ;
  57 ;  ($sort_dot_simplifications)
  58 ;  (loop for term in (cdr $dot_simplifications)
  59 ;       for i from 0
  60 ;       when (evenp i)
  61 ;       do
  62 ;       (let (($dot_simplifications
  63 ;               (append (firstn (1+ i) dot-simps)
  64 ;                       (nthcdr (+ 3 i) dot-simps))))
  65 ;         (cond (($must_replacep term)
  66 ;                 (setq repl (nth (1+ i)  (cdr dot-simps)))
  67 ;
  68 ;                 (setq relat
  69 ;                      (cond ($new_fast_dotsimp
  70 ;                             (header-poly (n- repl term)))
  71 ;                            (t (vsub* repl term))))
  72 ;                 (setq dot-simps $dot_simplifications)
  73 ;                (return 'start-over))))
  74 ;       when (oddp i)
  75 ;       do
  76 ;       (setq temp2 term)
  77 ;       (show temp2)
  78 ;;;;    (cond ((and
  79 ;;;;;            (>= ($nc_degree temp2) from-degree)
  80 ;;;;             ($must_replacep temp2))
  81 ;
  82 ;              (loop until  (null ($must_replacep temp2))
  83 ;                    do
  84 ;                    (break 'here)(setq temp1 temp2)
  85 ;                    (cond ($free_dot (setq temp2 ($dotsimp temp2)))
  86 ;                          (t
  87 ;                           (setq temp3 (meval* temp1))
  88 ;                           (cond ((not (equal temp1 temp3))
  89 ;                                  (setq temp2 ($ratsimp temp3))))))
  90 ;                    finally
  91 ;                    (setf (nth (1+ i) $dot_simplifications) temp2)))
  92 ;  (cond (relat (setq $dot_simplifications dot-simps)
  93 ;              (convert-relation-to-dot-simp ($dotsimp relat)))
  94 ;       (t
  95 ;        $dot_simplifications))
  96 ;   (cond ($fast_dotsimp ($rat_the_dot_simplifications))
  97 ;        (t $dot_simplifications)))
  98 ;
  99 ;;;second The above definition worked for not fast_dotsimp
 100 ;
 101 ;
 102 ;(defun $simplify_dot_simplifications (&optional (from-degree 0) &aux
 103 ;                                     relat
 104 ;                                     (dot-simps $dot_simplifications))
 105 ;   ($sort_dot_simplifications)
 106 ;   (loop for v on (cdr $dot_simplifications) by #'cddr
 107 ;        for i from 1 by 2
 108 ;        do
 109 ;         (setq dot-simps $dot_simplifications)
 110 ;        (setq $dot_simplifications (append (firstn i $dot_simplifications) (cddr v)))
 111 ;        (cond ((and (not ($zerop (second v)))
 112 ;                    (or ($must_replacep (car v))($must_replacep (second v))))
 113 ;               (setq relat (cond ($new_fast_dotsimp
 114 ;                                  (header-poly (n- (car v) (second v))))
 115 ;                                 (t (vsub* (car v) (second v)))))
 116 ;               (setq relat ($dotsimp relat))
 117 ;               (cond (($zerop relat))
 118 ;                     (t (convert-relation-to-dot-simp relat )))     ;; ;(copy-tree relat))))
 119 ;               ($simplify_dot_simplifications from-degree)
 120 ;               (return 'start-over))
 121 ;              (t (setq $dot_simplifications dot-simps))))
 122 ;   $dot_simplifications)
 123
 124
 125
 126 (defmfun ncmul2* (x y)
 127   (simplifya `((mnctimes) ,x ,y) nil))
 128
 129 (defun contains-a-zero-replacement (form &aux u)
 130   (or (in-nth-power-radical form $radical_nilpotent_of_order)
 131       (cond ((atom form)
 132              (loop for v on (cdr $dot_simplifications) by #'cddr
 133                    when (and (eql 0 (second v)) (eq (car v) form))
 134                    do (return t)))
 135             (t (loop for v on (cdr $dot_simplifications) by #'cddr
 136                      when (eql 0 (second v))
 137                      do (cond ((atom (setq u (car v)))
 138                                (cond ((member u form :test #'equal) (return t))))
 139                               (t (cond ((ordered-sublist (cdr u) form)(return t))))))))))
 140
 141 (defun  contains-a-replacement (form &aux tem1 u)
 142   "returns t if (ncmuln form nil) would cause a replacement  when simplified (form is
 143  not a product"
 144   (setq tem1 (cdr $dot_simplifications))
 145   (cond ((atom form)
 146          (loop for u in tem1
 147             for i from 1
 148             when (and (oddp i)(eq u form))
 149             do(return t)))
 150         (t
 151          (loop while tem1
 152                do (setq u (car tem1))
 153                (setq tem1 (cddr tem1))
 154                (cond ((atom u)(cond ((member u form :test #'equal) (return t))))
 155                      (t (cond ((ordered-sublist (cdr u) form )(return t)))))))))
 156
 157
 158
 159 (defvar *already-found* nil)
 160
 161 (defun $must_replacep (expression &aux(ptype (poly-type expression)) poly)
 162   (setq *already-found* nil)
 163   (case ptype
 164     (:rational-function (setq poly (num expression)))
 165     (:polynomial (setq poly expression))
 166     (:$rat (cond ($new_fast_dotsimp (setq poly (num (cdr expression)))))))
 167   (cond (poly
 168          (cond ((numberp poly) nil)
 169                (t
 170                 (loop while poly
 171                       when ($must_replacep (get (car poly) 'disrep))
 172                       do (return t)
 173                       when (numberp       (setq poly (fifth poly)))
 174                       do (return nil)))))
 175         (($ratp expression)
 176          (cond ((member (caadr expression) *genvar* :test #'eq)
 177                 ($must_replacep (cadr expression)))
 178                (t
 179                 (fsignal "this cre form should not be here, since its leading monom is not in *genvar*"))))
 180         (t
 181          (catch 'must-replace (must-replacep1 expression)))))
 182
 183 (defun must-replacep1-by-zero ( expression)
 184   (cond (*already-found* nil)
 185         ((atom expression) (cond ((contains-a-zero-replacement expression)
 186                                   (throw 'must-replace t))))
 187         ((atom (car expression))(must-replacep1-by-zero (cdr expression)))
 188         ((eq (caar expression) 'mnctimes)(cond ((contains-a-zero-replacement (cdr expression))
 189                                                 (setq *already-found* t)
 190                                                 (throw 'must-replace t) )))
 191         (t (must-replacep1-by-zero (car expression))(must-replacep1-by-zero
 192                                                       (cdr expression)))))
 193
 194 (defun must-replacep1 ( expression)
 195   (cond (*already-found* nil)
 196         ((atom expression) (cond ((contains-a-replacement expression)
 197                                   (throw 'must-replace t))))
 198         ((atom (car expression))(must-replacep1 (cdr expression)))
 199         ((eq (caar expression) 'mnctimes)(cond ((contains-a-replacement (cdr expression))
 200                                                 (setq *already-found* t)
 201                                                 (throw 'must-replace t) )))
 202         (t (must-replacep1 (car expression))(must-replacep1 (cdr expression)))))
 203
 204 ;(defun $check_associative (a b c &aux tem)
 205 ;  (setq tem (sub* ($ratsimp (ncmul*  (ncmul* a b) c)) ($ratsimp (ncmul* a (ncmul* b c))))))
 206 ;(defun $check_associative (a b c &aux tem)
 207 ;  (setq tem  ($ratsimp (sub*  ($simp_ncmul
 208 ;                              ($simp_ncmul a b) c)
 209 ;                            ($simp_ncmul a ($simp_ncmul b c))))))
 210
 211 ;(defun mread-noprompt (&rest read-args)
 212 ;  (let ((*mread-prompt* ""))
 213 ;  (caddr (apply #'mread read-args))))
 214
 215 (defun remove-header (x)
 216        (cond ((numberp x) x)
 217              (t (cdr x))))
 218
 219 (defun new-rat-ncmul1 (a b c &aux tem)
 220   (cond ((numberp b)
 221          (cond (($zerop (setq tem (ncmul* a b c))) 0)
 222                (t (new-rat tem))))
 223         ((affine-polynomialp b)(poly-ncmul1 a b c))
 224         ((rational-functionp b)(cons (poly-ncmul1 a (num b) c) (denom b)))
 225         (($ratp b)(new-rat-ncmul1 a (cdr b) c))
 226         (t (new-rat-ncmul a (new-rat b) c))))
 227 ;(defun $check_associative (a b c &aux tem term1 term2 answer)
 228 ;  (cond ($free_dot
 229 ;        (cond ($fast_dotsimp
 230 ;               (cond
 231 ;                 ($new_fast_dotsimp
 232 ;                  (setq term1 (new-rat-ncmul1 1 (remove-header ($dotsimp  (ncmul* a b))) c))
 233 ;                  (setq term2 (new-rat-ncmul1  a (remove-header ($dotsimp  (ncmul* b c))) 1))
 234 ;                  (setq tem  (n- term1 term2))
 235 ;;                 (check-rat-order tem)
 236 ;                  (setq answer ($dotsimp (header-poly tem)))
 237 ;                  (format t "~%Here is the associator")
 238 ;                  (displa answer)
 239 ;                  answer)
 240 ;                 (t
 241 ;
 242 ;                  (setq term1 (ncmul1 1 ($dotsimp ($vrat (ncmul* a b))) c))
 243 ;                  (setq term2 (ncmul1  a ($dotsimp ($vrat (ncmul* b c))) 1))
 244 ;                  (setq tem  (vsub* term1 term2))
 245 ;                  (check-rat-order tem)
 246 ;                  ($dotsimp tem))))))
 247 ;
 248 ;        (t
 249 ;          (setq tem  ($ratsimp
 250 ;                       ($dotsimp ($ratsimp
 251 ;                                   (sub*  (ncmul*
 252 ;                                            ($simp_ncmul a b) c)
 253 ;                                          (ncmul* a ($simp_ncmul b c))))))))))
 254
 255 ;;;works out a little simpler in the new notation. ;;see new-dotsimp.lisp
 256 ;(defun $check_associative (a b c &aux tem tem1 tem2)
 257 ;;  (show (dotsimp (n. a b)))
 258 ;;  (show (dotsimp (n. b c)))
 259 ;;  (show (n. a (dotsimp (n. b c))))
 260 ;;  (show (setq hee(n.  (dotsimp (n. a b)) c)))
 261 ;  (setq tem1 (dotsimp(n. (dotsimp (n. a b)) c)))
 262 ;  (setq tem2 (dotsimp(n. a (dotsimp (n. b c)))))
 263 ;;  (show tem1 tem2)
 264 ;  (setq tem (n- tem1 tem2))
 265 ;  (cond ((pzerop tem) tem)
 266 ;       (t (header-poly tem))))
 267
 268 ;(defmacro simp-zerop (n )
 269 ;  `(cond ((numberp ,n) (zerop ,n))
 270 ;        ((atom ,n) nil)
 271 ;        (t (cond ((member (caar ,n) '(mrat rat) :test #'eq)
 272 ;                  (cond ((eq (second (car ,n)) 'simp)
 273 ;                         (equal (cdr ,n) (rzero)))
 274 ;                        (t (zerop (setq ,n ($ratsimp ,n)))
 275 ;                           (format t "~having to Ratsimp ~A" ,n))))
 276 ;                  (t (and (numberp (setq ,n ($ratsimp ,n)))(zerop ,n)))))))
 277
 278 (defvar *previously-checked-pairs* nil)
 279 (defvar $global_dimension_three nil)
 280 (defvar $rank_function nil)
 281
 282 (defun two-times-n (n) (cond ((> n 4) (* 2 n))
 283                              (t 0)))
 284 (defun rank-dimension-three-modulo-cubic (n)
 285   (cond ((< n  4) 0)
 286         (t (- (rank-dimension-three n) (rank-dimension-three (- n 3))))))
 287
 288 (defun rank-1-1-1-9 (n)
 289   (cond ((< n 9) (polynomial-ring-1-1-1 n))
 290         (t
 291          (- (* n 9)  27))))
 292
 293 (defvar *all-rank-functions* '($global_dimension_3
 294                                rank-dimension-three
 295                                two-times-n rank-dimension-three-modulo-cubic
 296                                polynomial-ring-1-1-1 three-times-n
 297                                $standard_gorenstein rank-1-1-1-9 ))
 298
 299 (defun $check_overlaps (up-to-degree &optional (add-to-simps nil)
 300                                                (maybe-reset t) (from-degree nil)
 301                                                &aux tem to-replace test-list tem2 bef-overlap
 302                                                lowest-degree   deg
 303                                                tem1 assoc-list (ok t)
 304                                                list-of-degrees-to-mod-out
 305                                                dim reset-monomials)
 306   (or *previously-checked-pairs*
 307       (setq *previously-checked-pairs*
 308             (make-hash-table :test 'equal)))
 309   (cond (maybe-reset
 310          (cond ((y-or-n-p "Reset *previously-checked-pairs* ?")
 311                 (or *previously-checked-pairs*
 312                     (setq *previously-checked-pairs*
 313                           (make-hash-table :test 'equal)))
 314                 (clrhash *previously-checked-pairs*)
 315                 ))
 316          (cond ((y-or-n-p "use Hilbert for rank function")
 317                 (user-supply list-of-degrees-to-mod-out)
 318                 (setq $rank_function (hilbert-modulo list-of-degrees-to-mod-out)))
 319                ((y-or-n-p "Set the $Rank_function to a non nil value?")
 320                 (loop for v in *all-rank-functions*
 321                       do
 322                       (format t "~%Use ~A ?" v)
 323
 324                       (cond ((y-or-n-p ) (setq $rank_function v)(return 'done))
 325                             (t nil))
 326                       finally (format t
 327                                       "~%Enter a value for $rank_function:")
 328                       (setq $rank_function  (read))))
 329                (t (setq $rank_function nil)))
 330          (format t "~%~%The current variables are:")(displa $current_variables)
 331          (cond ((not (y-or-n-p "Keep the current variables?"))
 332                 (setq reset-monomials t)(format t "~%Enter a Macsyma list:")
 333                 (setq $current_variables  (mread-noprompt *standard-input*)))))
 334
 335         (t (format t "
 336 ~%Starting to check overlaps without resetting *previously-checked-pairs*")))
 337   (show $rank_function)
 338   (setq to-replace
 339         (loop for mon in $dot_simplifications
 340               for i from 0
 341               when (and (oddp i) (not (atom mon)))
 342               collecting mon))
 343   (loop for mon in to-replace
 344         minimize (setq deg ($nc_degree mon)) into the-min
 345         maximize deg into the-max
 346         finally (setq lowest-degree the-min)
 347         (cond ((not (numberp up-to-degree))
 348                (setq up-to-degree (1- (* 2 the-max))))))
 349
 350   (cond (from-degree (setq lowest-degree  from-degree)))
 351  (setq deg lowest-degree)
 352
 353
 354   (loop named angela
 355     while (<= deg up-to-degree)
 356     for deg from (1+ lowest-degree) to up-to-degree
 357
 358     do
 359     (cond ($rank_function
 360         (format t "~%current variables are")
 361         (displa $current_variables)
 362         (loop while (< deg (1+ up-to-degree))
 363               do
 364               (setq dim (1- (length
 365                                 ($mono $current_variables  deg reset-monomials))))
 366               (format t
 367                       "~%There are ~A independent monomials in degree ~A and rank function is ~A" dim deg (funcall $rank_function deg))
 368               (cond ((<= dim (funcall $rank_function  deg))
 369                      (setq lowest-degree  (min up-to-degree deg)) (setq deg (1+ deg)))
 370                     (t (return 'done))))))
 371     (cond ((<= deg up-to-degree)
 372     (loop
 373       for right1 in to-replace
 374       do
 375       (setq bef-overlap (- deg ($nc_degree right1)))
 376       (loop
 377         for left in to-replace
 378         do
 379         (loop while (< deg-so-far bef-overlap)
 380               for v in (cdr left)
 381               for ii from 2
 382               summing ($nc_degree v) into deg-so-far
 383               do
 384               (cond ((and (= deg-so-far bef-overlap)
 385                           (setq tem2 (nthcdr ii left))
 386                           (initial-equal tem2 (cdr right1)))
 387                      (setq tem (nthcdr (- (length left)
 388                                           ii) (cdr right1)))
 389                      (setq test-list (append (cddr left) tem))
 390                      (cond ((or
 391                               (contains-a-replacement
 392                                 (cddr (butlast test-list)))
 393                               (ordered-pair-in-list left right1
 394                                                     *previously-checked-pairs* ))
 395                             nil )
 396                            (t  (setf (gethash
 397                                               (list left right1) *previously-checked-pairs*)
 398                                     t)
 399                               (setq assoc-list
 400                                     (list (ncmuln (cdr (subseq left 0 ii)) t)
 401                                           (ncmuln (nthcdr ii left) t)
 402                                           (ncmuln tem t)))
 403                               (format t "~%Checking the overlap for")
 404                               (displa (list '(mlist simp) left right1))
 405                               (setq tem1 (apply  '$check_associative
 406                                                  assoc-list))
 407                               (cond (($zerop tem1)
 408                                      (format t "~%The overlap was associative."))
 409                                     (t
 410                                      (setq ok nil)
 411                                      (format t "~%The overlap was not associative.")
 412                                      (displa (cons '(mlist simp) assoc-list))
 413                                      (cond (add-to-simps
 414                                             (convert-relation-to-dot-simp tem1)))
 415                                      (cond ($fast_dotsimp ($rat_the_dot_simplifications)))
 416                                      ($check_overlaps
 417                                        up-to-degree
 418                                        add-to-simps nil lowest-degree)
 419                                      (return-from angela 'done)))))))))))))
 420   $dot_simplifications)
 421
 422 (defun rank-dimension-three (n)
 423   (cond ((oddp n)(div* (* (1+ n) (+ n 3)) 4))
 424         (t (div* (* (+ 2 n) (+ 2 n)) 4))))
 425
 426 ;  (cond ((null ok)
 427 ;       (cond ((y-or-n-p "Would you like to try again?")
 428 ;             (format t "~%Up to what degree would you like to test?
 429 ;  ~%Enter a number: ")
 430 ;                  (setq  up-to-degree (read))
 431 ;                  ($check_overlaps up-to-degree t)))))
 432
 433 (defvar *temp-pair* (list nil nil))
 434
 435 (defun ordered-pair-in-list (a b a-list-of-pairs)
 436   (let ((lis *temp-pair*))
 437     (setf (car lis) a (second lis) b)
 438     (gethash  lis a-list-of-pairs)))
 439
 440 (defmacro spsafe (&rest ll)
 441   (loop for u in ll
 442         do
 443         (cond ((get u 'special)(format t "~%Warning: ~A is a special" u))
 444               (t (format t "~%~A is not special" u)))))
 445
 446 ;;used &quote before but since never use old method of relations now.
 447
 448 (defun $set_up_dot_simplifications ( relations &optional check-thru-deg
 449                                     (ordering $order_function))
 450   (setq $dot_simplifications '((mlist simp)))
 451   (show $order_function)
 452   (cond (($listp (setq relations (meval* relations)))(setq relations (cdr relations))))
 453   (loop for relat in relations
 454         do
 455          (convert-relation-to-dot-simp relat ordering))
 456   (format t "~%The new dot simplifications are set up")
 457   (cond (check-thru-deg
 458          (displa $dot_simplifications)
 459          ($check_overlaps check-thru-deg t)))
 460   $dot_simplifications)
 461
 462 (defun $tes (relat) (convert-relation-to-dot-simp relat))
 463
 464
 465
 466
 467
 468 (defun convert-relation-to-dot-simp (relat &optional (ordering $order_function) &aux
 469                                      cof worst  )
 470   ordering
 471   (cond ($fast_dotsimp (setq relat ($dotsimp relat )))              ;;($vrat relat))))
 472         ((null $free_dot)(setq relat (meval* relat)))
 473         (t (setq relat ($dotsimp relat))))
 474   (cond
 475     (($zerop relat ) nil)
 476     ((numberp relat) (merror "Adding a constant relation  will result in the trivial algebra"))
 477     (t
 478      (cond
 479        ($new_fast_dotsimp
 480         (setq relat (num (cdr relat)))
 481         (setq worst (get (car relat) 'disrep)
 482               cof (third  relat))
 483         (setq relat  (nred relat cof))
 484        (setq relat (n-  worst relat))                ;;; should take cdddr etc. not subtractt
 485
 486         (setq relat (header-poly relat)))
 487 ;;;;    (setq relat (vsub* worst relat)))
 488        (t
 489         (multiple-value (worst cof )(find-worst-nc-monomial relat))
 490         (cond ((null $fast_dotsimp)
 491                (setq relat (div* relat cof))
 492                (setq relat
 493                      ($ratsimp (sub* worst relat))))
 494               (t
 495                (setq relat (vdiv* relat cof))
 496                (setq relat (vsub* ($vrat worst) relat))
 497 ;;;            (setq relat (check-rat-order relat))
 498                (setq relat (minimize-varlist relat))))))
 499      (format t "~%Adding  ")
 500      (displa (cons '(mlist simp) (list worst relat)))
 501      (format t "    to dot_simplifications")
 502      ($add_to_simps    (list  worst relat)))))           ;;;;(copy-tree relat))))))
 503
 504 (defvar $all_dotsimp_denoms nil)
 505
 506 (defun convert-relation-to-dot-simp (relat &optional (ordering $order_function) &aux
 507                                        tem tem1)
 508   ordering
 509   (cond ((or (null $fast_dotsimp) (null $free_dot))
 510          (fsignal "old code")))
 511   (setq relat ($dotsimp relat))
 512   (cond (($zerop relat) nil)
 513         (t (setq relat (cdr relat))
 514            (and (consp (num relat))
 515                 (consp (setq tem (third (num relat))))
 516                 (progn (format t "~%Denom : ")
 517                        (displa (setq tem1 (new-disrep tem)))
 518                        (if $all_dotsimp_denoms
 519                            (nconc $all_dotsimp_denoms (list tem1)))))
 520            ($add_to_simps (setq tem (make-simp (num relat))))
 521            (format t "~%Adding to simps: ")(displa (cons '(mlist) tem)))))
 522
 523
 524
 525 ;;       ;;;Cre form examples;;;
 526 ;;
 527 ;;   The global variables varlis and genvar are the third and fourth
 528 ;;items in the the (car expr) of expr if expr is in cre form.  The (cdr expr)
 529 ;;contains the actual information with its numerator being the car  and
 530 ;;the denominator the cdr.  When you add two expressions  expr1 and expr2 using add* the
 531 ;;current value of varlist and genvar are used.  If you put varlist equal to
 532 ;;nil then, ratrep* and its friend newvar will splice together the lists of
 533 ;;variables and set them up.  Prep1 does the actual calling of ratplus on the
 534 ;;cdr of the new expr1 and expr2.  They will have a common varlist which,
 535 ;;if the global varlist was nil, would be just the alphabetical splicing of
 536 ;;the two varlists.  Otherwise the global varlist would form the end of the
 537 ;;new varlist.  Since we usually just want to splice them, we introduce vadd*,etc.
 538 ;;which locally put the varlist to nil.  This method involves no ratdisrepping,
 539 ;;and so is quite fast.  The elements of genvar get reused again and again, unless
 540 ;;you specify genvar nil say, or somehow hide the old elements of genvar.
 541 ;;
 542 ;;   The elements of genvar sometimes get a 'disrep property, after disrepping.
 543 ;;This will not happen if you were not displaying intermediate results etc.
 544 ;;The elements of genvar get assigned the numerical values 1,2 ,3,4,...  as shown
 545 ;;below.
 546 ;;
 547 ;;   $gcd does not cause disrepping and can be used to find the common factor
 548 ;;of a numerator and denominator, so as to have the same effect as ratsimp
 549 ;;but without disrepping.
 550 ;;
 551 ;;   The order of varlist when splicing is alphabetical.  Thus if you want
 552 ;;to have zeta before a use %zeta or @zeta.  Otherwise you would have to add
 553 ;;a to varlist first, and that would give it precedence over b etc.  Note
 554 ;;that lists come after any atom in alphalessp order.
 555 ;;
 556 ;;   I add v to the name if we let varlist be nil.  Thus ($vrat '$a) will cause
 557 ;;the result to have only a varlist of ($a) .  The genvar will be the global
 558 ;;genvar.
 559 ;;
 560 ;;  The function minimize-varlist removes all elements from varlist and genvar
 561 ;;which don't appear in the expression.
 562 ;;
 563 ;;
 564 ;;
 565 ;;                                                     numerator     denominator:
 566 ;;((MRAT SIMP ($B) (#:G0558 #:G0557 #:G0556 #:G0555)) (#:G0558 1 1) . 1)  the cdr is 1
 567 ;;
 568 ;;             B
 569 ;;
 570 ;;             varlist       genvar                             numerator starts
 571 ;((MRAT SIMP ($A $B $C $D) (#:G0558 #:G0557 #:G0556 #:G0555))(#:G0556 1
 572 ;;                                                                     (#:G0557 1
 573 ;;                                                                              1
 574 ;;                                                                              0
 575 ;;                                                                              (#:G0558 1 1))
 576 ;;                                                                      0
 577 ;;                                                                      2)
 578 ;;                                                             #:G0555  <--denominator starts
 579 ;;                                                             2           and is cdr
 580 ;;                                                             1
 581 ;;                                                             0
 582 ;;                                                             1)
 583 ;;
 584 ;;          (B + A) C + 2
 585 ;;          -------------
 586 ;;              2
 587 ;;             D  + 1
 588 ;;
 589 ;;
 590 ;;   G0558 disreps $B value 1
 591 ;;   G0557 disreps $B value 2
 592 ;;   G0556 disreps $C value 3
 593 ;;   G0555 disreps $D value 4
 594 ;;
 595 ;;
 596 ;;   If one is going to do a lot of computations with entries of a matrix
 597 ;;say which all have more or less the same variables occurring one could
 598 ;;prepare them all to have the same varlist and genvar as follows.  Go
 599 ;;through them collecting the total varlist and sorting it into the right
 600 ;;order.  Call it new-varlist. You could at the same time collect
 601 ;;genvar's.  Now let varlist be the new collected one and genvar be the
 602 ;;collected one.  Then do (setq expr (cons (standard) (prep1 expr))) for
 603 ;;each entry expr, where standard is (list 'mrat simp new-varlist
 604 ;;new-genvar).  Then outside everything put a (let ((genvar new-genvar))
 605 ;;
 606 ;;   Prep1 does not disrep expr if the global varlist contains all
 607 ;;(varlist expr) in the right order, and if the global genvar
 608 ;;is as long as the varlist.  Note that rattimes will not automatically
 609 ;;simplifying such things as %i^2, so one may still want to just use vmul*.
 610 ;;Prep1 does not alter (cdr expr) if the varlist and genvar of expr are equal to
 611 ;;the global genvar and varlist.
 612 ;
 613 ;(defun vadd* (&rest llist)
 614 ;  (let ((varlist))
 615 ;    (apply 'add*  llist)))
 616 ;(defun vsub* (a b )
 617 ;  (setq b (vmul* -1 b))
 618 ;  (vadd* a b))
 619 ;(defun vsub* (a b )
 620 ;  (vadd* a (vminus* b)))
 621 ;
 622
 623 ;(defun vmul* (&rest a-list))
 624 ;
 625
 626
 627
 628 ;;;the following convert all to rat form using $new_rat unless given a affine-polynomialp or $ratp
 629
 630
 631
 632 ;(setq $new_fast_dotsimp t)
 633 (defun vadd* (&rest a-list)
 634   (cond ((null a-list) 0)
 635         ((eql (length a-list) 1) (car a-list))
 636         (t (header-poly (n+ (car a-list) (apply 'vadd* (cdr a-list)))))))
 637 (defun vmul* (&rest a-list)
 638   (cond ((null a-list) 1)
 639         ((eql (length a-list) 1) (car a-list))
 640         (t (header-poly (n* (car a-list) (apply 'vadd* (cdr a-list)))))))
 641 (defun vsub* (a b) (header-poly (n- a b)))
 642 (defun vdiv* (a b) (header-poly (nred  a b)))
 643 (defun vminus* (a)(header-poly (n- 0 a)))
 644 (defun $vrat (&rest a-list)
 645   (apply '$new_rat a-list))
 646
 647
 648
 649
 650 ;;second
 651 ;(setq $new_fast_dotsimp nil)
 652
 653
 654 (defun $vrat (&rest a-list)
 655   (let ((varlist))
 656     (apply '$rat  a-list)))
 657 (defun vadd* (&rest a-list)
 658   (let ((varlist ))
 659     (cond ((< (length a-list) 3)(apply 'add* a-list))
 660           (t (vadd* (car a-list) (apply 'vadd* (cdr a-list)))))))
 661 (defun vmul* (&rest a-list)
 662   (let ((varlist ))
 663     (cond ((< (length a-list) 3)(apply 'mul* a-list))
 664           (t (vmul* (car a-list) (apply  'vmul* (cdr a-list)))))))
 665 (defun vsub* (a b)
 666   (let ((varlist))
 667     (sub* a b)))
 668 (defun vdiv* (a b)
 669   "This needs to be fixed.  It still allows ratdisrep"
 670   (let ((varlist))
 671    (simplifya  (list '(mquotient) a b) nil)))
 672 (defun vminus* (a)
 673   (let ((varlist))
 674     (simplifya (list '(mminus) a) nil)))
 675
 676 ;(defun $add_relation_to_dot_simplifications( relat &optional (ordering $order_function))
 677 ;  (convert-relation-to-dot-simp relat ordering))
 678 (defun $add_to_simps (a-list &optional(ordering $order_function) )
 679   (cond (($listp a-list)(setq a-list (cdr a-list))))
 680   (cond ((and (zerop $expop)(null $fast_dotsimp))
 681          (setf (second a-list) ($multthru (second a-list)))))
 682   (loop for u in $dot_simplifications
 683         for i from 0
 684         when (and (oddp i) (not (funcall ordering  (car a-list) u)))
 685         do (setq $dot_simplifications (append (subseq $dot_simplifications 0 i)
 686                                                a-list
 687                                        (nthcdr i $dot_simplifications)))
 688         (return 'done)
 689         finally (setq $dot_simplifications (append $dot_simplifications a-list)))
 690   ($simplify_dot_simplifications ($nc_degree (car a-list))))
 691
 692 (defmacro $with_no_simp (form)
 693   (let (($dot_simplifications nil))
 694     `(progn  ' ,(meval* form))))
 695
 696 (defun $inverse_modulo (n modulo-p)
 697   "The inverse of any number n modulo-p.  ok if n negative, but n=0 gives 0
 698 and modulo-p not prime gives false answer"
 699   (mod (expt  n (- modulo-p 2)) modulo-p))
 700
 701 (defun $nu (i llist &aux
 702   (p1  (length llist)))
 703   (let ((fact ($inverse_modulo i p1)))
 704   (loop for ii from 1 below p1
 705         collecting
 706         (nth (mod (* ii fact) p1) llist) into tem
 707         do (show tem ii)
 708         finally (return (cons '(mlist simp) tem)))))
 709 (defun $phi (l )
 710   (check-arg l '$listp "macsyma list")
 711   (setq l (cdr l))
 712   (let ((p0 (length l)))
 713     (show p0)
 714     (add*  (mul* p0 (loop for i from 1 to (1- p0)
 715                           collecting (mul* (nth i l) (power '$sig (mod (- i) p0)))
 716                           into tem
 717                           finally (return (meval* (cons '(mplus) tem))))
 718                  (car l)
 719                  (mul* (- p0) (loop for i from 1 to (1- p0)
 720                                         collecting (nth i l) into tem
 721                                         finally (return (meval* (cons '(mplus) tem)))))))))
 722
 723
 724 (defvar $prime_order nil)
 725
 726 (defun $sig (n &optional (p0 $prime_order))
 727   (power '$sig (mod n p0)))
 728
 729 (defun $nu_poly (i poly)
 730   ($ratsimp (subst ($sig i $prime_order) '$sig poly)))
 731
 732 (defun $scalar_concat (&rest l)
 733   (let ((tem (apply '$concat l)))
 734     (putprop tem '(nil $scalar t) 'mprops)
 735     tem))
 736
 737 (defun $rel (i j)
 738   (let ((ii (max i j)) ans (jj (min i j)))
 739     (cond ((not (equal i j))
 740            (setq ans (list (ncmul* ($concat '$e ii) ($concat '$e jj))
 741                            (sub* ($scalar_concat '$a jj ii)
 742                                  (ncmul* ($concat '$e jj) ($concat '$e ii))))))
 743           (t (setq ans (list (ncmul* ($concat '$e ii) ($concat '$e jj))
 744                            ($scalar_concat '$a jj ii)))))
 745     (cons '(mlist simp) ans)))
 746
 747 (defun $clifford_dot_simplifications (n )
 748   (loop for i from 1 to n
 749         appending
 750         (loop for j from 1 to i
 751               appending (cdr ($rel i j)))
 752         into tem
 753         finally (return (cons '(mlist simp) tem))))
 754
 755
 756 (defun $replacements (&aux (tem (cdr $dot_simplifications)))
 757   (loop while tem
 758         collecting (car tem) into a-list
 759         do (setq tem (cddr tem))
 760         finally (return (cons '(mlist simp) a-list))))
 761
 762
 763 (defvar $list_of_zeroes nil)
 764
 765 (defun $earlier_mono (monom &optional (variables $current_variables) &aux answer)
 766   (setq answer ($mono variables ($nc_degree monom)))
 767   (loop for v in (cdr answer)
 768          when (funcall $order_function v monom)
 769          collecting v into tem
 770          finally (return (cons '(mlist simp) (sort tem $order_function)))))
 771
 772 (defun parse-string (a-string)
 773   (cond ((null (or (search "$" a-string :test #'char-equal) (search ";" a-string :test #'char-equal)))
 774          (setq a-string  (string-append a-string "$"))))
 775   (with-input-from-string (stream a-string)
 776     (let ((*standard-input* stream) answer)
 777          (setq answer (caddr (mread)))
 778            answer)))
 779 (defmacro   mac (a-string )
 780   "Parses a Macsyma string but does not call meval*"
 781       (list 'quote (parse-string  a-string)))
 782
 783 (defmacro   mac-eval-after (a-string )
 784   "Parses a Macsyma string but evaluates it at run time"
 785       (list 'quote `(meval* ,(parse-string  a-string))))
 786 (defmacro   mac-eval (a-string )
 787   "Reads a macsyma string and evaluates it at compile time returning the quoted form"
 788       (list 'quote (meval* (parse-string  a-string))))
 789
 790 (defun $hh (a)
 791   (cond ((atom a) nil)
 792         (t (equal (caar a) 'mtimes))))
 793
 794 (defun $monomial_dimensions (n &optional (order-weight nil)  &aux j  the-count)
 795   (cond ($current_variables
 796          (format t "~%The current variables are")
 797          (displa $current_variables)
 798          (cond ((y-or-n-p "Use them?"))
 799                 (t (setq $current_variables nil)))))
 800   (cond ($current_variables nil)
 801         (t (format t "~%Enter Macsyma list of variables to use:") (setq $current_variables
 802                                                       (caddr (mread-raw *standard-input*)))))
 803   (cond (order-weight
 804          (loop for j from 1 to n
 805                do (setq the-count 0)
 806                (loop for i from 1 to n
 807                      do
 808                      (loop for v in (cdr ($mono $current_variables i))
 809                            when (eq ($nc_degree v :order-weight) j)
 810                            do (incf the-count)))
 811                (format t "~%There are ~A monomials of order weight ~A and of degree less than ~A"
 812                        the-count j n)
 813                summing the-count into tem
 814                finally (format t "~%The total number of monomials of weight and degree less than ~A is ~A." n (1+ tem)))))
 815
 816   (loop for i from 1 to n
 817         do
 818         (format t "~%There are ~A independent monomials in degree ~A."
 819                 (setq j (length  (cdr ($mono $current_variables i)))) i)
 820                 collecting (list '(mlist) i j ) into tem
 821         summing j into the-sum
 822   finally (format t "~%The sum of the dimensions from dimension 0 through dimension ~A is ~A"
 823                   n (+ the-sum 1)) (return (cons '(mlist simp) tem))))
 824
 825
 826 (defun $mono_weighted (n &aux the-list)
 827   (loop for i from 1 to n
 828         do
 829         (loop for v in (cdr ($mono $current_variables i))
 830               when (eq ($nc_degree v :order-weight) n)
 831               do (setq the-list (cons v the-list)))
 832         finally (return (cons '(mlist) the-list))))
 833
 834 (defmacro subscript (d i)
 835   `((,d  :array) ,i))
 836 (defvar $last_equations nil)
 837 (defvar $last_solutions nil)
 838 (defun $find_relations_among (a-list &optional (dotsimp-first nil) (term-names nil)&aux eqns gen-sum answers)
 839   (check-arg a-list '$listp "macsyma list")
 840   (cond (dotsimp-first (setq a-list
 841                              (loop for v in (cdr a-list)
 842                              collecting ($dotsimp v) into tem
 843                              finally (return (cons '(mlist) tem))))))
 844   (setq gen-sum ($general_sum  a-list $aaaa))
 845   (cond ((null dotsimp-first) (setq gen-sum ($dotsimp gen-sum))))
 846   (setq $last_equations (setq eqns ($nc_coefficients ($ratsimp `((mlist) , gen-sum)))))
 847   (setq answers (apply '$fast_linsolve (cdr eqns)))
 848   (setq $last_solutions answers)
 849   (cond ((null term-names)(setq term-names (loop for i from 1 to (length (cdr a-list))
 850                                        collecting ($concat '$term i) into tem
 851                                        finally (return (cons '(mlist) tem))))))
 852   ($sublis answers ($general_sum term-names
 853                                  $aaaa)))
 854
 855
 856
 857
 858 (deff $te #'$power_series_monomial_alphalessp)
 859
 860 (defun $check_skew_relations (&optional(variables $current_variables) &aux tem answer)
 861   (loop for  u in (cdr variables)
 862         unless ($must_replacep u)
 863         do
 864         (loop for v in (cdr variables)
 865               when (and (funcall $order_function u v) (not (equal u v)))
 866               do
 867               (cond (($must_replacep (setq tem (ncmul* v u))) nil)
 868                     (t (format t "~%no replacement for ~A" tem)
 869                        (setq answer (cons tem answer)))))
 870         finally (cond ((null answer)(return t))
 871                       (t (cons '(mlist) answer)))))
 872
 873 (defun get-alt (item a-list)
 874   (loop while a-list
 875         when (equal item (car a-list))
 876         do(return (second a-list))
 877         else do (setq a-list (cddr a-list))))
 878
 879 (defun $collect_skew_relations (&optional(variables $current_variables) &aux tem  tem2)
 880   (loop for  u in (cdr variables)               ;a
 881         unless ($must_replacep u)
 882         appending
 883         (loop for v in (cdr variables)
 884               when (and (funcall $order_function u v) (not (equal u v))
 885                    (setq tem2 (get-alt (setq tem (ncmul* v u)) (cdr $dot_simplifications))))
 886               collecting tem
 887               and
 888               collecting tem2 )
 889         into a-list
 890         finally (return (cons '(mlist simp) a-list))))
 891
 892 (defun $relations_from_dot_simps(&optional (dot-simps $dot_simplifications))
 893   (let ((a-list (cdr dot-simps)))
 894         (loop
 895               while a-list
 896               collecting  (sub* (first a-list) ($totaldisrep (second a-list))) into tem
 897               do
 898               (setq a-list (cddr a-list))
 899               finally (return (cons '(mlist) tem)))))
 900 (defun repl-deg ()
 901   (sort (mapcar '$nc_degree (cdr ($replacements))) 'alphalessp))
 902
 903 (defun $pb (n)
 904   (mfuncall '$playback (list '(mlist) n 1000)))
 905
 906 (defvar $centrals_so_far nil)
 907
 908 (defun $find_central_thru_deg (begin to-n &aux tem)
 909   (displa $current_variables)
 910   (setq $centrals_so_far nil)
 911   (loop for i from begin  to to-n
 912         do
 913         ($check_overlaps (1+ i) t nil i)
 914         (setq $centrals_so_far
 915               (cons (list '(mlist) i
 916                           (setq tem (mfuncall '$central_elements $current_variables i)))
 917                     (cdr $centrals_so_far)))
 918         (maclist $centrals_so_far)
 919         (format t "~%Central in degree ~A " i)
 920         (displa tem))
 921   $centrals_so_far)
 922
 923 (defmacro maclist (x)
 924   `(setq ,x (cons '(mlist) ,x)))
 925
 926 (defmacro rat-variable-info (cre-form)
 927   `(car ,cre-form))
 928
 929 (defun fast-scalarp (x)
 930   (cond ((atom x)
 931          (cond ((numberp x) t)
 932                ((symbolp x) (mget x '$scalar))
 933                (t ($scalarp x))))
 934         (t
 935          (case (caar x)
 936            (mnctimes nil)
 937            (mrat
 938             (let ((gen-vars (fourth (rat-variable-info x))))
 939               (loop for vv in gen-vars
 940                      when (not (fast-scalarp (get vv 'disrep)))
 941                      do
 942                      (cond ((appears-in  (cdr x) vv) (return nil)))
 943                      finally (return t))))
 944            (otherwise ($scalarp x))))))
 945
 946 (defun my-ratcoeff (expr monom &optional (deg 1) &aux answer )
 947  (setq answer
 948        (catch 'bad-variable-order
 949     (cond ((atom expr)(cond ((eq expr monom) 1)
 950                           (t 0)))
 951         ((and (equal( caar expr) 'mrat) (eq (second (car expr)) 'simp))
 952          (my-ratcoeff1 expr monom deg))
 953         (t ($ratcoef expr monom)))))
 954  (cond ((eq answer 'bad-order)(format t "~%Variable ~A is out of order so using $ratcoef."
 955                                       monom)
 956         ($ratcoef expr monom deg))
 957        (t answer)))
 958
 959 (defun my-ratcoeff1 (expr monom &optional (deg 1))
 960   (let ((var-list (varlist expr))
 961         (gen-vars (genvar expr)))
 962     (loop for v in var-list
 963           for w in gen-vars
 964           when (equal v monom)                  ;
 965           do
 966           (cond ((not (appears-in (num (cdr expr)) w)) (return 0))
 967                 (t
 968                  (return (cons (rat-variable-info expr) (cons (my-ratcoeff2 (num (cdr expr))
 969                                                                             w
 970                                                                             deg)
 971                                                               (denom (cdr expr)))))))
 972           finally (return 0))))
 973 ;
 974 ;(defun my-ratcoeff2(rat-expr gen-var &optional (deg 1))
 975 ;  (cond ((atom rat-expr) (rzero))
 976 ;
 977 ;       (t (let ((rat-num rat-expr))
 978 ;                                               ; (gshow rat-num)
 979 ;            (cond ((eql (car rat-num) gen-var)
 980 ;                   (setq rat-num (cdr rat-num))
 981 ;                   (loop while rat-num
 982 ;                         when (eql (car rat-num) deg)
 983 ;                         do
 984 ;                         (return  (second rat-num) )
 985 ;                         when (eql (car rat-num) 0)
 986 ;                         do
 987 ;                         (return (my-ratcoeff2 (second rat-expr) gen-var))
 988 ;                         else
 989 ;
 990 ;                         do
 991 ;                         (setq rat-num (cddr rat-num))))
 992 ;                  (t (setq rat-num (cdr rat-num))
 993 ;;                    (break t)
 994 ;                   (loop while rat-num
 995 ;                           when (eql (car rat-num) 0)
 996 ;                           do
 997 ;                           (return (my-ratcoeff2 (second rat-num) gen-var))
 998 ;                           else
 999 ;
1000 ;                           do (setq rat-num (cddr rat-num))
1001 ;;                          (gshow rat-num)
1002 ;                           )))))))
1003
1004
1005 (defun my-ratcoeff2(rat-expr gen-var &optional (deg 1))
1006   (cond ((atom rat-expr) (rzero))
1007
1008         (t (let ((rat-num rat-expr))
1009                                                 ; (gshow rat-num)
1010              (cond ((eq (car rat-num) gen-var)
1011                     (setq rat-num (cdr rat-num))
1012                     (loop while rat-num
1013                           when (eql (car rat-num) deg)
1014                           do
1015                           (return  (second rat-num) )
1016                           when (eql (car rat-num) 0)
1017                           do
1018                           (return (my-ratcoeff2 (second rat-expr) gen-var))
1019                           else
1020
1021                           do
1022                           (setq rat-num (cddr rat-num))))
1023                    (t (setq rat-num (cdr rat-num))
1024 ;                     (break t)
1025                     (loop while rat-num
1026                             when (eql (car rat-num) 0)
1027                             do
1028                             (return (my-ratcoeff2 (second rat-num) gen-var))
1029                             else
1030                             do
1031                             (cond ((appears-in (cadr rat-num) gen-var)
1032                                    (throw 'bad-variable-order 'bad-order))
1033                                   (t(setq rat-num (cddr rat-num))))
1034 ;                           (gshow rat-num)
1035                             )))))))
1036 ;
1037 ;(defun $numerator (expr)
1038 ;;  (check-arg expr '$ratp "macsyma cre form")
1039 ;
1040 ;  (cond
1041 ;    ((numberp expr) expr)
1042 ;    ((mbagp expr)(cons (car expr) (mapcar '$numerator (cdr expr))))
1043 ;    (($ratp expr) (cons (car expr)(cons (cadr expr) 1)))
1044 ;        ((equal (caar expr) 'mtimes)
1045 ;                (loop for v in (cdr expr)
1046 ;                      unless  (and (not (atom v))
1047 ;                                 (equal (caar v) 'mexpt)
1048 ;;                                (show (third v))
1049 ;                                  (<  (third v) 0))
1050 ;                      collecting v into tem
1051 ;                      finally (return (cons '(mtimes) tem))))
1052 ;       (t ($numerator ($vrat expr)))))
1053
1054
1055 (defun $denominator (expr &aux tem1)
1056   (cond (($ratp expr)
1057          (cons (car expr)(cons (cddr expr) 1)))
1058         ((equal (caar expr) 'mtimes)
1059                  (loop for v in (cdr expr)
1060                        when  (and (not (atom v))
1061                                   (equal (caar v) 'mexpt)
1062 ;                                 (show (third v))
1063                                    (<  (third v) 0))
1064                        collecting (progn
1065                                     (setq tem1 (copy-list v))
1066                                     (setf (third tem1) (- (third v)))
1067                                     tem1)
1068                                   into tem
1069                        finally (return (simplifya (cons '(mtimes) tem) nil))))
1070         (t ($denominator ($vrat expr)))))
1071 (defun $te (x y)(my-ratcoeff x y))
1072 (defun variables-same-divide (f g)
1073   (cons (car f) (car (ratdivide (cdr f) (cdr g)))))
1074 (defun variables-same-divide (f g)
1075   (vdiv* f g))
1076 (defun simplify-rational-quotient (expr &aux the-gcd )
1077   "Divides top and bottom of a cre expression by their gcd"
1078   (check-arg expr '$ratp "Cre form")
1079   (let ((the-num ($numerator expr))
1080         (the-denom ($denominator expr)))
1081     (cond ((equal (cdr the-denom) (cons 1 1))
1082            expr)
1083           (t (setq the-gcd ($gcd the-num the-denom))
1084              (setq the-num (variables-same-divide the-num the-gcd))
1085              (setq the-denom (variables-same-divide the-denom the-gcd))
1086              (append (list (car expr) (num (cdr the-num))) (num (cdr the-denom)))))))
1087
1088
1089
1090 (defun two (x) x 2)
1091
1092 (defmvar $skew3_conditions)
1093 (defvar $bbbb)
1094 (defvar $cccc)
1095
1096 (defun $check_skew3_conditions ( relations  &aux monoms3 ($expop 100))
1097
1098   (setq relations
1099         (list '(mlist)  (div*  (nth 1 relations)
1100                               ($coeff (nth 1 relations) (ncmul* '$y '$x '$x)))
1101               (div*  (nth 2 relations)
1102                               ($coeff (nth 2 relations) (ncmul* '$y '$y '$x)))))
1103   (displa relations)
1104   (let (($dot_simplifications nil)(result t) coeff-values answer )
1105     (setq monoms3 ($mono '((mlist) $x $y) 3))
1106     (setq coeff-values
1107       (loop for relat in (cdr relations)
1108           for coeff-names in (list $bbbb $cccc)
1109           appending
1110           (loop for u in (cdr monoms3)
1111                 for j from 1
1112                 collecting `((mequal) ,(nth j coeff-names) ,($coeff  relat u)))
1113           into tem
1114           finally (return (cons '(mlist) tem))))
1115     (setq answer ($sublis coeff-values $skew3_conditions))
1116     (loop for v in (cdr answer)
1117           when (not ($zerop (sub*  (nth 1 v) (nth 2 v))))
1118           do (format t "~%Inconsistent condition ")
1119           (displa v)
1120           (setq result nil))
1121     (cond ((null result)(format t "~%Trying to interchange f and g...")
1122
1123            ($check_skew3_conditions (list '(mlist) (reverse (cdr relations)))))
1124           (t result))))
1125 (defun $list_variables (expr &rest key-strings &aux string-v tem)
1126   (cond ((null key-strings) ($listofvars expr))
1127         (t
1128          (loop for v in key-strings
1129                collecting (string-trim "&" v) into tem
1130                finally (setq key-strings tem))
1131          (loop for v in (cdr (setq tem ($listofvars expr)))
1132                do (setq string-v (string v))
1133                when (not (loop for key in key-strings
1134                                when  (search key string-v :test #'char-equal)
1135                                do (return t)))
1136                do (setq tem (delete v tem :test #'equal))
1137                finally (return tem)))))
1138 (defun $socle (variables deg &aux  f unknowns eqns tem1 parameters answer )
1139   (setq variables (cons '(mlist) (sort (cdr variables) $order_function)))
1140   (let* ((monoms ($mono variables deg))
1141         (monoms-higher ($mono variables (1+ deg))))
1142
1143     (setq f ($general_sum monoms $aaaa))
1144     (loop for v in (cdr variables)
1145           do
1146           (setq unknowns ($list_variables f "aa" "par"))
1147           (setq parameters (loop for vv in (cdr unknowns)
1148                                  when (search "par" (string vv) :test #'char-equal)
1149                                  collecting vv))
1150           (setq tem1 (cdr $aaaa))
1151           (loop for vv in parameters
1152                 do (loop while tem1
1153                          when (not (member (car tem1) unknowns :test #'eq))
1154                          do (setq f (subst (car tem1) vv f))
1155                          (setq unknowns (subst (car tem1) vv unknowns))
1156                          (format t "~%Replacing ~A by ~A in f." vv (car tem1))
1157                          (setq tem1 (cdr tem1)) (return 'done)
1158                          do (setq tem1 (cdr tem1))))
1159           (setq tem1 (list '(mlist) (ncmul* f v) (ncmul* v f)))
1160           (setq eqns ($extract_linear_equations tem1 monoms-higher))
1161           (show unknowns)
1162           (setq answer ($fast_linsolve eqns unknowns))
1163           (setq f (meval* ($sublis answer f)))
1164           (displa f)
1165           finally (return f))))
1166
1167
1168
1169 (defvar $commutators '((mlist)))
1170
1171 (defvar $centralizers  '((mlist)))
1172
1173 (defun $fast_normal_elements (deg &key (variables $current_variables) up_to_deg
1174                                   auto &aux  f unknowns eqns tem1 parameters answer tem )
1175   "Returns elements F in degree DEG in variables, which satisfy u.F=F.phi(u)
1176   where phi is the auto  which defaults to the identity and can be specified by
1177   the list of images of VARIABLES"
1178      (cond ((null auto) (setq auto variables)))
1179      (setq tem (sort (pairlis (cdr variables) (cdr auto)) #'$order_function :key #'car))
1180      (setq variables (cons '(mlist) (mapcar 'car tem)) auto (cons'(mlist)
1181                                                               (mapcar 'cdr tem)))
1182   (setq $commutators nil $centralizers nil)
1183   (let* ((monoms (cond (up_to_deg (cons '(mlist) (loop for i from 1 to deg appending (cdr  ($mono variables i)))))
1184                        (t ($mono variables deg))))
1185 ;       (monoms-higher ($mono variables (1+ deg)))
1186         )
1187     (setq f ($general_sum monoms $aaaa))
1188     (loop for v in (cdr variables)
1189           for im in (cdr auto)
1190           do
1191           (setq unknowns ($list_variables f "aa" "par"))
1192           (setq parameters (loop for vv in (cdr unknowns)
1193                                  when (search "par" (string vv) :test #'char-equal)
1194                                  collecting vv))
1195           (setq tem1 (cdr $aaaa))
1196           (loop for vv in parameters
1197                 do (loop while tem1
1198                          when (not (member (car tem1) unknowns :test #'eq))
1199                          do (setq f (subst (car tem1) vv f))
1200                          (setq unknowns (subst (car tem1) vv unknowns))
1201                          (format t "~%Replacing ~A by ~A in f." vv (car tem1))
1202                          (setq tem1 (cdr tem1)) (return 'done)
1203                          do (setq tem1 (cdr tem1))))
1204           (setq tem1 (list '(mlist)($totaldisrep ($dotsimp (sub* (ncmul* v f) (ncmul* f  im))))))
1205           (setq $commutators (append `((mlist) ,v ,(second tem1)) (cdr $commutators)))
1206 ;         (setq eqns ($extract_linear_equations tem1 monoms-higher))
1207           (setq eqns ($extract_linear_equations tem1 ($list_nc_monomials tem1)))
1208           (show unknowns)
1209           (setq answer ($fast_linsolve eqns unknowns))
1210           (setq f ($ratsimp ($sublis answer f)))
1211           (setq $centralizers (append `((mlist) ,v ,f) (cdr $centralizers)))
1212           (displa f)
1213           finally (return ($separate_parameters f)))))
1214 (defun $fast_central_elements (variables deg  &aux  f unknowns eqns tem1 parameters answer )
1215   (setq variables (cons '(mlist) (sort (cdr variables) $order_function)))
1216   (setq $commutators nil $centralizers nil)
1217   (let* ((monoms ($mono variables deg))
1218 ;       (monoms-higher ($mono variables (1+ deg)))
1219         )
1220
1221     (setq f ($general_sum monoms $aaaa))
1222     (loop for v in (cdr variables)
1223           do
1224           (setq unknowns ($list_variables f "aa" "par"))
1225           (setq parameters (loop for vv in (cdr unknowns)
1226                                  when (search "par" (string vv) :test #'char-equal)
1227                                  collecting vv))
1228           (setq tem1 (cdr $aaaa))
1229           (loop for vv in parameters
1230                 do (loop while tem1
1231                          when (not (member (car tem1) unknowns :test #'eq))
1232                          do (setq f (subst (car tem1) vv f))
1233                          (setq unknowns (subst (car tem1) vv unknowns))
1234                          (format t "~%Replacing ~A by ~A in f." vv (car tem1))
1235                          (setq tem1 (cdr tem1)) (return 'done)
1236                          do (setq tem1 (cdr tem1))))
1237           (setq tem1 (list '(mlist) ($com f v)))
1238           (setq $commutators (append `((mlist) ,v ,(second tem1)) (cdr $commutators)))
1239 ;         (setq eqns ($extract_linear_equations tem1 monoms-higher))
1240           (setq eqns ($extract_linear_equations tem1 ($list_nc_monomials tem1)))
1241           (show unknowns)
1242           (setq answer ($fast_linsolve eqns unknowns))
1243           (setq f ($ratsimp ($sublis answer f)))
1244           (setq $centralizers (append `((mlist) ,v ,f) (cdr $centralizers)))
1245           (displa f)
1246           finally (return ($separate_parameters f)))))
1247
1248 (defun $central_elements_in_degrees (&rest degrees &aux $centrals_so_far)
1249   (loop for i in degrees do
1250     ($check_overlaps (1+ i) t nil)
1251     (push ($fast_central_elements $current_variables i)
1252                             $centrals_so_far)
1253                       (mapcar 'displa $centrals_so_far))
1254   (cons '(mlist) $centrals_so_far))
1255
1256
1257 (defun $pbi (n &aux tem)
1258   (let ((linel (- linel 10)))
1259     (loop for i from n below $linenum
1260           do  (format t "~% ~4A: ~A" (string-trim "$" (string (setq tem ($concat '$c i))))
1261                       (string-grind (meval* tem))))))
1262
1263 (defvar *pbi-string* (make-array '(64) :element-type 'standard-char
1264                                  :fill-pointer 0 :adjustable t))
1265
1266 (defun new-concat (a &rest b &aux me)
1267   (loop for v in b collecting
1268         (format nil "~A" v) into tem
1269         finally (setq me (apply 'string-append tem))
1270         (return (intern (format nil "~A~A" a me) 'maxima))))
1271
1272
1273 (defun $list_sublis (a-list b-list expr)
1274   (cond (($listp a-list) (setq a-list (cdr a-list))))
1275   (cond ( ($listp b-list )(setq b-list (cdr b-list))))
1276   (loop for v in a-list
1277         for w in b-list
1278         collecting (cons v w) into tem
1279         finally (return (sublis tem expr))))
1280
1281
1282 (defun string-grind (x &key ($display2d) stream &aux answ)
1283   (setq answ(with-output-to-string (st)
1284               (cond ((null $display2d)  (mgrind x st))
1285                     (t  (let ((*standard-output* st))
1286                           (displa x))))))
1287   (setq answ (string-trim '(#\newline #\space  #\$) answ))
1288   (cond (stream  (format stream "~A" answ)) ;
1289         (t answ)))
1290
1291 (defun macsyma-typep (x)
1292   "acceptable to displa"
1293   (or (numberp x)
1294       (and (consp x)(consp (car x))(get (caar x) 'dimension))))
1295
1296 ;If foo is (#:X2 1 1 0 (#:X1 1 1)) for instance, then
1297 ;(format t "The functions  ~VQ  are the inverses" foo 'fsh)
1298 ;(format t "~%The functions  ~VQ  are  inverses" (st-rat #$[x1+1,1/ (x1+1)]$)  'fsh)
1299 ;The functions  X1+1 and 1/(X1+1)  are  inverses
1300
1301 (defun fsh ( x &optional (stream *standard-output*))
1302   (let ($display2d)
1303     (cond
1304       ((macsyma-typep x)(string-grind x :stream stream))
1305       ((or (affine-polynomialp x) (rational-functionp x))
1306        (string-grind (header-poly x) :stream stream))
1307       ((and ( listp x)(or (affine-polynomialp (car x)) (rational-functionp (car x))))
1308        (case (length x)
1309          (1 (fsh x stream))
1310          (2 (fsh (first x) stream) (format stream " and ") (fsh (second x) stream))
1311          (t
1312            (loop for v in  x
1313                  for i from 1 below (length x)
1314                  do (fsh v stream) (format stream ",")
1315                  finally (format stream " and ") (fsh (car (last x)) stream)))))
1316       ((and (listp x) (macsyma-typep (car x)))
1317        (format stream "(")
1318        (loop for v in x do (fsh v) (format t " ")
1319              finally (format t ")")))
1320       (t (des x)))))