share/affine/polysmp.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 get-altq (item a-list)
  11   (loop for v on a-list by #'cddr
  12         when (eq item  (car v))
  13         do (return (second v))))
  14
  15 (defun mdegree (monom var)
  16   "Degree of a monomial in var Returns nil if var not in monom"
  17   (do ((v monom ))
  18       (nil)
  19     (cond ((numberp v) (return nil)))
  20     (cond ((eq (car v) var)(return (second v))))
  21     (setq v (caddr v))))
  22
  23 (defun part-above-degree (poly  monomial-degs &aux tem)
  24   "The argument is  a Cre polynomial and the
  25   monomial degs is an alternating list of the variable and the degree
  26   Note these must be ordered in increasing degree for the degree sequence"
  27
  28   (cond ((null monomial-degs)poly)
  29         ((numberp monomial-degs) poly)
  30         ((numberp poly)(cond ((loop for w in (cdr monomial-degs) by #'cddr
  31                                     when (not (eql w 0))
  32                                     do(return t)) nil)
  33                              (t poly)))
  34         (t (let ((deg))
  35              (cond ((eq (car monomial-degs)(car poly))
  36                     (setq deg (1- (second monomial-degs)))
  37                     (setq monomial-degs (cddr monomial-degs)))
  38                    (t  (setq deg (get-altq (car poly) monomial-degs))
  39                        (cond (deg (setq deg (1- deg)))
  40                              (t (setq deg -1)))))
  41              (loop for v  on (cdr poly) by #'cddr
  42                    while (> (car v)  deg)
  43                    when (setq tem (part-above-degree (second v)  monomial-degs))
  44                    collecting (car v) into the-cof
  45                    and
  46                    collecting tem into the-cof
  47                    finally (cond (the-cof (setq the-cof  (cons (car poly) the-cof))
  48                                           (return the-cof))
  49                                  (t nil)))))))
  50
  51 (defun part-above (poly monom)
  52  "x^3*y+x^2*u, x^2 ==> x*y+u"
  53   (cond ((part-above1 poly monom))
  54         (t 0)))
  55
  56 (defun part-above1 (poly monom &aux tem (true-deg 0))
  57   (cond ((numberp monom) poly)
  58         ((numberp poly) nil)
  59         (t (let ((deg))
  60              (cond ((eq (p-var monom) (p-var poly))
  61                     (setq deg (1- (second monom)))
  62                     (setq true-deg (second monom))
  63                     (setq monom (caddr monom)))
  64                    (t (setq deg (mdegree monom (car poly)))
  65                       (cond (deg (setq deg (1- deg)))
  66                             (t (setq deg -1)))))
  67              (loop for v on (cdr poly) by #'cddr
  68                     while (> (car v) deg)
  69                     when (setq tem (part-above1 (second v) monom))
  70                     collecting (- (car v) true-deg)
  71                     into the-cof
  72                     and
  73                     collecting tem into the-cof
  74                     finally (cond (the-cof (setq the-cof (cons (car poly) the-cof))
  75                                            (return the-cof))
  76                                   (t nil)))))))
  77 ;
  78 ;(defmacro push-new (item llist)
  79 ;       `(cond ((not (member  ,item ,llist :test #'eq) )
  80 ;              (push ,item ,llist))))
  81
  82
  83 ;;see polyd
  84 ;(defun list-variables (u &aux a-list)
  85 ;  (declare (special a-list))
  86 ;  (list-variables2 u)
  87 ;  a-list)
  88 ;
  89 ;(defun list-variables2 (expr)
  90 ;  (cond ((atom expr) nil)
  91 ;       ((affine-polynomialp expr) (list-variables1 expr))
  92 ;       ((rational-functionp expr) (list-variables2 (car expr))
  93 ;        (list-variables2 (cdr EXPR)))
  94 ;       (t(loop for v on expr while (listp v) do (list-variables2 (car v))))))
  95 ;
  96 ;(defun list-variables1 (u)
  97 ;   (declare (special a-list))
  98 ;  (cond ((atom u) niL)
  99 ;       (t
 100 ;        (push-new (car u) a-list)
 101 ;        (loop for (def cof) on (cdr u) by #'cddr
 102 ;              do (list-variables1 cof)))))
 103 ;
 104
 105 (defun pcoeff (poly monom &optional
 106                (variables-to-exclude-from-cof (list-variables monom)))
 107
 108   "x^3*y+x^2*u, x^2 ==> u . Rule: if a variable appears in monom it must be to
 109 the exact power, and if it is in variables to exclude it may not appear unless it
 110 was in monom to the exact power.  (pcoeff pol 1 ..) will exclude variables like
 111 substituting them to be zero."
 112   (cond ((eql monom 1) nil)
 113         ((atom monom) (setq monom (list monom 1 1))))
 114   (check-arg monom (or (listp monom)(eql monom 1)) "1 or a monomial")
 115   (cond ((pcoeff1 poly monom variables-to-exclude-from-cof ))
 116         (t 0)))
 117 (defun pcoeff1 (poly monom &optional
 118                 (variables-to-exclude-from-cof (list-variables monom))
 119                 &aux tem (true-deg 0))
 120   (cond
 121         ((and (numberp monom) (numberp poly)) poly)
 122         ((numberp poly) nil)
 123         (t (let ( (first-case t) )
 124              (cond ((numberp monom)
 125                     (cond ((member (p-var poly) variables-to-exclude-from-cof :test #'eq)
 126                                     (setq true-deg 0))
 127                           (t (setq first-case nil))))
 128                    ((eq (p-var monom) (p-var poly))
 129                     (setq true-deg (second monom))
 130                     (setq monom (caddr monom)))
 131                    ((member (p-var poly) variables-to-exclude-from-cof :test #'eq)
 132                     (setq true-deg 0))
 133                    (t (setq first-case nil)))
 134              (cond (first-case
 135                     (loop for (deg cof) on (cdr poly) by #'cddr
 136                           when (eq deg true-deg)
 137                           do
 138                           (return (pcoeff1 cof monom
 139                                               variables-to-exclude-from-cof))))
 140                    ((or (numberp monom)(pointergp (car poly) (car monom)))
 141                     (loop for (deg cof) on (cdr poly) by #'cddr
 142                           when (setq tem (pcoeff1 cof monom variables-to-exclude-from-cof))
 143                           collecting deg into temm
 144                           and collecting tem into temm
 145                           finally (cond (temm
 146                                          (cond ((eql (car temm) 0)
 147                                                 (return (second temm)))
 148                                                (t (return (cons (car poly) temm)))))))))))))
 149
 150
 151 (defun replacements (&optional (simps *poly-simplifications*))
 152   (loop for v in simps by #'cddr collecting v))
 153
 154 ;;These check for  occurrence of seq in poly.  But using this seems to
 155 ;;be slower than just checking the part-above just in case we will need
 156 ;;to replace. The only advantage is this does no consing.  But
 157 ;;with our new storage techniques..
 158 (defun occurs-in (poly seq)
 159  (cond ((eq (catch 'its-there (occurs-in1 poly seq)) 'ok)
 160         t)
 161        (t nil)))
 162
 163 (defun occurs-in1 (poly seq)
 164   (catch 'not-in-this-one
 165     (cond ((null seq)(throw 'its-there 'ok))
 166           ((numberp poly))
 167           (t (let ((var (first seq))
 168                    (deg (second seq)))
 169                (cond ((pointergp (car poly) var)
 170                       (loop for (deg cof) on (cdr poly) by #'cddr
 171                             do (occurs-in1 cof seq)))
 172                      ((eq (car poly) var)
 173                       (cond ((< (second poly ) deg)
 174                              (throw 'not-in-this-one nil))
 175                             (t (loop for (degg cof) on (cdr poly) by #'cddr
 176                                      while (>= degg deg)
 177                                      do (occurs-in1 cof (cddr seq))))))
 178                      (t (throw 'not-in-this-one nil))))))))
 179
 180 (defun convert-deg-sequence-to-monomial (sequ &aux tem)
 181   (cond ((null sequ) nil)
 182         ((numberp  sequ) sequ)
 183         (t (setq tem (convert-deg-sequence-to-monomial (cddr sequ)))
 184            (cond (tem (cond ((not (eql 0 (second sequ))) (list (car sequ) (second sequ) tem))
 185                             (t tem)))
 186                  (t (list (car sequ) (second sequ) 1))))))
 187
 188 (defun constant-deg-sequencep (seq)
 189   (loop for w in seq by #'cddr
 190         when (not (eql w 0))
 191         do (return nil)
 192         finally (return t)))
 193
 194 (defun $declare_constant_list( ll)
 195   (apply #'declare-constant (cdr ll)))
 196
 197 (defun declare-constant (&rest l)
 198   (loop for v in l
 199          when (get v 'disrep)
 200          do (putprop v t 'constant)
 201          else do
 202          (putprop  (car (num (new-rat v)))  t 'constant)))
 203
 204 (defun show-constants ()
 205   (loop for v in *genvar* when (get v 'constant)
 206         do (format t "~A disreps to ~A and is constant." v (get v 'disrep))))
 207
 208 (defun remove-constant-property (&rest a-list)
 209   (loop for u in a-list do
 210   (loop for v in *varlist*
 211         for w in *genvar*
 212         when (or (eq v u) (eq w u))
 213         do (remprop w 'constant))))
 214
 215 (defvar *use-constants* t)
 216
 217 (defun constant-functionp( poly)
 218   "Assumes that the variables have been ordered so that
 219   the constant variables are the least main"
 220    (or (numberp poly)(numberp (car poly))
 221        (and *use-constants* (symbolp (car poly)) (get (car poly) 'constant))))
 222
 223
 224 (defun find-main-monomial (poly &aux cof seq)
 225   (cond ((numberp poly) (values nil poly))
 226         (t
 227          (setq seq (loop while poly
 228                          collecting (car poly)
 229                          collecting (second poly)
 230                          do (setq poly (third poly))
 231                          when (constant-functionp poly) do (setq cof poly poly nil)))
 232          (values seq cof))))
 233
 234 (defun find-least-main-monomial (poly &aux   tem mono-so-far cof)
 235   (loop while poly
 236         do
 237         (setq tem (last2 poly))
 238         (cond ((not  (and (null mono-so-far)(zerop (car tem)) (numberp (second tem))))
 239                (cond ((not (zerop (car tem)))
 240
 241                (setq mono-so-far (cons (car tem) (cons (car poly)
 242                                                        mono-so-far) )))))
 243               (t (setq tem (lastn 4 poly))(setq mono-so-far (cons (car tem)
 244                                                                      (cons (car poly)
 245                                                                            mono-so-far)))))
 246         (cond ((numberp (second tem))(setq poly nil)(setq cof (second tem)))
 247               (t (setq poly (second tem)))))
 248   (values (nreverse mono-so-far) cof))
 249
 250 ;;(defvar  *monomial-order-function* 'seq-lessp)
 251 ;;(defvar  *monomial-order-function* 'xyz-and-degree)
 252 (defvar  *monomial-order-function* 'xyz-order-sequence)
 253
 254 (defun deg-sequence (variables degrees)
 255   (loop for v in variables
 256         for w in degrees
 257         unless (eql w 0)
 258         collecting v
 259         and
 260         collecting w))
 261
 262 (defun factor-out-monomial (poly monomial-degs copy)
 263   (cond (copy (setq poly (copy-tree poly))))
 264   (cond ((or (numberp poly)     (numberp monomial-degs))
 265          poly)
 266         (t
 267          (let ((deg (get-altq (car poly)  monomial-degs)))
 268            (loop
 269              for v on poly by #'cddr
 270              while (third v)
 271              when deg
 272              do
 273              (setf (second v)(sub (second v) deg))
 274              do (setf (third v)(factor-out-monomial
 275                                  (third v)
 276                                  monomial-degs
 277                                  nil))
 278              finally (return poly))))))
 279
 280 (defun replace-monomial (poly replacement monomial-degs
 281                          &optional part-above &aux to-replace denom-replacement)
 282   "replacement is a poly not a rat-poly"
 283   (check-arg poly affine-polynomialp "polynomial")
 284   (check-arg replacement rational-functionp "rat'l function")
 285   (cond (part-above (setq to-replace part-above))
 286         (t  (setq to-replace (part-above-degree poly monomial-degs))))
 287   (let (
 288         answer tem )
 289     (cond
 290       ((and (numberp poly) (zerop poly))(rzero))
 291       (to-replace
 292 ;          (format t "~%Replacing occurrence of ~A in.." (disrep-list monomial-degs))
 293 ;          (sh poly)
 294 ;          (format t "~% by the replacement")
 295 ;          (cond (denom-replacement (sh (cons replacement denom-replacement)))
 296 ;                (t          (sh replacement)))
 297        (setq answer (pdifference poly to-replace))
 298        (setq denom-replacement (denom replacement))
 299        (setq answer (ptimes denom-replacement answer))
 300        (setq tem (ptimes
 301                    (num replacement)
 302                    (factor-out-monomial to-replace monomial-degs t)))   ;copy?
 303        (setq answer (pplus answer tem))
 304        (setq answer (maybe-ratreduce answer denom-replacement))
 305 ;          (format t "The result is..")
 306 ;          (sh  answer)
 307        answer)
 308       (t poly))))
 309
 310 (defun replace-monomial-rat  (poly replacement monomial-degs
 311                          &optional part-above &aux to-replace denom-replacement)
 312   "replacement is a poly not a rat-poly"
 313   (check-arg poly rational-functionp "rat'l function")
 314   (check-arg replacement rational-functionp "rat'l function")
 315   (cond (part-above (setq to-replace part-above))
 316         (t  (setq to-replace (part-above-degree (num poly) monomial-degs))))
 317   (let (
 318         answer tem )
 319     (cond
 320       ((and (numberp (num poly)) (zerop (num poly)))(rzero))
 321       (to-replace
 322 ;          (format t "~%Replacing occurrence of ~A in.." (disrep-list monomial-degs))
 323 ;          (sh poly)
 324 ;          (format t "~% by the replacement")
 325 ;          (cond (denom-replacement (sh (cons replacement denom-replacement)))
 326 ;                (t          (sh replacement)))
 327        (setq answer (pdifference (num poly) to-replace))
 328        (setq denom-replacement (denom replacement))
 329        (setq answer (ptimes denom-replacement answer))
 330        (setq tem (ptimes
 331                    (num replacement)
 332                    (factor-out-monomial to-replace monomial-degs t)))   ;copy?
 333        (setq answer (pplus answer tem))
 334        (cond ((eql 1 (denom poly))
 335               (setq answer (maybe-ratreduce answer  denom-replacement)))
 336              (t  (setq answer (maybe-ratreduce answer
 337                                                (ptimes (denom poly )
 338                                                        denom-replacement)))))
 339        answer)
 340       (t poly))))
 341
 342 (defvar *reduce* t)
 343 (defun maybe-ratreduce (a b)(cond (*reduce* (ratreduce a b))
 344                                   (t (cons a b))))
 345 (defun $add_to_poly_simplifications (rat-monom)
 346   (setq rat-monom ($new_rat rat-monom))
 347   (setq rat-monom (cdr rat-monom))
 348   (add-to-poly-simplifications rat-monom))
 349 (defvar *poly-simplifications* nil)
 350
 351 (defvar *degenerate-conditions* nil)
 352
 353
 354
 355 (defun add-to-poly-simplifications (rat-monom &optional (simplify t) &aux orig added
 356                                      replacement replacement-num)
 357   (setq orig rat-monom)
 358   (cond (*simplify-rhs* (setq rat-monom (polysimp rat-monom)))
 359         (t (setq rat-monom (new-rat rat-monom))))
 360   (cond ((numberp rat-monom)
 361          (cond ((or (not (zerop rat-monom))
 362                     (not  ($zerop orig))) (setq *poly-simplifications*
 363                                                 (list 1 (rzero))))
 364                (t *poly-simplifications*)))
 365         ((eql 0 (num rat-monom))  (princ ".")(setq simplify nil) *poly-simplifications*)
 366         ((or (numberp rat-monom) (and (listp rat-monom)
 367                                       (constant-functionp (num rat-monom))))
 368          (setq *poly-simplifications* (list 1 (rzero)))
 369          (let ((*print-level* 2))
 370          (format t  "There is a constant ~A in the *poly-simplifications*" rat-monom)))
 371         (t
 372          (multiple-value-bind (seq cof)(find-main-monomial (num rat-monom))
 373            (setq replacement-num (pdifference
 374                                    (ptimes cof
 375                                            (convert-deg-sequence-to-monomial
 376                                              seq))
 377                                    (num rat-monom)
 378                                    ))
 379            (cond ((and (not (numberp (denom rat-monom)))
 380                        (not (constant-functionp (denom rat-monom))))
 381                   (format t "Warning: a denominator was discared which might take the value 0" nil)))
 382            (setq replacement (ratreduce replacement-num cof))
 383 ;          (format t "~%Adding to simplifications..")
 384 ;         (sh (convert-deg-sequence-to-monomial seq))(format t "  ==>  ") (sh replacement)
 385            (cond ((not (numberp cof))
 386                   (pushnew cof *degenerate-conditions* :test 'equal)))
 387            (setq *poly-simplifications* (append *poly-simplifications*
 388                                                 (list seq replacement))))
 389         (setq added t)
 390         ))
 391   (cond (simplify  (simplify-poly-simplifications))
 392         (t *poly-simplifications*))
 393 ;  (cond (added (conclusionp)))      ;11/30/85, chou
 394   *poly-simplifications*
 395    )
 396
 397 (defun grobner-basis (list-polys &key homogeneous)
 398   (setq *poly-simplifications* nil)
 399   (setq *degenerate-conditions* nil)
 400   (loop for v in list-polys
 401         do (add-to-poly-simplifications v))
 402   (check-overlaps 40 :reset t :homogeneous homogeneous)
 403   *poly-simplifications*)
 404
 405
 406 ;(defun ideal-subsetp (id1 id2 &key (reset-basis t) &aux tem)
 407 ; (cond (reset-basis  (grobner-basis id2)))
 408 ;  (loop for v in id1
 409 ;       when (not (rzerop (setq tem  (polysimp v))))
 410 ;       do (format t "~A is not zero modulo the second ideal. It is congruent to ~A"
 411 ;                  v tem)
 412 ;       (return nil)
 413 ;       finally (return t)))
 414
 415 (defun ideal-subsetp (id1 id2 &key ( verbose t) (reset-basis t) &aux tem)
 416   (cond (reset-basis  (grobner-basis id2)))
 417   (loop for v in id1
 418         when (not (rzerop (setq tem  (polysimp v))))
 419         do
 420         (cond (verbose (format t "~A is not zero modulo the second ideal. It is congruent to ~A"
 421                                id2 tem)))
 422         (return nil)
 423         finally (return t)))
 424
 425 (defun $set_up_poly_simplifications(a-list)
 426   (check-arg a-list '$listp nil)
 427   (setq *poly-simplifications* nil)
 428   (loop for v in (cdr a-list)
 429         do
 430         (add-to-poly-simplifications (st-rat v))
 431   finally (return (convert-poly-simplifications-to-macsyma-form *poly-simplifications*))))
 432
 433 (defun $show_simps ()
 434   (convert-poly-simplifications-to-macsyma-form *poly-simplifications*))
 435
 436 ;;tried using occurs-in instead of part-above-degree but wasted time...
 437
 438 (eval-when
 439     (:load-toplevel :compile-toplevel)
 440   (fmakunbound 'must-replacep)
 441   (defun must-replacep ( poly &aux part-above monom  replace)
 442                                         ;  (declare (values part-above monom replace))
 443     (loop for (seq repl) on *poly-simplifications* by #'cddr
 444           when (setq part-above (part-above-degree poly seq))
 445           do (setq monom (convert-deg-sequence-to-monomial seq) replace repl)
 446           (return 'found))
 447     (setq *must-replace-data* (list part-above monom replace))
 448     (values part-above monom replace)))
 449
 450 (defvar *must-replace-data* nil)
 451 (defun monomial-finish-replace ( )
 452   (destructuring-let (((part-above mon repl) *must-replace-data*))
 453     (cons (ptimes (num  repl) (num (ratreduce part-above mon))) (denom repl))))
 454
 455 (defun sort-poly-simplifications ()
 456  (setq *poly-simplifications*
 457        (sort-grouped-list
 458          *poly-simplifications*
 459          2 #'(lambda (x y)
 460                (cond ((equal x y) nil)
 461                      (t (funcall *monomial-order-function* y x)))))))
 462
 463 (defun check-simp (a simped &optional (tag 'simps-not-equal) )
 464   (cond ((not (equal (polysimp a) simped))(break tag))))
 465 (defvar *show-grob* nil)
 466 (defvar *conclusion* nil)
 467 (defun conclusionp ()
 468   (cond (*conclusion*
 469          (format t "~%Checking conclusion... ")
 470          (cond ((zero-modulo-simplifications *conclusion*)
 471                 (format t "~%The conclusion already follows.")
 472                 (throw 'its-proven t))
 473                (t (format t "no good yet"))))))
 474
 475 (defun new-zerop (x)
 476   (or (eql x 0)
 477       (and (listp x) (eql (car x) 0))))
 478
 479 (defvar *simplify-rhs* t)
 480
 481 (defun simplify-poly-simplifications (&aux seq repl hi old-simps repl-simp try-conclusion
 482                                       resimplify mono dif leng replacement-of-mono)
 483   ;;  (if test-all-theorems (check-proof-time))
 484   (cond
 485     (*show-grob*
 486      (setq leng (length *poly-simplifications*))
 487      (format t "~%Starting to simplify ~A poly-simplifications " (truncate leng 2))
 488      (cond ((> 10 leng )
 489             (show-simps))
 490            (t (format t "~%Replacements are :") (show-replacements)))))
 491   (cond
 492     ((loop
 493         for v in *poly-simplifications* by #'cddr
 494         when (numberp v)
 495         do (setq *poly-simplifications*
 496                  (list 1 (rzero))) (return 'unit-ideal)))
 497     (t
 498      (sort-poly-simplifications)
 499      (loop
 500         ;;for (seq repl) on *poly-simplifications* by #'cddr
 501         for i from 2 by 2 while (<= i (length *poly-simplifications*))
 502         do
 503           (setq seq (nth (- i 2) *poly-simplifications*)
 504                 repl(nth (- i 1) *poly-simplifications*))
 505           (setq old-simps *poly-simplifications*)
 506           (setq *poly-simplifications* (append     ;;;nconac ok..
 507                                         (subseq *poly-simplifications* 0 (- i 2))
 508                                         (nthcdr i *poly-simplifications*)))
 509           (cond ((must-replacep (setq mono (convert-deg-sequence-to-monomial seq)))
 510                  (cond (*show-grob*  (sh mono) (format t "   is not  reduced")))
 511
 512                  (cond ((equal (setq replacement-of-mono
 513                                      ;;  (polysimp mono)
 514                                      (monomial-finish-replace))
 515                                repl)
 516                         (setq i (- i 2))
 517                         (format t "eliminated one")
 518                         ;;  (check-simp mono replacement-of-mono 'a)
 519                         )
 520                        (t (setq dif
 521                                 (pdifference      ;;;just need the numerator...
 522                                  (ptimes (num replacement-of-mono)
 523                                          (denom repl))
 524                                  (ptimes (num repl) (denom replacement-of-mono))))
 525
 526                           (cond (*simplify-rhs*
 527                                  (setq dif (polysimp dif))))
 528                           (cond ((not (new-zerop dif))
 529                                  (add-to-poly-simplifications   dif nil)
 530                                  (setq try-conclusion t)))
 531                           (setq       resimplify t)
 532                                         ; (format t "~%having to resimplify.." ) (sh dif)
 533                           (return 'added))))
 534                 ((and *simplify-rhs*
 535                       (must-replacep (num repl)))
 536                  ;;replace the replacement only so monom ok..
 537                  (setq repl-simp (polysimp (num repl)))
 538                  (setq repl-simp
 539                        (ratreduce (num repl-simp)
 540                                   (ptimes (denom repl-simp) (denom repl))))
 541                                         ;(check-simp (nth (1- i) old-simps) repl-simp 'c)
 542                  (setf (nth (1- i) old-simps ) repl-simp)
 543                  (setq *poly-simplifications* old-simps))
 544                 (t (setq *poly-simplifications* old-simps))))))
 545   ;;?? this was not in earlier version:
 546   (cond (resimplify (simplify-poly-simplifications)))
 547   ;;  (cond (try-conclusion (conclusionp)))  ;;moved to add-to-simps..
 548   *poly-simplifications*)
 549
 550
 551 (defun seq-op (seqa seqb &optional (op #'-))
 552   (cond ((null seqa) seqb)
 553         ((null seqb) seqa)
 554         ((eq (car seqa) (car seqb))
 555          (nconc (list (car seqa ) (funcall op (second seqa) (second seqb)))
 556                       (seq-op (cddr seqa) (cddr seqb) op)))
 557         ((pointergp (car seqa) (car seqb))
 558          (nconc (list (car seqa )  (second seqa))
 559                 (seq-op (cddr seqa)  seqb) op))
 560         (t (nconc (list (car seqb )  (second seqb))
 561                   (seq-op  seqa (cddr seqb) op)))))
 562 (defun seq-lessp (seqa seqb)
 563  (cond ((catch 'less (seq-op seqa seqb
 564                              #'(lambda (x y)
 565                                  (cond ((> x y)(throw 'less t))))))
 566         nil)
 567        (t t)))
 568
 569 (defun seq-minus (seqa)
 570   (setq seqa (copy-list seqa))
 571   (loop for v on (cdr seqa) by #'cddr
 572         do (setf (first v) (- (first v))))
 573   seqa)
 574 ;;see other one above
 575 ;(defun simplify-poly-simplifications ( &aux hi new monom new-simps)
 576 ;  "This works but is not too quick"
 577 ;  (setq *poly-simplifications* (sort-grouped-list
 578 ;                                *poly-simplifications*
 579 ;                                2 '(lambda (x y)
 580 ;                                     (cond ((equal x y) nil)
 581 ;                                           (t (funcall *monomial-order-function* y x))))))
 582 ;  (cond (*show-grob*  (format t "~%Starting to simplify ..") (show-simps)))
 583 ;  (loop for (deg repl) on *poly-simplifications* by #'cddr
 584 ;       for i from 2 by 2
 585 ;        do
 586 ;       (setq new-simps nil)
 587 ;       (let ((*poly-simplifications* (append     ;;;nconc ok..
 588 ;                                       (firstn (- i 2) *poly-simplifications*)
 589 ;                                       (nthcdr i *poly-simplifications*))))
 590 ;         (cond ((or (must-replacep (setq monom (convert-deg-sequence-to-monomial deg)))
 591 ;                    (must-replacep (num  repl)))
 592 ;                ;;should use values of must-replacep to save time
 593 ;                (setq new (pdifference  (num repl) (ptimes monom (denom repl))))
 594 ;                (setq new (polysimp new))
 595 ;                (setq new-simps (add-to-poly-simplifications new nil))
 596 ;                (show-simps new-simps))))
 597 ;       (cond (new-simps (setq *poly-simplifications* new-simps)
 598 ;                        (return 'found-one))))
 599 ;  (cond (new-simps (simplify-poly-simplifications)))
 600 ;  *poly-simplifications*)
 601
 602 (defun basis_from_simps(&optional(simps *poly-simplifications*))
 603   (cons '(mlist)
 604         (loop for (u v) on simps by #'cddr
 605                collecting (n- u v))))
 606
 607 (defun convert-poly-simplifications-to-macsyma-form (simps)
 608   (loop for (u v) on simps by #'cddr
 609         collecting (list '(mlist) (header-poly (convert-deg-sequence-to-monomial u))
 610                          (header-poly v)) into tem
 611         finally (return (cons '(mlist) tem))))
 612
 613
 614 (defun overlap-degree (seqa seqb &aux answer temm  tem)
 615  (setq answer (+ (loop for v on seqa by #'cddr
 616         when (and (setq tem (get-altq (car v) seqb))
 617                   (not (zerop (setq temm (min tem (second v))))))
 618         summing (- (second v) temm)
 619         else summing (second v))
 620       (degree-of-deg-sequence seqb))))
 621
 622 (defun overlap (seqa seqb &aux tem)
 623   (loop for v on seqa by #'cddr
 624         when (and (setq tem (get-altq (car v) seqb))
 625                   (not (zerop (setq tem (min tem (second v))))))
 626         collecting (car v)
 627         and
 628         collecting tem))
 629 (defun degree-of-deg-sequence (seq)
 630   (loop for v in (cdr seq)  by #'cddr
 631         summing v))
 632
 633 (defun replacement-sequences ()
 634   (loop for v in *poly-simplifications* by #'cddr collecting v))
 635 (defun show-replacements ()
 636   (displa ($poly_replacements)))
 637 (defun $polY_replacements()
 638   (loop for v in (replacements)
 639         collecting (header-poly (convert-deg-sequence-to-monomial v))
 640         into tem finally (return (cons '(mlist) tem))))
 641
 642 (defun degree-of-worst-monomial (poly)
 643   (degree-of-deg-sequence (find-main-monomial poly)))
 644 (defun show-simps (&optional (simps *poly-simplifications*))
 645   (displa (convert-poly-simplifications-to-macsyma-form simps)))
 646 (defvar  *overlaps-already-checked* nil)
 647 (defvar *homogeneous-relations* nil)
 648 (defvar *good-thru-degree* 0)
 649
 650 (defvar *timed-grobner-basis* nil)
 651 (defun timedbasis (initial-time)
 652   (and  *timed-grobner-basis*
 653         (<  *timed-grobner-basis*
 654             (- (get-internal-run-time) initial-time)
 655             )
 656         (progn (if-verbose (format t "~%Grobner basis took too long, didn't finish")
 657                (show-simps))
 658
 659                (throw 'took-too-long 'took-too-long))))
 660
 661 (defvar *show-complexity* t)
 662
 663 (defun check-overlaps (deg  &key reset  homogeneous (from-degree 1)
 664                        (add-to-simps t) (maybe-reset t)
 665                        (initial-time (get-internal-run-time))
 666                        &aux tem associator v)
 667   (cond (*homogeneous-relations*
 668          (show *homogeneous-relations*)
 669          (format t "Assuming homogeneous relations")
 670          (setq homogeneous t)))
 671
 672   (cond (reset (setq *overlaps-already-checked* nil))
 673         (maybe-reset
 674          (cond ((y-or-n-p "Reset  *overlaps-already-checked* to nil")
 675                 (setq  *overlaps-already-checked* nil)))))
 676   (loop for i from from-degree to deg do
 677        (block big-sue
 678          (loop for subl on (nreverse (replacement-sequences))
 679             do (setq v (car subl))
 680               (loop for (w) on (cdr subl)
 681                  when
 682                    (and w
 683                         (= (overlap-degree v w) i)
 684                         (setq tem (overlap v w));;could check without consing
 685                         (not (member (list v w) *overlaps-already-checked* :test #'equal)))
 686                  do
 687                    (timedbasis initial-time)
 688                    (setq associator (check-associative v tem w))
 689                    (push (list v w) *overlaps-already-checked* )
 690                    (cond ((rzerop associator)
 691                           (setq initial-time (+ initial-time 30))
 692                           (if-verbose (show (- (get-internal-run-time) initial-time))
 693                                       (format t " zero associator, ")))
 694                          (add-to-simps
 695                           (add-to-poly-simplifications associator)
 696                           (cond (*show-complexity*
 697                                  (terpri)
 698                                  (loop for v in (cdr *poly-simplifications*) by #'cddr
 699                                     do (format t " . ~a" (+ (pcomplexity (num v))
 700                                                             (pcomplexity (denom v)))))))
 701                           (conclusionp)
 702                           (cond ((null homogeneous)
 703                                  (setq from-degree (min from-degree (degree-of-worst-monomial (car associator))))))
 704                           (check-overlaps deg :homogeneous homogeneous
 705                                           :from-degree from-degree :initial-time initial-time
 706                                           :add-to-simps add-to-simps :maybe-reset nil)
 707                           (return-from big-sue 'done))
 708                          (t (break 'found-non-associative))))
 709             finally (return *poly-simplifications*)))
 710        (cond (homogeneous (setq from-degree i))))
 711   *poly-simplifications*)
 712
 713
 714 (defun $polysimp (f &aux answer)
 715   (setq answer  (case (poly-type f)
 716                          (:number f)
 717                          (:polynomial (polysimp f))
 718                          (:rational-function  (ratquotient (polysimp  (num f))
 719                                                            (polysimp (denom f))))
 720                          (t nil)))
 721   (cond (answer
 722          (header-poly answer))
 723         (t  (cond ((mbagp f)(cons (car f) (mapcar '$polysimp (cdr f))))
 724                   (t ($polysimp (new-rat ($totaldisrep f))))))))
 725
 726 (defun check-associative (a overlap b &aux answer simpa simpb)
 727   (when *show-grob*
 728     (format t "Checking overlap ~A for ~A ~A "(disrep-list overlap) (disrep-list a) (disrep-list b)))
 729   (let ((ma (convert-deg-sequence-to-monomial a))
 730         ( mo (convert-deg-sequence-to-monomial overlap))
 731         (mb (convert-deg-sequence-to-monomial b)))
 732                                         ;   (setq simpa (polysimp (cons ma 1)))
 733     (setq simpa  (get-altq a *poly-simplifications*))
 734     (setq simpb  (get-altq b *poly-simplifications*))
 735                                         ;    (setq simpb (polysimp (cons mb 1)))
 736
 737     (setq simpa (cons (ptimes (num simpa) (num (ratreduce mb mo))) (denom simpa)))
 738     (setq simpb (cons (ptimes (num simpb) (num (ratreduce ma mo))) (denom simpb)))
 739
 740     ;;everything you're just shifting a numerator there
 741     (cond ((equal simpa simpb)
 742            (princ ".")
 743            (cond (*show-grob*   (format t"..The overlap is associative")(rzero)))
 744            (rzero))
 745           (t
 746            (setq answer (ratdifference simpa simpb))
 747            (cond ((null *simplify-rhs* ) (setq plain-answer answer)
 748                   (setq answer (nred (polysimp (function-numerator answer))
 749                                      (function-denominator answer)))))
 750            (cond ((null *simplify-rhs*)
 751                   (cond ((new-zerop answer))
 752                         (t (setq answer plain-answer)))))
 753            (cond (*show-grob*
 754                   (format t "~%The difference  is ..")(sh answer)))
 755            (cond ((new-zerop answer) (princ ".")(rzero))
 756                  (t answer))))))
 757
 758 (defvar *simplify-simplifications* nil)
 759
 760 (defun reverse-z-order-sequence (seqa seqb)
 761   (and (not (equal seqa seqb))(not (z-order-sequence seqb seqa))))
 762
 763 (defun z-order-sequence (seqa seqb)
 764    (loop for v on seqa by #'cddr
 765          for w on seqb by #'cddr
 766          when (or (< (symbol-value (car v)) (symbol-value (car w)))
 767                   (< (second v)(second w)))
 768          do (return nil)
 769          finally (return (not (equal v w)))))
 770
 771 (defun n** (&rest l)
 772   (case (length l)
 773     (0 1)
 774     (1 (car l))
 775     (2 (n* (car l) (cadr l)))
 776     (t (n* (car l) (apply 'n** (cdr l))))))
 777
 778 ;(setq me #$[x,x^2,z,z^2,z^2*x,z^2*x*t]$)
 779 (defun xyz-order-sequence (seqa seqb)
 780   " so x<x^2<z<z^2<z^2*x<z^2*x*t the worst monomial  is the most main monomial ie
 781   highest degree in z then highest in y etc. "
 782   (loop for (var deg) on seqb by #'cddr while seqa
 783      do
 784        (cond ((eq (car seqa) var)
 785               (cond ((> (second seqa) deg) (return nil))
 786                     ((< (second seqa) deg) (return t))
 787                     (t (setq seqa (cddr seqa)))))
 788              ((pointergp (car seqa) var) (return nil))
 789              (t (return t)))
 790      finally (cond ((null seqa) (return (and seqb)))
 791                    (t (return nil)))))
 792
 793
 794 (defun degree-of-sequence (seq)
 795   (loop for v on (cdr seq) by #'cddr
 796         summing (car v)))
 797 (defun xyz-and-degree (seqa seqb &aux (d1 (degree-of-sequence seqa))
 798                           (d2 (degree-of-sequence seqb)))
 799   (cond ((eq d1 d2)(xyz-order-sequence seqa seqb))
 800         ((< d1 d2) t)
 801         (t nil)))
 802
 803 (defun remove-0 (lis)
 804   (loop for (var deg) on lis by #'cddr
 805         unless (zerop deg)
 806         collecting var and collecting deg))
 807
 808 ;(defun xyz-order-sequence (seqa seqb)
 809 ;  " so x<x^2<z<z^2<z^2*x<z^2*x*t the worst monomial  is the most main monomial ie
 810 ;  highest degree in z then highest in y etc. "
 811 ;  (loop for v on seqa by #'cddr
 812 ;       for w on seqb by #'cddr
 813 ;       do
 814 ;       (cond ((eq (car v) (car w))
 815 ;              (cond ((eq (second w)(second v))nil)
 816 ;                    ((< (second v)(second w))(return t))
 817 ;                    (t (return nil))))
 818 ;             ((< (symbol-value (car v)) (symbol-value  (car w)))(return t))
 819 ;             (t (return nil)))))
 820
 821 (defun reverse-xyz-order-sequence (seqa seqb)
 822   (and (not (equal seqa seqb))(not (xyz-order-sequence seqb seqa))))
 823
 824 ;(defun polysimp (poly &aux (changes t) changed tem repl (replaced-th 0))
 825 ;  (cond ((affine-polynomialp poly)(setq poly (cons poly 1)))
 826 ;       ((rational-functionp poly) nil)
 827 ;       (t (setq poly (new-rat poly))))
 828 ;  (cond ((constant-functionp poly) (values  poly nil))
 829 ;       (t
 830 ;        (loop while changes
 831 ;              do (setq changes nil)
 832 ;              (loop for v on *poly-simplifications* by #'cddr
 833 ;                    for i from 0 by 2
 834 ;                    when (setq tem (part-above-degree (num poly) (car v)))
 835 ;                    do (setq changes t changed t)
 836 ;                    (cond ((< i replaced-th)
 837 ;                           (format t "~%Replacing by earlier ~A before ~A"
 838 ;                                 (disrep-list  (nth i *poly-simplifications* ))
 839 ;                                (disrep-list   (nth replaced-th *poly-simplifications* ))
 840 ;                                (break t))))
 841 ;
 842 ;                     (setq replaced-th i)
 843 ;
 844 ;                    (setq repl (replace-monomial-rat  poly (second v) (car v) tem ))
 845 ;                    (cond ((eql (num repl) 0)(setq changes nil)))
 846 ;                    (setq poly repl)
 847 ;                    (return 'changed))
 848 ;              finally (return 'done))
 849 ;        (values poly changed))))
 850
 851 (defun zero-modulo-simplifications (poly )
 852   (pzerop (polysimp-ignore-denominator poly)))
 853
 854
 855 (defun polysimp-ignore-denominator  (poly &aux (changes t) *reduce* changed tem repl)
 856   "Returns something G where G*unit= poly modulo simplifications"
 857   (setq poly
 858         (cond ((affine-polynomialp poly)(cons poly 1))
 859               ((rational-functionp poly) (cons (num poly ) 1))
 860               (t  (cons (num (new-rat poly)) 1))))
 861   (cond ((constant-functionp poly) nil)
 862         (t
 863          (loop while changes
 864                do (setq changes nil)
 865                (loop for v on *poly-simplifications* by #'cddr
 866                      when (setq tem (part-above-degree (num poly) (car v)))
 867                      do (setq changes t changed t)
 868                      (setq repl (replace-monomial-rat  poly (second v) (car v) tem ))
 869                      (cond ((eql (num repl) 0)(setq changes nil)))
 870                      (setq poly (cons (num repl) 1))
 871                      (return 'changed))
 872                finally (return 'done))))
 873   (function-numerator poly))
 874
 875 (defun polysimp (poly &aux (changes t) changed tem repl)
 876   (cond ((affine-polynomialp poly)(setq poly (cons poly 1)))
 877         ((rational-functionp poly) nil)
 878         (t (fsignal "bad-poly ~a" poly) (setq poly (new-rat poly))))
 879   (cond ((constant-functionp poly)
 880          (cond ((get-altq 1 *poly-simplifications*)(values 0 nil))
 881                (t (values  poly nil))))
 882         (t
 883          (loop while changes
 884                do (setq changes nil)
 885                (loop for v on *poly-simplifications* by #'cddr
 886                      when (setq tem (part-above-degree (num poly) (car v)))
 887                      do (setq changes t changed t)
 888                      (setq repl (replace-monomial-rat  poly (second v) (car v) tem ))
 889                      (cond ((eql (num repl) 0)(setq changes nil)))
 890                      (setq poly repl)
 891                      (return 'changed))
 892                finally (return 'done))
 893          (values poly changed))))
 894 ;(defun polysimp (poly &aux (changes t) changed tem repl)
 895 ;  (cond ((affine-polynomialp poly)(setq poly (cons poly 1)))
 896 ;       ((rational-functionp poly) nil)
 897 ;       (t (setq poly (new-rat poly))))
 898 ;  (cond ((constant-functionp poly) (values  poly nil))
 899 ;       (t
 900 ;        (loop while changes
 901 ;              do (setq changes nil)
 902 ;              (loop for v on *poly-simplifications* by #'cddr
 903 ;                    when (setq tem (part-above-degree (num poly) (car v)))
 904 ;                    do (setq changes t changed t)
 905 ;                    (setq repl (replace-monomial (num poly) (second v) (car v) tem ))
 906 ;                    (cond ((eql (num repl) 0)(setq changes nil)))
 907 ;                    (setq poly (ratreduce (num repl) (ptimes (denom poly) (denom repl))))
 908 ;                    (return 'changed))
 909 ;              finally (return 'done))
 910 ;        (values poly changed))))
 911 ;
 912 ;(defun polysimp (poly &aux (changes t) tem)
 913 ;  (cond ((affine-polynomialp poly)nil)
 914 ;       ((rational-functionp poly) nil)
 915 ;       (t (setq poly (new-rat poly))))
 916 ;  (loop while changes
 917 ;       do (setq changes nil)
 918 ;       (loop for v on *poly-simplifications* by #'cddr
 919 ;             when (setq tem (part-above-degree poly (car v)))
 920 ;             do (setq changes t)
 921 ;             (setq poly (replace-monomial poly (second v) (car v)))
 922 ;             (return 'changed))
 923 ;       finally (return poly)))
 924
 925
 926 (DEFUN CONVERT-MONOMIAL-TO-DEG-SEQUENCE (MONOM &AUX TEM)
 927     (COND ((OR (ATOM MONOM)
 928                (NULL MONOM))
 929            NIL)
 930           (T
 931            (SETQ TEM (CONVERT-MONOMIAL-TO-DEG-SEQUENCE (THIRD MONOM)))
 932            (APPEND (LIST (CAR MONOM) (SECOND MONOM)) TEM))))
 933
 934 ;
 935 ;;;could introduce *remaining-poly-simplifications* which is set by rat-polysimp
 936 ;;;initially  to all *poly-simplifications* and then gradually reduced by rat-polysimp-once
 937 ;(defun rat-polysimp-once (poly &aux new-denom tem)
 938 ;  (cond ((not (numberp poly))
 939 ;        (loop for (deg-seq repl) on *poly-simplifications* by #'cddr
 940 ;              when (setq tem (part-above-degree poly deg-seq))
 941 ;              do
 942 ;              (setq new-denom (denom repl))
 943 ;              (setq poly (replace-monomial poly  repl deg-seq))
 944 ;                    (return 'changed))
 945 ;        (values poly new-denom))
 946 ;       (t poly)))
 947 ; ;;new-denom will be nil unless there was a changed
 948 ;
 949 ;
 950 ;;;;we just use polysimp now.  it accepts a rat or non rat f'n checking...
 951 ;
 952 ;(defun rat-polysimp (rat-poly &aux (changes t)  a b changed-since-reduction
 953 ;                    changed-at-least-once answer)
 954 ;;  (format t "~%beginning to rat-polysimp..") (sh rat-poly)
 955 ;  (setq a (num rat-poly) b (denom rat-poly))
 956 ;  (setq answer
 957 ;       (loop while changes
 958 ;             do (setq changes nil)
 959 ;             (multiple-value-bind (new-poly changed-denom)
 960 ;                 (rat-polysimp-once  a)
 961 ;               (cond (changed-denom (setq changes t)(setq changed-since-reduction t)
 962 ;                                    (setq changed-at-least-once t)
 963 ;                                    (setq b (ptimes changed-denom b))
 964 ;                                    (setq a new-poly)
 965 ;                                    (cond ((numberp b) nil)
 966 ;                                          (t (setq rat-poly (ratreduce a b))
 967 ;                                             (setq changed-since-reduction nil)
 968 ;                                             (setq a (num rat-poly)  b (denom rat-poly)))))))
 969 ;             finally(cond (changed-since-reduction (return (ratreduce a b)))
 970 ;                          (t (return rat-poly)))))
 971 ;;  (format t "~%ending rat-polysimp with answer...")  (sh answer)
 972 ;  (values answer changed-at-least-once))
 973
 974 (defun $current_grobner_basis (&optional (simps *poly-simplifications*))
 975   (loop for (seq repl) on simps by #'cddr
 976         when (numberp seq) do (return '((mlist) 1))
 977         collecting (header-poly
 978                      (cons  (pdifference
 979                               (num repl)
 980                               (ptimes (convert-deg-sequence-to-monomial seq)
 981                                       (denom repl)))
 982                             1))
 983         into tem
 984         finally (return (cons '(mlist) tem))))
 985
 986 (defun $grobner_basis (ideal &optional no-concl homogeneous (upto 100))
 987   (cond (no-concl (setq *conclusion* nil)))
 988   (check-arg ideal $listp "macsyma list")
 989
 990   ($set_up_poly_simplifications ideal)
 991 ;  (show-simps)
 992   (conclusionp)
 993   (check-overlaps upto :reset t :homogeneous homogeneous)
 994   ($current_grobner_basis))
 995
 996 (defun delete-pair (first-item a-list )
 997   "deletes first occurrence"
 998   (loop
 999      for w on a-list by #'cddr
1000      when (not (equal (car w) first-item))
1001      collecting (car w) into tem
1002      and
1003      collecting (second w) into tem
1004      else do (return  (nconc tem (cddr w)))))
1005
1006 (defun seq-subset (seqa seqb &aux answer)
1007   (setq answer (catch 'not-in
1008                  (seq-op seqa seqb #'(lambda (.a. b) (cond ( (<= .a. b)
1009                                                             (- b .a.))
1010                                                            (t (throw 'not-in 'no)))))))
1011   (cond ((eq answer 'no) nil)
1012         (t answer)))
1013
1014
1015
1016 (defvar *old-genvar* nil)
1017
1018 ;;this is dangerous for some reason...
1019 (defun substitute-real-symbols-for-genvar (&optional a-list)
1020   (cond ((null a-list)
1021         (setq a-list (loop for v in *genvar* collecting
1022                            (intern (string-trim "$" (get v 'disrep)))))))
1023   (check-arg a-list  (eql (length a-list) (length *genvar*)) "not right length")
1024   (loop for v in *genvar*
1025         for w in a-list
1026         do
1027         (setf (symbol-value w) (symbol-value v))
1028         (setf (symbol-plist w) (symbol-plist v))
1029         finally (setq *genvar* a-list *old-genvar* *genvar*)))
1030
1031 (defun $generate_matrix (funct rows cols)
1032   (loop for i from 1 to rows
1033         collecting
1034         (loop for j from 1 to cols
1035               collecting (mfuncall funct i j)
1036               into tem
1037               finally (return (cons '(mlist) tem)))
1038         into temm
1039         finally (return (cons '($matrix) temm))))
1040
1041 (defun $nilpotent ($f &aux (old-simps *poly-simplifications*) )
1042  (unwind-protect
1043    (let (*conclusion* (me (sub* (mul* $f '$ggggg) 1)) mee)
1044     (setq mee (new-rat me))
1045     (add-to-poly-simplifications mee)
1046     (check-overlaps 10 :add-to-simps t :reset t)
1047       (cond ((eql 0 ($polysimp '$ggggg)) t)
1048             (t nil)))
1049       (setq *poly-simplifications* old-simps)))
1050
1051 (defun rat-variable( a-var)
1052   (loop for v in *varlist*
1053         for w in *genvar*
1054         when (eq v a-var)
1055         do (return (list w 1 1))
1056         finally (show a-var  v w)(new-rat (mul* a-var 2))
1057         (return (rat-variable a-var))))
1058
1059 (defun must-replace-monom (monom repl)
1060   (loop for seq in repl
1061         when (occurs-in monom seq)
1062         do (return t)))
1063
1064 (defun sort-list-of-monomials (a-list &aux grouped)
1065  (setq grouped
1066        (loop for v in a-list collecting (convert-monomial-to-deg-sequence v)
1067         collecting v))
1068        (setq grouped (sort-grouped-list grouped 2 *monomial-order-function*))
1069        (loop for w in (cdr grouped) by #'cddr collecting w))
1070 ;
1071 ;(user:defremember grobner-monomials(variables n &optional (replacements
1072 ;                                                            (replacements))
1073 ;                                              reset
1074 ;                         )
1075 ;  "Note that in this code the replacements are mereley for the remember feature,
1076 ; and *poly-simplifications* are the ones used for the must-replacep.  The version
1077 ; below is similar but just uses occurs-in instead of must-replacep"
1078 ;  (let (rat-vars tem vv)
1079 ;  (check-arg variables $listp "macsyma list")
1080 ;  (cond ((= n 0)(cond ((member 1 replacements) nil)
1081 ;                     (t (list 1))))
1082 ;       ((= n 1)
1083 ;        (setq rat-vars (loop for v in (cdr variables)
1084 ;                             do (setq vv (rat-variable v))
1085 ;                             when (not (must-replacep  vv))
1086 ;                             collecting vv ))
1087 ;        (setq rat-vars (sort rat-vars #'(lambda (x y) (pointergp (car x) (car  y)))))
1088 ;        rat-vars)
1089 ;       (t (loop for w in (grobner-monomials variables (1- n) replacements reset)
1090 ;                appending
1091 ;                (loop for uu in (grobner-monomials variables 1 replacements reset)
1092 ;                when (and
1093 ;                       (>= (symbol-value (car uu))
1094 ;                           (symbol-value (car w)))
1095 ;                       (not  (must-replacep (setq tem (ptimes w uu)) )))
1096 ;                collecting tem))))))
1097
1098
1099 (defun $poly_monomial_dimensions (variables to-deg &aux tem)
1100   (loop for i from 0 to to-deg
1101         do (format t "~%There are ~A independent monomials in degree ~A"
1102                    (setq tem (length (grobner-monomials variables i))) i)
1103                 summing tem into tot
1104         finally (format t "~% The total is ~A." tot)(return tot)))
1105
1106 (defremember grobner-bi-monomials(variables-y variables-x n-y n-x &optional (replacements
1107                                                              (replacements))
1108                                                reset
1109                           )
1110   (cond (reset (clear-memory-function 'grobner-monomials)))
1111   (let ( tem )
1112   (check-arg variables-y $listp "macsyma list")
1113   (cond ((= n-x 0)
1114          (grobner-monomials variables-y n-y replacements reset))
1115         ((= n-y 0)
1116          (grobner-monomials variables-x n-x replacements reset))
1117         (t (loop for w in (grobner-bi-monomials variables-y variables-x  (1- n-y) n-x replacements reset)
1118                  appending
1119                  (loop for uu in (grobner-monomials variables-y 1 replacements reset)
1120                  when (and
1121                         (>= (symbol-value (car uu))
1122                             (symbol-value (car w)))
1123                         (not  (must-replace-monom (setq tem (ptimes w uu)) replacements)))
1124                  collecting tem))))))
1125
1126 (defremember grobner-monomials(variables n &optional (replacements
1127                                                              (replacements))
1128                                                reset
1129                           )
1130   (cond (reset (clear-memory-function 'grobner-monomials)))
1131   (let (rat-vars tem vv)
1132   (check-arg variables $listp "macsyma list")
1133   (cond ((= n 0)(cond ((member 1 replacements) nil)
1134                       (t (list 1))))
1135         ((= n 1)
1136          (setq rat-vars (loop for v in (cdr variables)
1137                               do (setq vv (rat-variable v))
1138                               when (not (must-replace-monom  vv replacements))
1139                               collecting vv ))
1140          (setq rat-vars (sort rat-vars #'(lambda (x y) (pointergp (car x) (car  y)))))
1141          rat-vars)
1142         (t (loop for w in (grobner-monomials variables (1- n) replacements reset)
1143                  appending
1144                  (loop for uu in (grobner-monomials variables 1 replacements reset)
1145                  when (and
1146                         (>= (symbol-value (car uu))
1147                             (symbol-value (car w)))
1148                         (not  (must-replace-monom (setq tem (ptimes w uu)) replacements)))
1149                  collecting tem))))))
1150
1151 (defun $grobner_monomials (variables n &optional (replacements
1152                                                              (replacements))
1153                                                reset)
1154   (cons '(Mlist) (mapcar 'new-disrep (grobner-monomials variables n replacements reset))))
1155
1156
1157 (defun $list_relations (&rest l)
1158   (loop for v in l
1159         when ($listp v)
1160         appending (cdr (apply '$list_relations (cdr v))) into tem
1161         else
1162          when (equal (caar v) 'mequal)
1163          collecting (sub* (second v) (third v))into tem
1164           else collecting v into tem
1165           finally (return (cons '(mlist) tem))))
1166
1167 (defun rat-degree (poly var)
1168   (cond ((numberp poly)0)
1169         ((eq (car poly) var)(second poly))
1170         (t (loop for (deg cof) on (cdr poly) by #'cddr
1171                  maximize (rat-degree cof var)))))
1172
1173 (defun rational-sub  (poly var new-denom &aux numerator)
1174   "to do y-->y/u where u is earlier in the order and u is a poly"
1175   (let ((deg (rat-degree poly var)))
1176     (setq poly (copy-tree poly))
1177     (setq numerator (rational-sub1 poly var new-denom deg))
1178     (values poly new-denom deg)))
1179
1180 (defun rational-sub1 (poly var new-denom degree)
1181   (cond ((atom poly)(setq poly (ptimes poly (pexpt new-denom degree) )))
1182         ((eq (car poly) var)
1183          (loop for tail  on (cdr poly) by #'cddr
1184                 for (deg cof) on (cdr poly) by #'cddr
1185                 do (setf (second tail)
1186                          (ptimes cof (pexpt new-denom (- degree deg))))))
1187         (t (loop for tail on (cdr poly) by #'cddr
1188                   do (setf (second tail)
1189                            (rational-sub1 (second tail) var new-denom degree)))))
1190   Poly)
1191
1192 (defun $te (poly var repl)
1193   (header-poly (rational-sub (num (new-rat poly)) (caar (new-rat var)) (num (new-rat repl)))))
1194
1195 (defun $te (poly var)
1196   (rat-degree(num (new-rat poly )) (caar (new-rat var))))
1197
1198 (defun $Contract_to_variables (a-list upto-variable &optional (before-var nil) &aux tem)
1199    (check-arg a-list $listp "Macsyma list")
1200    (cond ((null before-var)(setq before-var (+ (symbol-value (get-genvar upto-variable)) .5)))
1201          (t (setq before-var  (symbol-value  (get-genvar before-var)))))
1202    (loop for v in (cdr a-list)
1203          do
1204          (setq tem nil)
1205          (case (poly-type v)
1206            (:polynomial (setq tem (car v)))
1207            (:rational-function (setq tem (caar v)))
1208            (:$rat (setq tem (caadr v))))
1209          (cond ((> before-var (symbol-value tem))
1210                 (setq tem v))
1211                (t (setq tem nil)))
1212          when tem
1213          collecting tem into llist
1214          finally (return (cons '(mlist ) llist))))
1215 ;
1216 ;(defun st-rat (x &aux answer)
1217 ;  "Returns a a number or cre polynomial or rational-function corresponding to x
1218 ;where x may be a polynomial"
1219 ;  (cond ((and (listp x) (mbagp x))
1220 ;        (cons (car x) (mapcar 'st-rat (cdr x)))))
1221 ;  (cond ((case (poly-type x)
1222 ;          (:number t)
1223 ;          (:polynomial (member (car x) *genvar* :test #'eq))
1224 ;          (:rational-function (member (caar x) *genvar* :test #'eq))
1225 ;          (:$rat (member (caadr x) *genvar* :test #'eq)(setq x (cdr x)))
1226 ;          (t nil))
1227 ;        (cond ((numberp x) x)
1228 ;              ((and (rational-functionp x)(eql (denom x) 1))(num x))
1229 ;              (t x)))
1230 ;       (t (cond ((member (poly-type x) '(:rational-function :polynomial :number) :test #'eq)
1231 ;                 (merror "not in standard rat form ~A " x))
1232 ;                (t (setq answer (new-rat x))(show (denom answer))
1233 ;                   (cond ((and (rational-functionp answer)
1234 ;                               (eql (denom answer) 1)) (num answer))
1235 ;                         (t answer)))))))
1236
1237 (defun st-rat (x)
1238   (cond ((and (listp x)(listp (car x)) (mbagp x))
1239           (mapcar 'st-rat1 (cdr x)))
1240         (t (st-rat1 x))))
1241
1242 (defun st-rat1 (x)
1243   (cond ((affine-polynomialp x) x)
1244         ((rational-functionp x) (cond ((eql (denom x) 1)(num x))
1245                                       (t x)))
1246         (($ratp x) (cond ((numberp (num (cdr x))) (st-rat (cdr x)))
1247                          ((member  (caadr x) *genvar* :test #'eq)(st-rat (cdr x)))
1248                          (t (st-rat ($totaldisrep x)))))
1249         (t (st-rat (new-rat x)))))
1250
1251 (defun get-genvar (var)
1252   (caar (new-rat var)))
1253
1254 ;;Thus  h is in ideal1 ^R ideal2 iff h is in ideal1*zz+ideal2*(1-zz)
1255 ;;the ==> direction is trivial and
1256 ;; for <==  h(x,y)=i1zz+i2(1-zz) and h is free of zz then i1=i2 and h=i2.
1257 (defun $ideal_intersection (ideal1 ideal2)
1258    "Returns the grobner basis of the intersection of the ideals generated by the
1259   two lists ideal1 and ideal2"
1260    (check-arg ideal1 '$listp "macsyma list")
1261   (let ((new-var (intern "ZZZZZ1"))
1262         gen     basis)
1263     (setq gen (get-genvar new-var))
1264     (setq basis (append
1265                   (loop for w in (cdr ideal1)
1266                         collecting (n* new-var w))
1267                   (loop for w in (cdr ideal2)
1268                         collecting (n* (n- 1 new-var ) w))))
1269     (setq basis ($grobner_basis (cons '(mlist)basis)))
1270     ($contract_to_variables basis nil new-var)))
1271
1272 (defun $grobner_basis_of_slice (slices original-relations &aux relats)
1273   (check-arg slices '$listp "macsyma list")
1274   (check-arg original-relations '$listp "macsyma list")
1275   (setq relats
1276         (append
1277           (loop for u in (cdr slices )
1278                 collecting (bring-to-left-side u))
1279           (cdr original-relations)))
1280   (displa (setq relats (cons '(mlist) relats)))
1281   ($grobner_basis  relats))
1282
1283 ;
1284 ;(DEFUN TE ($X)
1285 ;    (MEVAL `((MPLUS) $X 3)))
1286 ;
1287 ;
1288 ;(defmacro vb (&rest l)
1289 ;  (loop for v in l
1290 ;       collecting (list 'eval? v) into tem
1291 ;    finally (return (cons 'list  tem))))
1292 ;
1293 ;
1294 ;(defmacro eval? (ssymbol)
1295 ;  `(cond ( ,(list 'variable-boundp ssymbol) ,ssymbol)
1296 ;        (t (print 'special) ',ssymbol)))
1297 ;(defmacro teee (x)
1298 ;  (funcall 'variable-boundp x)
1299 ;  '(cons 1 1))
1300 ;
1301 ;
1302 ;(hii 3)
1303 ;(defun te (b)
1304 ;  (teee b)
1305 ;  (vb b x))
1306 ;
1307 ;si:
1308 ;(defun describe-region (region &key octal  (stream standard-output) &aux bits)
1309 ;  (cond (octal (let ((sstring  (format nil "~A" region))
1310 ;                    (ibase #10r8))
1311 ;                (setq region (read-from-string sstring)))))
1312 ;
1313 ;            (SETQ BITS (REGION-BITS REGION))
1314 ;            (FORMAT stream
1315 ;                    "~%  Region #~O: Origin ~:O, Length ~:O, Free ~:O, GC ~:O, Type ~A ~A, ~[~;Read only, ~]~[~;Temp, ~]~@[Swapin ~O, ~]~[NoScav~;Scav~]~[~; Safeguard~]~[~; NoCons~]~[~; Ephemeral~]~:[ Level=~D~]~%"
1316 ;                    REGION (REGION-ORIGIN REGION) (REGION-LENGTH REGION)
1317 ;                    (REGION-FREE-POINTER REGION) (REGION-GC-POINTER REGION)
1318 ;                    (NTH (LDB %%REGION-REPRESENTATION-TYPE BITS)
1319 ;                         '(LIST STRUC "REP=2" "REP=3"))
1320 ;                    (OR (NTH (LDB %%REGION-SPACE-TYPE BITS)
1321 ;                             '(FREE OLD NEW COPY STATIC STACK))
1322 ;                        (FORMAT NIL "TYPE=~O" (LDB %%REGION-SPACE-TYPE BITS)))
1323 ;                    (LDB %%REGION-READ-ONLY BITS)
1324 ;                    (LDB %%REGION-TEMPORARY BITS)
1325 ;                    (AND (PLUSP (LDB %%REGION-SWAPIN-QUANTUM BITS))
1326 ;                         (LDB %%REGION-SWAPIN-QUANTUM BITS))
1327 ;                     (LDB %%REGION-SCAVENGE-ENABLE BITS)
1328 ;                    (LDB %%REGION-SAFEGUARD BITS)
1329 ;                    (LDB %%REGION-NO-CONS BITS)
1330 ;                    (LDB %%REGION-EPHEMERAL BITS)
1331 ;                    (ZEROP (LDB %%REGION-LEVEL BITS))
1332 ;                    (LDB %%REGION-LEVEL BITS)))
1333 ;
1334 ;
1335 ;si:
1336 ;(DEFUN DESCRIBE-AREA (AREA &key
1337 ;                     (stream standard-output) &AUX  LENGTH USED N-REGIONS)
1338 ;  (AND (NUMBERP AREA) (SETQ AREA (AREA-NAME AREA)))
1339 ;  (DO AREA-NUMBER 0 (1+ AREA-NUMBER) (= AREA-NUMBER (length (the lisp::array #'AREA-NAME)))
1340 ;    (COND ((EQ AREA (AREA-NAME AREA-NUMBER))
1341 ;          (MULTIPLE-VALUE (LENGTH USED N-REGIONS) (ROOM-GET-AREA-LENGTH-USED AREA-NUMBER))
1342 ;          (FORMAT stream "~%Area #~O: ~S has ~D region~P, max size ~:O, region size ~:O (octal):~%"
1343 ;                    AREA-NUMBER AREA N-REGIONS N-REGIONS
1344 ;                    (AREA-MAXIMUM-SIZE AREA-NUMBER) (AREA-REGION-SIZE AREA-NUMBER))
1345 ;          (DO ((REGION (AREA-REGION-LIST AREA-NUMBER) (REGION-LIST-THREAD REGION))
1346 ;               (BITS))
1347 ;              ((MINUSP REGION))
1348 ;            (SETQ BITS (REGION-BITS REGION))
1349 ;            (FORMAT stream "  Region #~O: Origin ~:O, Length ~:O, Free ~:O, GC ~:O, Type ~A ~A, ~[~;Read only, ~]~[~;Temp, ~]~@[Swapin ~O, ~]~[NoScav~;Scav~]~[~; Safeguard~]~[~; NoCons~]~[~; Ephemeral~]~:[ Level=~D~]~%"
1350 ;                    REGION (REGION-ORIGIN REGION) (REGION-LENGTH REGION)
1351 ;                    (REGION-FREE-POINTER REGION) (REGION-GC-POINTER REGION)
1352 ;                    (NTH (LDB %%REGION-REPRESENTATION-TYPE BITS)
1353 ;                         '(LIST STRUC "REP=2" "REP=3"))
1354 ;                    (OR (NTH (LDB %%REGION-SPACE-TYPE BITS)
1355 ;                             '(FREE OLD NEW COPY STATIC STACK))
1356 ;                        (FORMAT NIL "TYPE=~O" (LDB %%REGION-SPACE-TYPE BITS)))
1357 ;                    (LDB %%REGION-READ-ONLY BITS)
1358 ;                    (LDB %%REGION-TEMPORARY BITS)
1359 ;                    (AND (PLUSP (LDB %%REGION-SWAPIN-QUANTUM BITS))
1360 ;                         (LDB %%REGION-SWAPIN-QUANTUM BITS))
1361 ;                     (LDB %%REGION-SCAVENGE-ENABLE BITS)
1362 ;                    (LDB %%REGION-SAFEGUARD BITS)
1363 ;                    (LDB %%REGION-NO-CONS BITS)
1364 ;                    (LDB %%REGION-EPHEMERAL BITS)
1365 ;                    (ZEROP (LDB %%REGION-LEVEL BITS))
1366 ;                    (LDB %%REGION-LEVEL BITS)))
1367 ;          (RETURN t)))))
1368 ;
1369 ;si:
1370 ;(defvar *bad-length* nil)
1371 ;si:
1372 ;(DEFUN RESET-AREA (AREA &OPTIONAL FREE-REGIONS &aux (reset t) nex)
1373 ;  (DO ((REGION (AREA-REGION-LIST AREA) NEXT-REGION)
1374 ;       (NEXT-REGION))
1375 ;      ((MINUSP REGION))
1376 ;    (SETQ NEXT-REGION (REGION-LIST-THREAD REGION))
1377 ;    (cond
1378 ;      ((> (region-free-pointer region) (region-length region))
1379 ;       (setq reset nil)
1380 ;       (push
1381 ;        (List
1382 ;          (time:print-current-time nil)
1383 ;          (region-free-pointer region)
1384 ;          (region-length region)
1385 ;          (describe-region region :stream nil)
1386 ;          (describe-region (%region-number (setq nex (+ (region-origin region)
1387 ;                                                        (region-free-pointer region))))
1388 ;                           :stream nil)
1389 ;          (describe-area 31 :stream nil)
1390 ;          region (region-bits region) area) *Bad-length*)
1391 ;       (beep)
1392 ;       (format tv:who-line-documentation-window
1393 ;              "~A****The region too long***" (describe-region region :stream nil))
1394 ;       (with-open-file
1395 ;        (str  "isaac:>wfs>bad-length" (list :out))
1396 ;        (describe-area (%area-number nex) :stream str)
1397 ;        (describe-area 31 :stream str)
1398 ;        (format str "~A" *bad-length*))
1399 ;       (beep))))
1400 ;  (cond (reset
1401 ;        (without-interrupts
1402 ;          (prog1
1403 ;        (DO ((REGION (AREA-REGION-LIST AREA) NEXT-REGION)
1404 ;             (NEXT-REGION))
1405 ;            ((MINUSP REGION))
1406 ;          (SETQ NEXT-REGION (REGION-LIST-THREAD REGION))
1407 ;          (GC-RESET-FREE-POINTER REGION 0)
1408 ;          (COND (FREE-REGIONS
1409 ;                 (SETF (AREA-REGION-LIST AREA) NEXT-REGION)
1410 ;                 (%GC-FREE-REGION REGION))))
1411 ;        (cond ((eq area macsyma:*temporary-polynomial-area*)
1412 ;               (Si:clear-cons-caches))))
1413 ;        ))))
1414
1415 (defun display (expr)
1416   (cond ((atom expr)(displa expr))
1417         ((or (affine-polynomialp expr)(rational-functionp expr))(sh expr))
1418         ((or (affine-polynomialp (car expr))(rational-functionp (car expr)))
1419          (loop for v in expr do (display v)))
1420         (t (displa expr))))
1421
1422
1423 (defun ldata-simplifications-write (ldata &key (open-g 1)
1424                                     (pathname "haskell:>wfs>answer.lisp") &aux answ)
1425   (setq answ (ldata-simplifications ldata :open-g open-g))
1426   (with-open-file (st pathname :direction :output)
1427     (let ((*standard-output* st) *print-radix*)
1428       (for-editor (des answ))
1429       (format st "~%(setq (answ (rerat '~A)))" answ))))