share/affine/sheafa.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 ;(defstruct ( bil (:type list) :named  (:conc-name go-))
  11 ;  head
  12 ;  tail)
  13
  14
  15 (defun update (ldata)
  16   (cond ((< (length ldata) 6)
  17          (nconc ldata (loop for i from (length ldata) below 6
  18                             collecting nil)))))
  19 ;
  20 ;(defstruct  (subscheme (:type list) :named )
  21 ;  s-var
  22 ;  ;;    for each open of the s-var need
  23 ;  ;;the defining equations
  24 ;  ideal-sheaf)
  25 ;
  26 ;
  27 ;(defstruct (open :named (:conc-name open-))
  28 ;    ;; a list of variables
  29 ;    variables
  30 ;    ;;a polynomial in prev variables to invert
  31 ;    inequality
  32 ;    ;;record the ordered pair of opens this came from
  33 ;    intersection)
  34 ;
  35 ;(defstruct (map :named (:conc-name map-))
  36 ;  ;;an open
  37 ;  domain
  38 ;  ;;an open
  39 ;  range
  40 ;  ;;a list of polynomial functions( maybe the same denominator)
  41 ;  substitutions
  42 ;  denom)
  43 ;
  44 ;(defstruct (morphism :named (:conc-name morphism-))
  45 ;  ;;two s-varieties
  46 ;  domain
  47 ;  range
  48 ;  ;;for each open in (sv-opens domain) a map  to one in
  49 ;  ;; (sv-opens range)
  50 ;  list-of-maps)
  51 ;;
  52 ;;(defun sv-find-intersections (s-var &key (set-them t) &aux tem answ)
  53 ;;  (setq answ (loop for v in (sv-opens s-var)
  54 ;;                 appending
  55 ;;                 (loop for w   in (sv-opens s-var)
  56 ;;                       when (not (eql v w))
  57 ;;                       collecting
  58 ;;                       (setq tem (make-open))
  59 ;;                       and
  60 ;;                       do
  61 ;;                       (setf (open-variables tem) (open-variables w))
  62 ;;                       ;;ought to set the inequality here but with out the map how??
  63 ;;                       (setf (open-inequality tem) (open-inequality v))
  64 ;;                       (setf (open-intersection tem )(list v w)))))
  65 ;;  (cond (set-them (setf (sv-intersections s-var) answ)))
  66 ;;  answ)
  67 ;(defun set_up_s_var(a-list birat-string &aux map answ amap subl subs intop )
  68 ;  "takes alternating list of coord inequal equations "
  69 ;  (setq answ (make-s-var))
  70 ;  (setf (sv-opens answ)
  71 ;       (loop for (coord inequal ident) on (cdr a-list) by #'cdddr
  72 ;             with op
  73 ;             collecting (setq op (make-open))
  74 ;             do (setf (open-variables op) (loop for v in (cdr coord)
  75 ;                                                collecting (car (st-rat v))))
  76 ;             (setf (open-inequality op) (st-rat inequal))
  77 ;             (setf (open-intersection op) ident)))
  78 ;  (loop for op1 in (sv-opens answ)
  79 ;        appending
  80 ;       (loop for op2 in (sv-opens answ)
  81 ;             for (coord inequal ident) on (cdr a-list) by #'cdddr
  82 ;             when (not (eql op1 op2))
  83 ;             collecting (setq intop (make-open))
  84 ;             and
  85 ;             do
  86 ;              (setf (open-variables intop) (open-variables op1))
  87 ;             (setq map ($find_birational_map (open-intersection op1)
  88 ;                                            (open-intersection op2)
  89 ;                                            birat-string
  90 ;                                            coord
  91 ;                                            ))
  92 ;             ;;;temporarily put it in the inequal slot
  93 ;             (setf (open-inequality intop) map)
  94 ;             (setf (open-intersection intop) (list op1 op2)))
  95 ;
  96 ;       into ints
  97 ;       finally (setf (sv-intersections answ) ints))
  98 ;  (loop for intop in (sv-intersections answ)
  99 ;       with op1 with op2
 100 ;       do (setq op1 (first (open-intersection intop)))
 101 ;        (show intop op1 op2)
 102 ;       (setq op2 (second (open-intersection intop)))
 103 ;       collecting (setq amap (make-map domain op1
 104 ;                                       range op2 ))
 105 ;       into glues
 106 ;       do
 107 ;       (setq map (open-inequality intop))
 108 ;       (setq subs
 109 ;             (loop for v in (cdr map)collecting
 110 ;                   (third v)))
 111 ;       (show subs)
 112 ;       (setq subl (loop for v in (open-variables op1)
 113 ;                        for vv in *vvv*
 114 ;                        collecting (cons vv  (get v 'disrep) )))
 115 ;       (setq subs (sublis subl subs))
 116 ;       (show subs)
 117 ;       (setq subs (mapcar 'new-rat subs))
 118 ;       (loop for v in subs
 119 ;             with denom = 1
 120 ;             do (setq denom (nplcm denom (denom v)))
 121 ;             finally (loop for v in subs collecting
 122 ;                            (ptimes (num v)(pquotient denom
 123 ;                                                          (denom v)))
 124 ;                            into sub
 125 ;                           finally (setf (map-substitutions amap)  sub)
 126 ;                           (show sub)
 127 ;                           (setf (map-denom amap) denom)))
 128 ;       (setf (open-inequality intop)
 129 ;             (apply-map amap (open-inequality op1)))
 130 ;       finally (setf (sv-glueing answ)glues))
 131 ;  answ)
 132 ;
 133 ;(defun find-map (list-maps dom rang)
 134 ;  (loop for ma in list-maps
 135 ;       when (and (eq (map-domain ma) dom)
 136 ;                 (eq (map-range ma) rang))
 137 ;       do (return ma)))
 138
 139 ;(defun check-compatibility (svar)
 140 ;  (loop for v in (sv-intersections svar)
 141 ;       with rev
 142 ;       do (setq rev (reverse (open-intersection v)))
 143 ;       (loop for w in (sv-intersections svar)
 144 ;             when (equal (rev (open-intersection w)))
 145 ;             do
 146 ;             (loop for va in (open-variables v) with var
 147 ;                   do
 148 ;             (setq var (list va 1 1))
 149 ;             (setq ma (find-map (sv-glueing svar)v w))
 150 ;             (setq rev-ma (find-map (sv-glueing svar)w v))
 151 ;       (setq answ (ratrquotient (apply-map rev-ma (num (apply-map ma var)))a
 152 ;                           (apply-map rev-ma (denom (apply-map ma var)))))
 153 ;               (cond ((not (and (equal (num answ) var)) (denom answ) 1))
 154 ;                      (merror "not the identity map")))))))
 155 ;
 156 ;(defun describe-s-var (svar)
 157 ;  (mapcar 'describe-open (sv-opens svar))
 158 ;  (mapcar 'describe (sv-glueing svar))
 159 ;  (mapcar 'describe (sv-intersections svar)))
 160 ;(setq hi
 161 ;#$ sublis([a1=z1,a2=z2,a0=1],[ [x1,x2],[x1-a1/a0,x2-a2/a0],[x1,x2],[x1-a0/a1,x2-a2/a1],[x1,x2],[x1-a0/a2,x2-a1/a2]])$)
 162 ;
 163 ;[[X1,X2],1,[X1-Z1,X2-Z2],[X1,X2],1,[X1-1/Z1,X2-Z2/Z1],[X1,X2],1,[X1-1/Z2,X2-Z1/Z2]]$
 164
 165 (defvar *xxx*
 166   (let (*print-radix*)
 167     (loop for i from 1 to 30 collecting (add-newvar (intern (format nil "$X~A" i))))))
 168
 169 ;(defun describe-open (open &aux tem)
 170 ;  (format t "~%It is ~D space coordinatized by ~A "
 171 ;         (length (setq tem (open-variables open))) tem)
 172 ;  (cond ((not (numberp  (open-inequality open)))
 173 ;        (format t "less the locus where ")
 174 ;        (sh (open-inequality open))(format t "vanishes.")))
 175 ;  (cond ((Open-intersection open)
 176 ;        (format t "~%It is an intersection of ~A" (open-intersection open)))))
 177 ;;
 178 ;(defun $find_birational_map (eqn1 eqn2 birat-sstring coord2 &aux sol1 birat-vars
 179 ;                            answ eqns)
 180 ;  (setq eqn1  ($numerator eqn1))
 181 ;  (setq eqn2 ($numerator eqn2))
 182 ;;  (mshow eqn1)
 183 ;  (setq birat-vars ($list_variables eqn1 birat-sstring))
 184 ;  (setq sol1 ($fast_linsolve ($expand eqn1) birat-vars))
 185 ;  (setq eqn2 (sublis *vforx* eqn2 ))
 186 ;  (setq coord2 (sublis *vforx* coord2))
 187 ;    (setq eqns ($sublis sol1 eqn2))
 188 ;;    (mshow eqns coord2)
 189 ;  (setq answ ($fast_Linsolve  eqns coord2))
 190 ;  (sublis answ (sublis *xforv* answ))
 191 ;  answ)
 192
 193 (defun zl-pairlis (a b)
 194        (loop for v in a for w in b
 195              collecting (cons v w)))
 196
 197 (defvar *vvv* (loop for i from 0 for u in *xxx* collecting (intern (format nil "$VV~A" i))))
 198 (defvar *vforx* (zl-pairlis *xxx* *vvv*))
 199 (defvar *xforv* (zl-pairlis *vvv* *xxx*))
 200
 201 ;(defun construct-birational-glueing-data (s-var ring-data birat-string)
 202 ;  "ring-data is a list one for each open of equations like x0=x1*z1,x2=x1*z2
 203 ; where the x0,x1,x2 are the open coordinates, and the zi are coordinates of the
 204 ; function field"
 205 ;  (check-arg $listp ring-data "macsyma list")
 206 ; (loop for int in (s-var intersections)
 207 ;       do (setq map
 208 ;               (make-map))
 209 ;       (int-opens (open-intersection int))
 210 ;       (setf (map-domain map) int)
 211 ;       (setf (map-range map)
 212 ;            (loop for w in int
 213 ;                             with rev = (reverse
 214 ;                                          (open-intersection int))
 215 ;                             when (equal (open-intersection w)
 216 ;                                         rev)
 217 ;                             do (return w)))
 218 ;       (setf (map-substitutions)
 219 ;            (progn
 220 ;              (loop for v in (cdr  ring-data)
 221 ;                  for w in (sv-opens s-var)
 222 ;                  when (equal w (first int-opens))
 223 ;                  do(setq eqn1 v)
 224 ;                  when (equal w (second int-opens))
 225 ;                  do (setq eqn2 w))
 226 ;                  (setq eqn2 (sublis *vforx* eqn2)
 227 ;                   (setq answ
 228 ;                        ($find_birational_map eqn1 eqn2 birat-string
 229 ;                                        (firstn (length eqn2) *vvv*)))
 230 ;                (setq substit
 231 ;                      (loop for w in answ collecting (new-rat (third w))))
 232 ;                (loop for w in substit
 233 ;                      with denom =1
 234 ;                      do (setq denom (plcm (denom w)))
 235 ;                      finally
 236 ;                      (setq substit (loop for w in substit
 237 ;                                          collecting (num (rattimes w
 238 ;                                                                    denom))))
 239 ;                      (setf (map-denom map) denom)
 240 ;                      (setf (map-substitutions map) substit)))))))
 241
 242 (defvar *p3* (make-s-var))
 243 ;
 244 ;(setf (sv-opens *p3*) (loop for i from 0 to 2
 245 ;                           collecting (setq tem (make-open))
 246 ;                           do
 247 ;                           (setf (open-variables tem)
 248 ;                                 (loop for j from 0 to 2
 249 ;                                       when (not (eql j i))
 250 ;                                       collecting (car (st-rat (nth j *xxx* )))))
 251 ;                           (setf (open-inequality tem) (st-rat (nth i *xxx*)))))
 252
 253
 254 ;;(setf (sv-intersections *p3*
 255 ;(defun apply-map (map poly &aux answ subl)
 256 ; (setq subl (pairlis (open-variables (map-range map)) (map-substitutions map)))
 257 ;  (setq answ  (psublis subl (map-denom map) poly))
 258 ;  (remove-common-factors answ (map-denom map)))
 259 ;
 260 ;(defun apply-map (map poly &aux answ subl)
 261 ; (setq subl (pairlis (open-variables (map-range map)) (map-substitutions map)))
 262 ;  (setq answ  (psublis subl (map-denom map) poly))
 263 ;  (remove-common-factors answ (map-denom map)))
 264
 265 (defun remove-common-factors (from-poly divisor &aux answ)
 266   (setq answ (pgcdcofacts from-poly divisor))
 267   (if (eql 1 (car answ))
 268       (second answ)
 269       (remove-common-factors (second answ) (first answ))))
 270
 271 (defun strings-search (keys sstring &aux tem)
 272   (loop for v in keys when (setq tem (search v sstring :test #'char-equal))
 273         do (return tem)))
 274
 275 ;; will replace the genvar with more readable symbols
 276 (defun replace-genvar(&rest strings &aux tem hi)
 277   (loop for v in *genvar*
 278         when (and (atom  (setq tem  (get v 'disrep)))
 279                   (eql (aref (string v) 0) #\G)
 280                   (cond (strings (strings-search strings (string tem)))
 281                         (t t)))
 282         do
 283         (setq hi (make-symbol (string-trim "$" tem)))
 284         (setf (symbol-value hi) (symbol-value v))
 285         (setf (symbol-plist hi) (copy-list (symbol-plist v)))
 286         and
 287         collecting hi into answ
 288         else
 289         collecting v into answ
 290         finally (setq *genvar* answ)))
 291
 292 (defun describe-components (comp &aux answ)
 293   (check-arg comp (and (listp comp) (eq (car comp) 'components)) "components")
 294 ;  (setq comp (cdr comp))
 295   (loop for v in (car (cdr comp))
 296         for i from 0
 297         do (format t "~%Component ~A occurs on opens ~A on ldata ~A" i
 298                    (loop for  w in v  when (cdr w) collecting (car w)
 299                          into tem
 300                          and
 301                          collecting (cdr w) into tem2
 302                          finally (setq answ (list tem tem2))
 303                          (return (car answ)))
 304                    (second answ)))
 305   (cond ((third comp) (format t "~%There were bad components so the above must be taken with a grain of salt.  They were ~A" (third comp)))
 306         (t (format t "~%There were no bad components" nil))))
 307
 308
 309
 310 (defun tableaux-lessp (a b)
 311   (let ((stringa (string a))
 312         (stringb (string b)))
 313     (loop for j from 3  below (string-length a)
 314           when (> (aref stringa j) (aref stringb j))
 315           do (return nil)
 316           finally (return t))))
 317
 318 (defun last-string (x)
 319   (aref x (1- (string-length x))))
 320
 321 (defremember list-tableaux ( row-length dimension )
 322   (case row-length
 323     (0  nil)
 324     (1 (loop for i from 1 to dimension
 325              collecting (list i)))
 326     (t (loop for j from 1 to dimension
 327              appending
 328              (loop for u in (list-tableaux (1- row-length ) dimension)
 329                    when (> j (car (last u)))
 330                    collecting (append u (list j)))))))
 331
 332 (defun $list_tableaux (length dimension)
 333   (cons '(mlist) (mapcar #'(lambda (x)
 334                             (apply #'$concat '$ta x))
 335                          (list-tableaux length dimension))))
 336
 337 (defun $sub_matrix (mat  rows-to-take cols-to-take &aux row)
 338   (loop for ii in rows-to-take
 339          do (setq row (nth ii mat))
 340          collecting
 341          (loop for i from 1
 342                 for v  in (cdr  row)
 343                 when (member i cols-to-take)
 344                 collecting v into new-row
 345                 finally (return (cons '(mlist) new-row)))
 346          into matrix
 347          finally (return (cons '($matrix simp) matrix))))
 348
 349 (defun $sub_matrix_columns (mat &rest cols-to-take)
 350   (loop for row in (cdr mat)
 351          collecting
 352          (loop for i from 1
 353                 for v  in (cdr  row)
 354                 when (member i cols-to-take)
 355                 collecting v into new-row
 356                 finally (return (cons '(mlist) new-row)))
 357          into matrix
 358          finally (return (cons '($matrix simp) matrix))))
 359
 360 (defun $plucker (vector-list)
 361   (check-arg vector-list '$listp "macsyma list")
 362   (let ((mat (cons '($matrix) (cdr vector-list))))
 363     (cons '(mlist)
 364           (loop for
 365                 u in
 366                 (list-tableaux ($length vector-list) ($length (second vector-list)))
 367                 collecting ($determinant (apply '$sub_matrix_columns mat u))))))
 368
 369 (defun $remove_number (llist)
 370   (loop for v in  llist when (not (numberp v))
 371         collecting v))
 372
 373 (defmacro from-case (a-list)
 374   `(and (listp ,a-list) (cdr (member '$case ,a-list :test #'eq))))
 375
 376 (defvar $solution_tree '(mlist))
 377 (defun $sort_solution_tree ()
 378  (setq $solution_tree
 379        (cons '(mlist)(sort (cdr $solution_tree)
 380                            #'(lambda (x y )
 381                                (let ((fromx (from-case x))
 382                                      (fromy (from-case y)))
 383                                  (cond ((null fromx) t)
 384                                        (t
 385
 386                                 ($monomial_alphalessp fromx
 387                                                       fromy)))))))))
 388 (defmacro type? (a form &optional descrip &aux (ctrl-string "not a ~A"))
 389   (cond (descrip nil)
 390         (t (setq descrip form) (setq ctrl-string "Doesn't satisfy ~A")))
 391   (cond ((functionp form)
 392          `(cond ( (, form ,a) nil)
 393                 (t (merror (format nil ,ctrl-string ,descrip)))))
 394         (t `(cond (,form nil)
 395                   (t        (t (merror (format nil ,ctrl-string ,descrip))))))))
 396
 397 (defvar $nonzero_factors nil)
 398
 399 (defun $dichotomy (polynomials dich &optional add-to-tree  &aux answer)
 400   (check-arg polynomials $listp  "macsyma list")
 401   (cond ($nonzero_factors
 402          (format t "~%Assuming the following factors are nonzero and cancelling")
 403          (displa $nonzero_factors)
 404          (setq polynomials ($Eliminate_nonzero_factors polynomials))))
 405 ;  (cond (add-to-tree (type? $solution_tree (and ($listp $solution_tree)
 406 ;                                               (or (null (second $solution_tree))
 407 ;                                                   ($listp (second $solution_tree)))))))
 408
 409   (cond (($listp dich)(setq dich(second dich))))
 410   (cond ((numberp dich)(setq dich (nth  dich polynomials))))
 411   (setq dich
 412         (loop for v in (cdr ($list_irreducible_factors dich))
 413               when (not (numberp v))
 414               collecting v))
 415   (format t "We have the dichotomy among..")
 416     (displa (cons '(mlist) dich))
 417   (loop for w on dich
 418         for v in dich
 419         collecting
 420         (let (($nonzero_factors (append (cons '(mlist) (cdr  $nonzero_factors))
 421                                         (cdr w))))
 422           ( $factor($eliminate_nonzero_factors ($simplify_relative  polynomials v))))
 423         into tem
 424         finally (setq answer (cons '(mlist) tem)))
 425   (cond (add-to-tree ($add_to_solution_tree answer)))
 426   answer)
 427
 428 (defun constant-term (poly)
 429   (cond ((atom poly) poly)
 430         (t (constant-term
 431              (let ((leng (length poly)))
 432                (cond ((eql (nth (- leng 2) poly) 0)
 433                       (constant-term (nth (1- leng) poly) ))
 434                      (t 0)))))))
 435
 436
 437 (defun gone-prepared (poly &key (inequal 1) linear-variables &aux tem)
 438    (cond ((null linear-variables) (setq linear-variables
 439                                       (degree-one-variables poly))))
 440   (cond ((numberp  inequal )
 441          (one-prepared poly :linear-variables linear-variables))
 442         (t
 443          (cond  ((atom poly) nil)
 444                 (t (loop for v in linear-variables
 445                          when (eql (pdegree poly v) 1)
 446                          do (setq tem  (zero-sublis poly v))
 447                          (cond ((may-invertp tem inequal)
 448                                       (return v)))))))))
 449
 450 (defun degree-one-variables (poly)
 451   (loop for v in (list-variables poly)
 452         when (eql (pdegree poly v) 1)
 453         collecting v ))
 454
 455 (defun one-prepared (poly &key linear-variables &aux tem)
 456  (cond ((null linear-variables) (setq linear-variables
 457                                       (degree-one-variables poly))))
 458   (cond  ((atom poly) nil)
 459          ((zerop (constant-term poly)) nil)
 460          (t (loop for v in  linear-variables
 461                   when (eql (pdegree poly v) 1)
 462                   do (setq tem  (zero-sublis  poly v))
 463                   (cond ((numberp tem) (return v)))))))
 464
 465 ;;the following was gm-prepared by x4,x3 !!
 466 ;;2*X3*X6+2*X3^2*X4-1
 467 ;;(SETQ F (RERAT (QUOTE (X6 1 (X3 1 2) 0 (X4 1 (X3 2 2) 0 -1)))))
 468 ;(defun gm-prepared (poly &key  m (inequal 1) linear-variables &aux lins tem answ)
 469 ;  (cond (m
 470 ;  (cond ((atom poly ) nil)
 471 ;       ((eql 1 m)(cond ((setq tem (gone-prepared poly inequal :linear-variables linear-variables)))
 472 ;                       (list tem))))
 473 ;       (t (loop for v in (list-variables poly)
 474 ;                when
 475 ;                (eql (pdegree poly v) 1)
 476 ;                do (setq tem (zero-sublis poly v))
 477 ;                (cond ((setq answ (gm-prepared tem :m (1- m) :inequal inequal))
 478 ;                       (return (cons v answ))))))))
 479 ;       (t (setq lins (loop for v in (list-variables poly)
 480 ;                           when (eql (pdegree poly v) 1)
 481 ;                           count 1))
 482 ;          (loop for i from 1 to lins
 483 ;                when (setq answ (gm-prepared poly :m i :inequal inequal))
 484 ;                do (return answ))))
 485 ;;;failed to check gm-prepared on this.
 486 ;(des-editor badg)1/10/85 18:02:59
 487 ;(AA1*AA2-AA1)*BB1*CC6+(AA1^2*AA2-AA1^2)*CC5+BB1*CC3+AA1*CC2+((AA1*AA2^2-AA1*AA2)*BB6+AA2*BB3-1)
 488 ;*CC1+(AA1^2*AA2-AA1^2)*BB1*BB6+(AA1^3*AA2-AA1^3)*BB5+AA1*BB1*BB3+AA1^2*BB2-AA1*BB1
 489 ;
 490 ;(SETQ BADG (RERAT (QUOTE (CC6 1 (BB1 1 (AA2 1 (AA1 1 1) 0 (AA1 1 -1))) 0 (CC5 1 (AA2 1 (AA1 2 1) 0 (AA1 2 -1)) 0 (CC3 1 (BB1 1 1) 0 (CC2 1 (AA1 1 1) 0 (CC1 1 (BB6 1 (AA2 2 (AA1 1 1) 1 (AA1 1 -1)) 0 (BB3 1 (AA2 1 1) 0 -1)) 0 (BB6 1 (BB1 1 (AA2 1 (AA1 2 1) 0 (AA1 2 -1))) 0 (BB5 1 (AA2 1 (AA1 3 1) 0 (AA1 3 -1)) 0 (BB3 1 (BB1 1 (AA1 1 1)) 0 (BB2 1 (AA1 2 1) 0 (BB1 1 (AA1 1 -1))))))))))))))
 491 ;NIL
 492
 493
 494
 495 (defun gm-all-prepared (poly &key m (inequal 1) linear-variables &aux tem)
 496   (cond ((null linear-variables) (setq linear-variables
 497                                        (degree-one-variables poly))))
 498   (cond (m
 499          (cond ((eql m 0)(cond ((may-invertp poly inequal)(list 'ok))
 500                                (t nil)))
 501                (t
 502                 (loop for v in (list-variables poly)
 503                       when (and  (member v linear-variables :test #'eq)
 504                                  (setq tem (gm-all-prepared
 505                                              (zero-sublis poly v)
 506                                              :linear-variables linear-variables
 507                                              :m (1- m)
 508                                              :inequal   inequal)))
 509
 510                       appending (loop for ww in tem collecting (cons v ww))))))
 511         (t (setq linear-variables (degree-one-variables poly))
 512            (loop for i from 1 to (length linear-variables)
 513                  when (setq tem (gm-all-prepared poly
 514                                                  :linear-variables linear-variables
 515                                                  :m i
 516                                                  :inequal       inequal))
 517                  do (return tem)))))
 518
 519 (defvar *maximum-size-for-m-prepared* 1)
 520
 521 (defun gm-prepared (poly &key  m (inequal 1) linear-variables &aux lins tem answ)
 522      (cond ((null linear-variables) (setq linear-variables
 523                                       (degree-one-variables poly))))
 524   (cond (m
 525          (cond ((atom poly ) nil)
 526                ((eql 1 m)(cond ((setq tem (gone-prepared poly :inequal inequal  :linear-variables
 527                                                         linear-variables))
 528                                (list tem))))
 529                (t (loop for v in (list-variables poly)
 530                         when
 531                           (member v linear-variables :test #'eq)
 532                           do (setq tem (zero-sublis poly v))
 533                              (cond ((setq answ (gm-prepared tem :m (1- m) :inequal inequal
 534                                                             :linear-variables linear-variables))
 535                                     (return (cons v answ))))))))
 536
 537         (t (setq lins (loop for v in linear-variables
 538                             when (eql (pdegree poly v) 1)
 539                               count 1))
 540            (loop for i from 1 to (min lins  *maximum-size-for-m-prepared* )
 541                  when (setq answ (gm-prepared poly :m i :inequal inequal :linear-variables
 542                                               linear-variables))
 543                    do (return answ)))))
 544
 545
 546 (defun m-prepared (poly &key m linear-variables &aux  tem answ)
 547   (cond ((null linear-variables) (setq linear-variables
 548                                        (degree-one-variables poly))))
 549   (cond (m
 550          (cond ((atom poly ) nil)
 551                ((zerop (constant-term poly)) nil)
 552                ((eql 1 m)(cond ((setq tem (one-prepared poly))
 553                                (list tem))))
 554                (t (loop for v in linear-variables
 555                         do (setq tem  (zero-sublis poly v))
 556                         (cond ((setq answ (m-prepared tem :m (1- m)))
 557                                (return (cons v answ))))))))
 558         (t
 559          (loop for i from 1 to (length linear-variables)
 560                when (setq answ (m-prepared poly :m i))
 561                do (return answ)))))
 562
 563 (defun non-constant-factors (poly &optional invert &aux tem (genvar *genvar*))
 564   (setq tem (npfactor poly))
 565   (loop for (v deg) on tem by #'cddr
 566      when (and (not (numberp v))
 567                (or (null invert)(null (may-invertp v invert))))
 568      collecting v
 569      and
 570      collecting deg))
 571
 572 (defvar *factored-list* nil)
 573 (defvar *all-factors* nil)
 574
 575 (defun gen-ptimes (&rest l)
 576   (cond ((null l) 1)
 577         ((eql (length l) 1) (car l))
 578         (t (ptimes (car l) (apply 'gen-ptimes (cdr l))))))
 579
 580
 581 (defun find-good-dichotomy (ldata &key ( open-g 1) &aux tem list-factors answ
 582                             (list-polys (ldata-eqns ldata))
 583                             (gg (ldata-inequality ldata)))
 584   (setq gg (nplcm gg open-g))
 585   (cond ((member nil list-polys :test #'eq) (break 'here)))
 586   (setq list-factors
 587         (loop for v in list-polys
 588               ;;only collect non trivial factors
 589               when (> (length  (setq tem      (non-constant-factors v gg)))
 590                       2)
 591               collecting  tem into multiple-factors
 592               ;;when have a constant factor
 593               when (null tem)
 594               do (show tem) (setf (ldata-eqns ldata) '(1))
 595               (setf (ldata-usedup ldata ) 1)
 596               (return 'done-unit-ideal)
 597               collecting tem into all-factors
 598               collecting
 599               (apply 'gen-ptimes (loop for (fac deg) on tem by #'cddr
 600 ;                                      when (not (may-invertp fac gg))
 601                                        collecting fac))
 602               into terms
 603
 604               finally
 605               (cond ((not (member 1 terms))
 606                      (setf (ldata-eqns ldata) terms))
 607                     (t (setf (ldata-usedup ldata) 1)
 608                        (setf (ldata-eqns ldata) '(1))))
 609               (setq *all-factors* all-factors)
 610               (return multiple-factors)))
 611
 612   (cond ((member '(nil 1) *all-factors* :test #'equal) (break 'hii)))
 613   (cond
 614     ((eq list-factors 'done-unit-ideal) nil)
 615     (t
 616      (setq *factored-list* list-factors)
 617      (setq list-factors
 618            (sort list-factors #'(lambda (u v)(cond ((atom u) t)
 619                                                    ((atom v) nil)
 620                                                    (t (< (length u) (length v)))))))
 621      (cond
 622        ((setq answ
 623               ;; grabs the shortest product a^2*b^3 where a,b are variables
 624               (loop for v in list-factors
 625                     when (loop for (fac deg) on v by #'cddr
 626                                when (not (or (atom fac)
 627                                              (< (length fac) 4)))
 628                                do
 629                                (return nil)
 630                                finally (return t))
 631                     do (return v))))
 632
 633        ((setq answ (loop named kay
 634                          for i from 1 to (length (list-variables list-polys))
 635                          do
 636                          (loop for v in list-factors
 637                                when (loop for (fac deg) on v by #'cddr
 638                                          when (not (or (atom fac)
 639                                                        (< (length fac) 4)
 640                                                        (gm-prepared fac :m i :inequal gg)))
 641                                          do
 642                                          (return nil)
 643                                          finally (return t))
 644                                do (return-from kay v)))))
 645        ;;grab the shortest one with a variable factor
 646        ((setq answ
 647               (loop named pat
 648                     for
 649                     v in list-factors
 650                     do (loop for (fac deg) on v by #'cddr
 651                              when  (< (length fac ) 4)
 652                              do (return-from pat v)))))
 653        ((setq answ
 654               (loop named sue
 655                     for
 656                     v in list-factors
 657                     do (loop for (fac deg) on v by #'cddr
 658                              when  (or (one-prepared fac)
 659                                        (any-linearp fac
 660                                                     (ldata-inequality ldata)))
 661                              do (return-from sue v))))))
 662      answ)))
 663
 664 (defun order-dichotomy ( dich  &aux mult)
 665   (setq dich (loop for (v deg) on dich by #'cddr
 666                    collecting v))
 667   (setq mult
 668         (loop for v in dich
 669               collecting
 670               (loop for u in *all-factors*
 671                     when (member v u :test #'equal)
 672                     count 1 )))
 673   (setq dich
 674         (loop for i downfrom (length *all-factors*) to 0
 675               when  (member i mult)
 676               appending
 677               (loop for v in mult
 678                     for u in dich
 679                     when (eq v i)
 680                     collecting u))))
 681
 682
 683
 684 ;;LDATA : ((F1 F2 F3 ..) G)
 685 ;;LDATA MEANS HAVE FI=0 AND G NOT ZERO
 686 ;;SIMPLIFICATION STRATEGY FOR A LIST OF LDATUM
 687 ;;IF U^2+5*G*X IS IN THE LIST OF EQUATIONS ELIMINATE X FROM ALL OF THE EQUATIONS
 688 ;;(EXCEPT ITSELF)
 689 ;;IF 5+U^2*Y*X APPEARS SOLVE FOR X ELIMINATING IT FROM THE OTHERS AND SAVING THAT
 690 ;;EQUATION BUT MODIFYING THE INEQUALITY BY  G:LCM (G,U^2*Y)
 691 ;;MUST ELIMINATE FACTORS OF G FROM THE Fi
 692 ;;If 5+u^2*x+v^2*y occurs in the fi we create two ldata in place of the one ldata
 693 ;;   must solve for x=.. and add u^2 to the G making one ldata
 694 ;;   must solve for y=.. and add v^2 to the G making another ldata
 695 ;;Also if 5+u*v is an fi can set G: (lcm(g, u*v) and can eliminate u
 696 ;;from all equations
 697 ;;
 698 ;
 699 ;(defun dichotomy (list-eqns inequality)
 700 ;  (setq dich (find-good-dichotomy list-eqns))
 701 ;  (setq dich  (order-dichotomy dich))
 702 ;  (loop for v in dich
 703 ;;list of simplification types and their tests: where g is the inequal
 704 ;;and f is the equation to use in the replacement.
 705 ;;;   f                  test            replacement of h
 706 ;;1 u^2+5*g*x  (poly-linearp f x g)         (gen-prem h f x)
 707 ;;2 g1+u*x     (gone-prepared f g) (progn (setq g (lcm g u) ) (gen-prem h f x))
 708 ;;3 u1+u2*x+g*x^2  (invertible-leading-coefficient f x g)   (gen-prem h f x))
 709
 710 (defun any-linearp (f g &key variables-to-exclude among-variables)
 711   (cond ((null among-variables)(setq among-variables (list-variables f))))
 712   (loop for v in among-variables
 713
 714         when (and (not (member v variables-to-exclude :test #'eq)) (poly-linearp f v g))
 715         do (return v)))
 716
 717 (defvar *clear-above* nil)
 718
 719 (defun any-invertible-leading-coefficient
 720        (f g &aux deg (varl (list-variables f)))
 721   (loop for v in varl
 722          when  (may-invertp (list v 1 1) g)
 723          do (setq varl (delete v varl :test #'equal)))
 724   (loop for v in varl
 725          when (and (or (null *clear-above*) (not (<= (loop for w in *clear-above*
 726                                                             when (eq (car w) v)
 727                                                             minimize (cdr w))
 728                                                      (setq deg (pdegree f v)))))
 729                    (may-invertp (pcoeff f (list v deg 1)) g))
 730          do (return v)))
 731
 732 (defun invertible-leading-coefficient (poly var inequal)
 733   (may-invertp (pcoeff poly (list var (pdegree poly var) 1)) inequal))
 734
 735 (defun plength-order (u v)
 736   (cond ((atom u) t)
 737         ((atom v) nil)
 738         (t (< (length u) (length v)))))
 739
 740 ;(defun replace-functions (list-to-replace f var &aux tem)
 741 ;  (loop for h in list-to-replace
 742 ;       when (not (eql f h))
 743 ;       when (not (pzerop (setq tem
 744 ;                           (square-free (gen-prem h f var)))))
 745 ;       collecting
 746 ;       tem))
 747 ;
 748 ;(defun replace-functions (list-to-replace f var &AUX original &aux tem)
 749 ;  (loop for h in list-to-replace
 750 ;       when (not (eql f h))
 751 ;       when (not (pzerop (setq tem
 752 ;                               (gen-prem h f var))))
 753 ;       when (not (eq h tem))
 754 ;       collecting (square-free tem)
 755 ;       else collecting h))
 756
 757 (defun replace-functions (list-to-replace f var &key  general-leading-cof (invertible-g 1) &aux c-reqd remaind tem)
 758   (cond ((null general-leading-cof)
 759          (loop for h in list-to-replace
 760             when (not (eql f h))
 761             when (not (pzerop (setq tem
 762                                     (gen-prem h f var))))
 763             when (not (eq h tem))
 764             collecting (square-free tem)
 765             else collecting h))
 766         (t
 767          (loop for h in list-to-replace
 768             when (eq f h) do (setq h nil)
 769             when h
 770             do
 771             (multiple-value
 772              (remaind c-reqd)
 773              (gen-prem h f var))
 774             when  h
 775             when (eq h remaind) collecting h
 776             else collecting (square-free remaind)
 777             ;;collect h if the c-reqd is not a unit.
 778             when  h
 779             when (not (may-invertp c-reqd invertible-g))
 780             collecting h))))
 781
 782 (defun any-gm-prepared (poly gg)
 783   (gm-prepared poly :inequal gg))
 784
 785 (defvar *inside-simplify-svar-ldata* nil)
 786 (defvar   *stop-simplify* nil)
 787
 788
 789
 790 (defun check-for-gm-prepared (list-eqns gg &aux tem)
 791   (loop for eqn in list-eqns
 792         when (and (not (any-linearp eqn gg))
 793                   (setq tem (gm-prepared eqn :inequal gg)))
 794         do (setq *stop-simplify* (list eqn (length tem)))
 795         (return 'done))
 796   *stop-simplify*)
 797
 798 ;;if we have a gm-prepared poly then we can simplify the
 799 ;;ldata by moving to m+1 opens where we are able to eliminate
 800 ;;variables:
 801 ;;eg.if  x*u+g=0  is an equation in ldata then we can assume u is
 802 ;;generically nonzero so we can on the open set u not zero solve
 803 ;;for x=-g/u.  The u=0 locus does not meet the x*u+g=0 locus so as far
 804 ;;as calculation of components goes, we need not keep track.
 805 ;;If however we wish to make a change of coordinates and look at some
 806 ;;other data then we must cover ourselves with two open sets:
 807 ;;u not zero and x*u+g not zero.
 808 ;;criterion for f going in used-up should be that there is
 809 ;;no possibility of using f any more.  This will be the case
 810 ;;for example if there is a variable in f and no other polynomials in the ldata
 811 ;;have that variable occurring. This happens after a linear is used eg. Otherwise
 812 ;;there is the possibility of dividing one by the other or taking a resultant.
 813
 814 ;;note the possibility of x1+x2+x3*x4 = 0 being an equation.  It
 815 ;;might be that the replacement of x1 in the other equations is preferable
 816 ;;to the replacement of x2 or vice versa.eg in the following x7 replacement is better
 817 ;;since it leads to splitting into two components.
 818 ;;2*X6*X7-2*X6^2+X4*X5
 819 ;;X7-X6-X1*X5
 820
 821 (defmacro check-containments (ldata lis-ldat &optional gg)
 822   `(cond (error-check-containments
 823          (check-component-containment ,ldata ,lis-ldat ,gg))))
 824 (Defun Ldata-Simplifications (ldata &key (open-g 1) error-check-containments
 825                               recursive-p
 826                               &aux tem *clear-above* simplif answ unused stop-simplify
 827                               used-up var (changed t)  orig-ldata)
 828   (declare (special *in-linear-dich*))
 829   (check-arg ldata (eq (car ldata) 'ldata) "not an ldata")
 830   (setq orig-ldata ldata)
 831   (setq ldata (copy-list ldata))
 832
 833   (let ((fns (ldata-eqns ldata))
 834         (gg (ldata-inequality ldata)))
 835     (setq gg (nplcm open-g gg))
 836     ;;note that using this gg we may produce I1 with I ^P I1R[gg-1] but
 837     ;;unless I1 is prime we will not have I ^P I1
 838     (setq fns  (loop for f in fns
 839                      when (not (and (numberp f) (zerop f)))
 840                      do (setq tem (square-free (remove-common-factors f gg)))
 841                         and
 842                      when (numberp tem) do (setf (ldata-usedup ldata) 1)
 843                      (return (setq fns '(1)))
 844                      else
 845                      collecting  tem))
 846     (setq used-up (subseq fns 0 (ldata-usedup ldata)))
 847     (setq fns (nthcdr (ldata-usedup ldata) fns))
 848     (setq fns  (sort (copy-list fns) 'plength-order))
 849     (setq vars (list-variables fns))
 850     (loop while changed
 851           do (setq changed nil)
 852           (loop named sue for test in '(any-linearp
 853                                          any-invertible-leading-coefficient
 854 ;                                        any-irreducible
 855                                          )
 856                 do
 857
 858                 (loop for f in fns
 859                       when (setq var (funcall test f gg))
 860                       do
 861 ;                      (cond ((eq test 'any-irreducible)
 862 ;
 863 ;                            (setq f var) (setq var (find-variable-with-simple-lc
 864 ;                                                     f 1 1))
 865 ;                            (sh f) (show  var)))
 866                       (setq fns
 867                             (replace-functions fns f var))
 868                       (setq used-up (replace-functions used-up f var))
 869                       (cond ((or (eq test 'any-linearp)
 870                                  (variable-doesnt-occur var fns used-up))
 871                              (push f used-up)
 872                              )
 873                             ;;put f last so that the lower degree ones will come first.
 874                             (t (setq fns (nconc fns (list f)))))
 875
 876                       (cond ((eq test 'any-invertible-leading-coefficient)
 877                              (push (cons var (pdegree f var)) *clear-above*)
 878                              (show *clear-above*)))
 879
 880                       (setq changed t)
 881                       (return-from sue))))
 882
 883     ;;if contain number 'done else try to fix leading cofs
 884     (check-containments orig-ldata (list(make-ldata :eqns (zl-UNION fns used-up)
 885                                                     :inequality (ldata-inequality
 886                                                                  orig-ldata))))
 887
 888     (cond ((loop for f in (setq unused  (zl-union (delete 0 fns)))
 889                  when (numberp f) do (setq used-up '(1) unused nil)
 890                 (return 'done)))
 891           (t
 892 ;          nil))
 893            (setq simplif
 894 ;                unused)
 895                  (simp-lead unused :open-g gg))
 896            (cond ((eq simplif 'try-dichotomy) nil)
 897                  (t (setq unused  (zl-UNION simplif))))))
 898     (setf (ldata-inequality ldata) gg)
 899     (setf (ldata-eqns ldata)
 900           (setq fns  (append  used-up unused)))
 901     (setf (ldata-usedup ldata) (length used-up))
 902     ;;the following may make  *stop-simplify* and ask to refine.
 903     ;;this check is done in MAKE-DICHOTOMY (check-for-gm-prepared  unused open-g)
 904     (setq stop-simplify *stop-simplify*)
 905     (check-containments orig-ldata (list ldata) (ldata-inequality ldata))
 906     (cond ((null *stop-simplify*)  (setq answ (MAKE-DICHOTOMY ldata :open-g open-g)))
 907           (t (setq answ (list ldata))))
 908     ;;should make this so it won't check for redundant ldata more than once.
 909     (cond ((> (length answ ) 1)
 910            (setq answ (delete-redundant-ldata answ :gg open-g)))))
 911
 912   (check-containments orig-ldata answ)
 913   (cond ((and (null *stop-simplify*)
 914               ;;this makes the divide dichotomy only apply after no more prod. dichot.
 915               (eql (length answ) 1))
 916          (setq answ
 917                (divide-dichotomy (car answ) :open-g open-g))))
 918   (check-containments orig-ldata answ)
 919 ;  (cond ((not (equal (length answ) 1))  (mshow answ)))
 920   (cond ((and (null *stop-simplify*)
 921               (eql (length answ) 1))
 922          (setq answ (try-factor-irreducible-ldata (car answ) open-g))))
 923    (cond ((and (null *stop-simplify*) (not recursive-p)
 924
 925           (let ((*in-linear-dich* t))
 926             (setq answ
 927                   (loop for ld in answ
 928                         appending (linear-dichotomy ld :open-g open-g )))))))
 929 ;   (cond ((and (null *stop-simplify*)(not (variable-boundp *in-linear-dich*))
 930 ;              (not recursive-p)
 931 ;              (eql (length answ) 1))
 932 ;         (let ((*in-linear-dich* t))
 933 ;           (setq answ (linear-dichotomy (car answ) :open-g open-g )))))
 934 ;    (format t "~%Verifying the ~A component contain the original" (length answ))
 935 ;    (check-components-contain-original ldata answ)
 936   (check-containments orig-ldata answ)
 937 ;  (cond ((null recursive-p)
 938 ;        (setq answ (loop for v in answ
 939 ;                         appending
 940 ;                      (jacobian-dichotomy ldata :open-g open-g)))))
 941   (setq answ (delete-redundant-ldata answ :gg open-g))
 942   (check-containments orig-ldata answ)
 943   answ)
 944
 945
 946 (defun any-irreducible (f gg &aux fac)
 947   (setq fac (non-constant-factors f gg))
 948   (cond ((> (length fac) 2)nil)
 949         (t (car fac))))
 950
 951 (defun leading-cof (poly var &aux (deg (pdegree poly var)))
 952   (cond ((zerop deg) 0)
 953         (t (pcoeff poly (list var (pdegree poly var)  1)))))
 954
 955 (defun find-variable-with-simple-lc (f  fns gg)
 956   fns gg
 957   (let ((varl (list-variables f)))
 958     (loop for v in varl collecting (pcomplexity (leading-cof f v)) into tem
 959           finally(return (nth  (find-position-in-list (apply 'min tem) tem)
 960                                varl)))))
 961
 962 (defun ldata-unused (ld)
 963   (nthcdr (ldata-usedup ld)(ldata-eqns ld)))
 964 (defvar *refine-opens* t)
 965
 966 ;;not sure about the *refine-opens* = nil mode working.
 967 #+old
 968 (defun MAKE-DICHOTOMY (ldata &key (open-g 1)&aux int-open-g all-facs stop-simplify  eqns-modv ld  answ gg dich tem lin-dich)
 969   "If stop-simplify is true then it only works if have linear dichotomy.  It returns
 970   a list of  ldata "
 971   (cond (*stop-simplify* (merror "how did *stop-simplify* get  here")))
 972   ;;(setq *stop-simplify* nil)
 973   (setq dich (find-good-dichotomy  ldata))
 974   (setq all-facs *all-factors*)
 975   (setq dich (order-dichotomy dich))
 976   (show dich)
 977   (setq gg (nplcm open-g (ldata-inequality ldata)))
 978   (cond ((null *refine-opens*) (setq int-open-g gg) )
 979         (t (setq int-open-g open-g)))
 980   (setq lin-dich
 981         (loop for v in dich when (not (any-linearp v gg)) do (return nil)
 982               finally (return (and dich t))))
 983   (show lin-dich)
 984 ;; I think the gm-prepared business should all be done after.
 985 ;; this may not be true.  The simplifications from finding a gm-prepared
 986 ;; and performing the elimination of variables might be necessary.
 987   (cond ((null lin-dich)
 988          (check-for-gm-prepared (ldata-eqns ldata) int-open-g)
 989          (show *stop-simplify*)))
 990   (setq stop-simplify *stop-simplify*)
 991   (cond ((and (null *refine-opens*) *stop-simplify* )
 992          (multiple-value-bind
 993            (lds ggs) (ldata-refinement ldata (car *stop-simplify*)
 994                                        (second *stop-simplify*) :inequality
 995                                        int-open-g)
 996            (setq *stop-simplify* nil)
 997            ;;this does not take into account the fact we should be contracting
 998            ;;the ideal back to the current open..
 999            (setq answ (loop for ld in lds
1000                             for ggi in ggs
1001                             do (setq tem (ldata-simplifications
1002                                         ld
1003                                         :open-g (sftimes int-open-g ggi)
1004                                         :recursive-p t))
1005                             appending
1006                             (loop for v in tem do
1007                                  (zl-copy-structure v ldata- open-inequality
1008                                                  (sftimes ggi (ldata-open-inequality
1009                                                                 v))))))))
1010
1011 ;        (setq answ (ldata-refinement ldata (car *stop-simplify*)
1012 ;                                     (second *stop-simplify*) :inequality int-open-g))
1013 ;;       (mshow answ)  (format t "**Is the refinement ")
1014 ;        (setq *stop-simplify* nil)
1015 ;        (setq answ (loop for v in answ
1016 ;                                         appending (ldata-simplifications v :open-g open-g :recursive-p t)))
1017
1018   ;;priority 1 Linear-dichotomy
1019   ;;         2 gm-prepared equation
1020   ;;         3 any dichotomy
1021   ;;proceed with dich if dich is linear or if found no gm-prepared
1022         ((or lin-dich (null *stop-simplify*))
1023          (cond (dich
1024                 (loop for v in dich
1025                       with so-far = 1
1026                       appending
1027                         (progn
1028                           (setq eqns-modv
1029                                 (loop for facs in all-facs
1030                                       when (not (member v facs :test #'equal))
1031                                         collecting
1032                                           (apply #'gen-ptimes (loop for ter in facs by #'cddr
1033                                                                    when (not
1034                                                                           (may-invertp
1035                                                                             ter so-far))
1036                                                                      collecting ter))))
1037                           (cond ((member nil eqns-modv :test #'eq) (merror "nil should not be here")))
1038                           (cond ((eq v nil) (merror "nil should not be here")))
1039                           (setq ld (make-ldata))
1040                           (setf (ldata-eqns ld) (cons v eqns-modv))
1041                           (setf (ldata-inequality ld)(nplcm gg so-far))
1042                           (setf so-far (nplcm so-far v))
1043                           (cond ((Null *stop-simplify*)
1044                                  (LDATA-SIMPLIFICATIONS ld :open-g open-g
1045                                                         :recursive-p t))
1046                                 (t (list  ld))))
1047                         into list-of-ld
1048                       finally
1049                         (cond ((null list-of-ld)
1050                                (setq answ (list (make-ldata :eqns '(1) :inequality 1))))
1051                               (t (setq answ list-of-ld)))))
1052                (t (setq answ (list ldata))))))
1053
1054   (cond ((null answ) (setq answ (list ldata))))
1055 ;;  (cond ((null answ) (setq answ (list (make-ldata :eqns '(1) :inequality 1)))))
1056 ;  (setq answ (delete-redundant-ldata answ)))
1057   answ)
1058 ;;not sure about the *refine-opens* = nil mode working.
1059 (defun MAKE-DICHOTOMY (ldata &key (open-g 1)&aux int-open-g all-facs stop-simplify  eqns-modv ld  answ gg dich lin-dich)
1060   "If stop-simplify is true then it only works if have linear dichotomy.  It returns
1061   a list of  ldata "
1062   (cond (*stop-simplify* (merror "how did *stop-simplify* get  here")))
1063   ;;(setq *stop-simplify* nil)
1064   (setq dich (find-good-dichotomy  ldata))
1065   (setq all-facs *all-factors*)
1066   (setq dich (order-dichotomy dich))
1067   (show dich)
1068   (setq gg (nplcm open-g (ldata-inequality ldata)))
1069   (cond ((null *refine-opens*) (setq int-open-g gg) )
1070         (t (setq int-open-g open-g)))
1071   (setq lin-dich
1072         (loop for v in dich when (not (any-linearp v gg)) do (return nil)
1073               finally (return (and dich t))))
1074   (show lin-dich)
1075 ;; I think the gm-prepared business should all be done after.
1076 ;; this may not be true.  The simplifications from finding a gm-prepared
1077 ;; and performing the elimination of variables might be necessary.
1078   (cond ((null lin-dich)
1079          (check-for-gm-prepared (ldata-eqns ldata) int-open-g)
1080          (show *stop-simplify*)))
1081   (setq stop-simplify *stop-simplify*)
1082   (cond ((and (null *refine-opens*) *stop-simplify* )
1083          (multiple-value-bind
1084            (lds ggs) (ldata-refinement ldata (car *stop-simplify*)
1085                                        (second *stop-simplify*) :inequality
1086                                        int-open-g)
1087            (setq *stop-simplify* nil)
1088            ;;this does not take into account the fact we should be contracting
1089            ;;the ideal back to the current open..
1090            (setq answ (loop for ld in lds
1091                             for ggi in ggs
1092                             appending (ldata-simplifications
1093                                         ld
1094                                         :open-g (sftimes int-open-g ggi)
1095                                         :recursive-p t))))
1096
1097 ;        (setq answ (ldata-refinement ldata (car *stop-simplify*)
1098 ;                                     (second *stop-simplify*) :inequality int-open-g))
1099 ;;       (mshow answ)  (format t "**Is the refinement ")
1100 ;        (setq *stop-simplify* nil)
1101 ;        (setq answ (loop for v in answ
1102 ;                                         appending (ldata-simplifications v :open-g open-g :recursive-p t)))
1103          )
1104   ;;priority 1 Linear-dichotomy
1105   ;;         2 gm-prepared equation
1106   ;;         3 any dichotomy
1107   ;;proceed with dich if dich is linear or if found no gm-prepared
1108         ((or lin-dich (null *stop-simplify*))
1109          (cond (dich
1110                 (loop for v in dich
1111                       with so-far = 1
1112                       appending
1113                         (progn
1114                           (setq eqns-modv
1115                                 (loop for facs in all-facs
1116                                       when (not (member v facs :test #'equal))
1117                                         collecting
1118                                           (apply 'gen-ptimes (loop for ter in facs by #'cddr
1119                                                                    when (not
1120                                                                           (may-invertp
1121                                                                             ter so-far))
1122                                                                      collecting ter))))
1123                           (cond ((member nil eqns-modv :test #'eq) (merror "nil should not be here")))
1124                           (cond ((eq v nil) (merror "nil should not be here")))
1125                           (setq ld (make-ldata))
1126                           (setf (ldata-eqns ld) (cons v eqns-modv))
1127                           (setf (ldata-inequality ld)(nplcm gg so-far))
1128                           (setf so-far (nplcm so-far v))
1129                           (cond ((Null *stop-simplify*)
1130                                  (LDATA-SIMPLIFICATIONS ld :open-g open-g
1131                                                         :recursive-p t))
1132                                 (t (list  ld))))
1133                         into list-of-ld
1134                       do (mshow list-of-ld)
1135                       finally
1136                         (cond ((null list-of-ld)
1137                                (setq answ (list (make-ldata :eqns '(1) :inequality 1))))
1138                               (t (setq answ list-of-ld)))))
1139                (t (setq answ (list ldata))))))
1140
1141   (cond ((null answ) (setq answ (list ldata))))
1142 ;;  (cond ((null answ) (setq answ (list (make-ldata :eqns '(1) :inequality 1)))))
1143 ;  (setq answ (delete-redundant-ldata answ)))
1144   answ)
1145 (defun order-variables-by-occurence (list-eqns &aux var-lists all-vars)
1146   (setq var-lists
1147         (loop for v in list-eqns
1148         collecting (list-variables v) ))
1149  (setq all-vars (apply #'zl-UNION var-lists))
1150  (loop for v in all-vars
1151        collecting v
1152        collecting
1153        (loop for lis in var-lists
1154              when (member v lis :test #'eq)
1155              count 1 )))
1156
1157
1158 ;;there is an error in simp-lead: it produced some equations which did not
1159 ;;give ideal containment.
1160 (defun simp-lead (eqns &key (open-g 1) variables &aux fac answ)
1161   (setq eqns (copy-list eqns))
1162   (show (length eqns))
1163   (setq answ (loop for v on eqns
1164         do
1165         (setq fac (non-constant-factors (car v) open-g))
1166         (cond ((and (numberp (car v)) (zerop (car v))) nil)
1167               ((> (length fac ) 2)
1168                (return 'try-dichotomy))
1169               ((null fac)  (return '(1)))
1170               (t (setf (car v) (car fac))
1171                  ))
1172         finally
1173          (return (simp-lead1 eqns :open-g open-g :variables variables))))
1174     (show (or (symbolp answ) (length answ)))
1175   answ)
1176
1177 (defun test-simp-lead (eqns)
1178   (shl eqns)
1179   (let ((answ (simp-lead1 eqns)))
1180     (cond ((atom answ) (show answ))
1181           (t (cond ((not (ideal-subsetp answ eqns))
1182                     (break 't1))
1183                    ((not (ideal-subsetp eqns answ))
1184                     (break 't2)))
1185     (format t "both ok")
1186     (shl answ) answ))))
1187
1188 (defun simp-lead1 (eqns &key (open-g 1) variables &aux  deg-vector cof tem
1189                   answ changed fac varl)
1190   (setq eqns (union-equal eqns))
1191     (cond (variables (setq varl variables ))
1192           (t (setq varl (list-variables eqns))))
1193     (loop named sue for var in varl
1194           for vari on varl
1195           do (setq variables vari)
1196        (setq tem var)
1197           (setq deg-vector (loop for u in eqns
1198                                   collecting (pdegree u var)))
1199           (loop for i in (sort (zl-union (delete 0 (copy-list  deg-vector))) 'alphalessp)
1200                 do
1201                 (loop for j in deg-vector
1202                       for f in eqns
1203                       when (eql j i)
1204                       do
1205                       (or (eq tem var) (error "bad"))
1206                       (setq cof (pcoeff f (list var j 1)))
1207                       (loop for ff in eqns
1208                             for jj in deg-vector
1209                             ;;ff=a*y^i+b f=c*y^j+d
1210                             ;;if j<=i  want to replace ff by
1211                             ;;ff-f*a/c*y^(i-j)  (i.e. c*ff-a*f*y^i-j )
1212                             ;;and this will be ok if
1213                             ;; (denom (a/c) is invertible )
1214                             when (and (>= jj j)(not (eql f ff)))
1215                             when (may-invertp (denom (setq fac
1216                                                            (ratreduce
1217                                                              (pcoeff ff (list var jj 1))
1218                                                              cof)))
1219                                               open-g)
1220                             do (setq eqns (delete ff eqns :test #'equal))
1221 ;                           (format t "replacing for variable ~A .."var)
1222 ;                           (sh ff)
1223 ;                           (format t "using " ) (sh f)
1224                             (setq answ (pdifference
1225                                          (ptimes (denom fac)
1226                                                  ff)
1227                                          (gen-ptimes f (pexpt (list var 1 1)
1228                                                               (- jj j)) (num fac))))
1229 ;                           (format t "   by") (sh answ)
1230
1231                             (cond ((not ($zerop answ))
1232
1233                                    (setq eqns (cons answ eqns))))
1234 ;                           (show (length eqns))
1235                                    (setq changed t)
1236
1237 ;                           (cond ((and changed (equal eqns orig))
1238 ;                                  (merror "the equations did not change")))
1239                             (return-from sue 'start-over)))))
1240     (cond (changed (simp-lead1 eqns :open-g open-g :variables variables))
1241           (t eqns)))
1242
1243
1244 (defun affine-open (list-vars &optional (inequal 1) &aux coords fns)
1245   (cond ((atom (car list-vars))
1246          (cond ((get  (car list-vars) 'disrep)
1247          (setq fns (loop for v in list-vars
1248                collecting (list v 1 1 ))))))
1249         ((mbagp list-vars) (setq fns (st-rat list-vars)))
1250         (t (setq fns list-vars)))
1251    (setq coords (construct-rmap fns))
1252    (make-zopen :coord coords :inv coords :inequality inequal))
1253 ;  (make-zopen :coord (make-rmap  fns fns
1254 ;                               denom 1)
1255 ;             inv      (make-rmap  fns fns
1256 ;                               denom 1)
1257 ;             inequality (st-rat inequal)))
1258 (defun construct-pre-ldata-sheaves (&key s-var data opens)
1259   (check-arg data (or (null data)(null (car data)) (eq (caaar data) 'ldata)) "list of lists of ldata")
1260   (check-arg s-var (or (null s-var) (eq (car s-var) 's-var)) "s-var")
1261   (cond (opens (setq s-var (make-s-var :zopens opens))))
1262   (make-pre-ldata-sheaves :s-var s-var :data data))
1263
1264 ;;;old forms worked
1265 ;(defremember grobner-basis-remember (basis)
1266 ;  (setq *poly-simplifications* (grobner-basis basis)))
1267 ;(defun grobner-remember (basis)
1268 ;    (setq *poly-simplifications* (grobner-basis-remember basis)))
1269
1270
1271 ;;to try to implemement not checking twice and timing out.
1272 (defremember grobner-basis-remember (basis &aux tem answ)
1273   (setq   answ (catch 'took-too-long
1274                 (setq tem  (grobner-basis basis))))
1275   (cond (tem tem)
1276         (t (list 'took-too-long *timed-grobner-basis* *poly-simplifications*))))
1277
1278 (defun grobner-remember (basis &aux simps)
1279   (setq  simps (grobner-basis-remember basis))
1280   (cond ((eq (car simps ) 'took-too-long)
1281          (cond ((eq (second simps) *timed-grobner-basis*)
1282                 (setq *poly-simplifications* (third simps))(throw 'took-too-long 'took-too-long))
1283                (t (setf (gethash  (get 'grobner-basis-remember :memory-table) basis) nil)
1284                   (grobner-remember basis))))
1285         (t (setq *poly-simplifications* (grobner-basis-remember basis)))))
1286
1287 ;(defun delete-redundant-ldata (list-ld)
1288 ;  (loop for v in list-ld
1289 ;       do
1290 ;       (loop for w in (ldata-eqns v)
1291 ;             when (numberp w)
1292 ;             do (setq list-ld (delete v list-ld))))
1293 ;  (cond ((>= (length list-ld) 2)
1294 ;        (loop for v on list-ld
1295 ;              do (grobner-basis (ldata-eqns (first  v)))
1296 ;              (loop for u in list-ld
1297 ;                    when(and (not (eql u (first v)))
1298 ;                             (ideal-subsetp (ldata-eqns u)
1299 ;                                        (ldata-eqns (first v))
1300 ;                                        :reset-basis nil :verbose nil))
1301 ;                    do (format t "~%deleting redundant component..")
1302 ;                    (setq list-ld (delete (first v) list-ld))
1303 ;                    (return t)))))
1304 ;
1305 ;  list-ld)
1306
1307 (defun variety-ldata-subset (ld1 ld2 &key (open-g 1) ignore-ldata-inequalities)
1308   (cond (ignore-ldata-inequalities
1309         (multiple-value-bind (cont unit)
1310           (grobner-subset (ldata-eqns ld2) (ldata-eqns ld1) open-g)
1311           (values cont (and unit 'empty))))
1312         (t
1313   (let ((gg (nplcm open-g (ldata-inequality ld1))))
1314     (and (unit-idealp (cons (ldata-inequality ld2) (ldata-eqns ld1) )
1315                       gg)
1316          ;;return only one value since we don't want to test if ld1 is empty on open-g
1317          (values(grobner-subset (ldata-eqns ld2) (ldata-eqns ld1)
1318                          (nplcm gg (ldata-inequality ld2)))))))))
1319
1320
1321 ;;;version 4.0 of delete has check of
1322 ;; V(I1,c1)^PV(I2,c2) where I1,c1 has inequality
1323 ;;if I2^PI1[c1^-1,c2^-1] and <c2,I1>[c1^-1] = <1>
1324 ;;ie if where c1 and c2 nonzero get containment, and also insist that V(I1,c1) does
1325 ;;not meet the c2=0 set.
1326
1327 ;(defun delete-redundant-ldata (list-ld &key     (gg 1)
1328 ;                              &aux cc cint use-inverse (complexity-for-inverse 50))
1329 ;  (loop for v in list-ld
1330 ;       do
1331 ;       (loop for w in (ldata-eqns v)
1332 ;             when (numberp w)
1333 ;             do (setq list-ld (delete v list-ld))))
1334 ;  (cond ((>= (length list-ld) 2)
1335 ;        (loop for v on list-ld
1336 ;              do
1337 ;              (setq cc (nplcm gg  (ldata-inequality (first v))))
1338 ;              (loop for u in list-ld
1339 ;                       when (not (eql u  (first v)))
1340 ;                       do
1341 ;                       (setq cint (nplcm cc (ldata-inequality u)))
1342 ;
1343 ;                       (multiple-value-bind (cont unit)
1344 ;                           (grobner-subset (ldata-eqns u)
1345 ;                                           (ldata-eqns (first  v))
1346 ;                                           cint)
1347 ;                         (cond ((and (or unit
1348 ;                                    cont) (unit-idealp (cons (ldata-inequality u)
1349 ;                                                             (ldata-eqns (first v)))
1350 ;                                                       cc))
1351 ;                                (setq list-ld (delete (first v) list-ld))
1352 ;                                (format t "~%Deleting redundant component")
1353 ;                                (des (first v))
1354 ;                                (format t "~%Because of ")
1355 ;                                (show (list unit cont))
1356 ;                                (cond ((null unit)
1357 ;                                       (des u)))
1358 ;                                (return t))))))))
1359 ;  list-ld)
1360 ;
1361 ;;;putting use-inverse = t can be dangerous I am moving it to the aux
1362 ;;;this version is not correct but involves less calculation.
1363 ;;;maybe its ok if the resultant things are prime.
1364 ;
1365 ;(defun delete-redundant-ldata (list-ld &key     (gg 1)
1366 ;                              &aux use-inverse (complexity-for-inverse 50))
1367 ;  (loop for v in list-ld
1368 ;       do
1369 ;       (loop for w in (ldata-eqns v)
1370 ;             when (numberp w)
1371 ;             do (setq list-ld (delete v list-ld))))
1372 ;  (cond ((>= (length list-ld) 2)
1373 ;        (loop for v on list-ld
1374 ;              do (loop for u in list-ld
1375 ;                       when (not (eql u  (first v)))
1376 ;                       do
1377 ;                       (cond (use-inverse (setq gg
1378 ;                                                (nplcm gg (ldata-inequality (first v))))))
1379 ;;                      (cond ( (> (pcomplexity gg) complexity-for-inverse)
1380 ;;
1381 ;;                             (format t "%Not using inverse for deletion :.."
1382 ;;                                     (des gg))
1383 ;;                                                             (setq gg 1)))
1384 ;                       (multiple-value-bind (cont unit)
1385 ;                           (grobner-subset (ldata-eqns u)
1386 ;                                           (ldata-eqns (first  v))
1387 ;                                           gg)
1388 ;                         (cond ((or unit
1389 ;                                    cont)
1390 ;                                (setq list-ld (delete (first v) list-ld))
1391 ;                                (format t "~%Deleting redundant component")
1392 ;                                (des (first v))
1393 ;                                (format t "~%Because of ")
1394 ;                                (show (list unit cont))
1395 ;                                (cond ((null unit)
1396 ;                                       (des u)))
1397 ;                                (return t))))))))
1398 ;  list-ld)
1399
1400
1401 ;;timed
1402 (defun delete-redundant-ldata (list-ld &key (open-g 1)
1403                                (gg 1)   ignore-ldata-inequalities
1404                                &aux redundant )
1405   (setq gg (sftimes open-g gg))
1406   (setq list-ld (loop for v in list-ld
1407         when (not (unit-idealp (ldata-eqns v) (nplcm (ldata-inequality v) gg)))
1408         collecting v
1409         else do (format t "~2%Deleting the empty ldata ")
1410         (fsh v)))
1411   (loop for v in list-ld
1412         for vi from 0
1413         do
1414         (loop for w in list-ld
1415               for wi from 0
1416               when (not (or (member vi redundant :test #'equal)
1417                             (member wi redundant :test #'equal)
1418                             (eql vi wi)))
1419               when
1420               (variety-ldata-subset v w :open-g gg :ignore-ldata-inequalities
1421                                     ignore-ldata-inequalities)
1422               do (push vi redundant)
1423               (format t "~2%Deleting redundant component")
1424               (fsh v)
1425               (format t "~%contained in ") (fsh w)))
1426   (loop for v in list-ld
1427         for vi from 0
1428         when (not (member vi redundant :test #'equal))
1429         collecting v))
1430
1431 (defun describe-ldata (ld)
1432   (cond ((eq (first ld) 'ldata)
1433          (format t "~%The ~A  equations are .." (length (ldata-eqns ld)))
1434          (shl (ldata-eqns ld))
1435          (format t "~%The inequality for the component calculation was ..")
1436          (sh (ldata-inequality ld)))
1437         (t (mapcar 'describe-ldata ld))))
1438 ;
1439 ; (setq mons
1440 ;      (loop for v in (cdr mat)
1441 ;           collecting
1442 ;           (loop for u in '($x1 $x2 $x3 $x4 $x5)
1443 ;                 for w in (cdr v)
1444 ;                 with answ = 1
1445 ;                 do (setq tem (st-rat u))
1446 ;                 (setq answ (ptimes answ (pexpt tem w)))
1447 ;                 finally (return answ))))
1448
1449 ;;will return a list of ldata
1450
1451 (defun $add_to_solution_tree (a-list &aux tem tem1)
1452   (check-arg a-list $listp "macsymya list")
1453   ;;check this case had not been done yet:
1454   (loop named pat for u in (cdr a-list)
1455         do
1456         (check-arg u $listp "macsyma list")
1457         (setq tem (from-case  u))
1458         (loop  for v in $solution_tree
1459               when (equal (from-case a-list) tem)
1460               do (cond ((yes-or-no-p "Remove from $solution_tree case ~A" tem)
1461                         (loop for v in $solution_tree
1462                               when (not (and (initial-equal
1463                                                (setq tem1(from-case v)) tem)
1464                                              (>= (length tem1)  tem)))
1465                               collecting v into tempp
1466                               finally (setq $solution_tree tempp)
1467                               (return-from pat  'fixed))))))
1468   (setq $solution_tree
1469         (cons '(mlist) (append (cdr a-list) (cdr $solution_tree))))
1470   ($sort_solution_tree))
1471
1472
1473 (defun $simplify_relative (system poly)
1474   ($grobner_basis (list '(mlist) poly))
1475    (append  (delete 0 ($totaldisrep ($polysimp system))) (list poly)))
1476
1477
1478 (defun $list_irreducible_factors(poly &aux answer)
1479   (cond
1480     ((atom poly ) (list '(mlist) poly))
1481     ((MBAGP poly) (CONS (CAR poly) (MAPCAR #'$list_irreducible_factors (CDR poly))))
1482     ((eq (caar poly) 'mtimes)
1483      (cond ((member 'factored (car poly) :test #'eq)
1484             (loop for u in (cdr poly)
1485                   do
1486                   (cond ((atom u)
1487                          (push u answer))
1488                         ((eq (caar u) 'mplus)
1489                          (push u answer))
1490                         ((eq (caar u) 'mexpt)
1491                          (push (second u) answer )))
1492                   finally (return ($sort (cons '(mlist) (zl-UNION  answer))))))
1493            (t ($list_irreducible_factors ($factor poly)))))
1494     ((eq (caar poly) 'mexpt) ($list_irreducible_factors ($factor (second poly))))
1495     ((eq (caar poly) 'mplus)
1496      (cond ((member 'irreducible (car poly) :test #'eq) (list '(mlist) poly))
1497                                    (t ($list_irreducible_factors ($factor poly)))))
1498     (t (merror "how did you get here"))))
1499
1500
1501 (defun $show_solution_tree (&optional full-format &aux tem)
1502   (loop for v in (cdr $solution_tree)
1503         do
1504         (setq tem (from-case v))
1505         (case (length tem)
1506           (0 (format t "~%Cases for the solution tree")
1507           (cond (full-format  (i-$grind  v))))
1508           (1 (format t "~%  Case: ~A"(string-grind (car tem)))
1509                   (cond (full-format  (i-$grind v))))
1510           (2 (format t "~%    Subcase ~A"(string-grind (second tem))))
1511           (3 (format t "~%       Subsubcase ~A") (string-grind (third tem)))
1512           (t (format t "~%         Sub......~A" ) (string-grind (nthcdr 3 tem)))))
1513   '$done)
1514
1515 (defun $eliminate_nonzero_factors (system &rest l)
1516   "rather crude see cancel below"
1517   (cond ($nonzero_factors
1518          (setq l (append l (cdr $nonzero_factors)))))
1519   (cond (L
1520   (setq system (div* system (power (apply 'mul*  l) 10)))
1521      ($numerator system))
1522         (t system)))
1523
1524
1525 (defun $cancel_factors_and_denominators (variety factors-to-elim
1526                                          &optional (homogeneous nil)
1527                                          &aux system rat-variety type)
1528   "tries to return same type"
1529   (cond ((atom factors-to-elim) (setq factors-to-elim (list factors-to-elim)))
1530         (($listp factors-to-elim) (setq factors-to-elim (cdr factors-to-elim)))
1531         (t factors-to-elim))
1532   (check-arg variety  $listp  "macsyma list")
1533   (setq type (cond ((affine-polynomialp (second variety)) 'polynomial)
1534                    ((rational-functionp (second variety)) 'rational-function)
1535                    (($ratp (second variety)) '$rat)
1536                    (t 'general)))
1537   (setq rat-variety
1538         (loop for v in (cdr variety)
1539                collecting (new-rat v)))
1540   (setq system rat-variety)
1541   (cond (homogeneous
1542          (loop for w in rat-variety
1543                 with tem = 1
1544                 do (setq tem (nplcm tem (denom w)))
1545                 finally
1546                 (setq tem (cons tem 1))
1547                 (setq system
1548                       (loop for w in system
1549                              collecting
1550                              (num (rattimes tem
1551                                             w t))))))
1552         (t (setq system
1553                  (loop for w in rat-variety
1554                         collecting (num w)))))
1555   (show factors-to-elim)
1556   (loop for v in factors-to-elim
1557          with tem
1558          do
1559          (setq tem   (pexpt (st-rat v) 7))
1560          (cond (homogeneous
1561                 (loop for w in system
1562                        do (setq tem (npgcd tem w))
1563                        (cond ((numberp tem) (return 'done)))
1564                        finally
1565                        (setq system
1566                              (loop for w in system
1567                                     collecting
1568                                     (pquotient w tem)))))
1569                (t (setq system
1570                         (loop for w in system
1571                                collecting
1572                                (num (ratreduce w tem)))))))
1573   (case type
1574     (general
1575      (loop for w in system
1576             collecting
1577             (new-disrep w)
1578             into tem
1579             finally (return (cons '(mlist) tem))))
1580     ($rat (cons '(mlist) (mapcar 'header-poly system)))
1581     (t system)))
1582
1583 (defun $blowup_chart (chart variety &optional homogeneous
1584                       &aux divisor transform)
1585   "Performs the substitutions of chart and cancels the exceptional divisor
1586  and gets rid of denominators.  It does this homogeneously if the value
1587  of homogeneous is not nil."
1588
1589   (check-arg chart $listp "macsyma list")
1590   (check-arg variety $listp "macsyma list")
1591   (loop for w in (cdr chart)
1592         with tem = (third (second chart))
1593         do
1594         (show tem )
1595         (setq tem ($gcd tem (third w)))
1596         finally (setq divisor tem)
1597         (format t "~%The exceptional divisor is ~A" tem))
1598   (setq transform  ($sublis chart variety))
1599   ($cancel_factors_and_denominators transform (list '(mlist)
1600                                                 divisor)
1601                                     homogeneous))
1602
1603 ;;;want this to blowup where the locus is given by a complete intersection.
1604 ;(defun $blowup (variety locus-to-blowup new-var-prefix)
1605 ;  (check-arg $listp locus-to-blowup "macsyma list")
1606 ;  (loop for u in (cdr locus-to-blowup)
1607 ;       collecting))
1608
1609 (defun my-minor (mat  rows-to-take cols-to-take)
1610   ($determinant ($sub_matrix mat rows-to-take cols-to-take)))
1611
1612 (defun $wedge_matrix (mat n)
1613  (let ((tabl (list-tableaux n ($length (second mat)) )))
1614    (show tabl)
1615    (loop for u in tabl
1616         collecting
1617         (loop for v in tabl
1618               collecting (my-minor mat u v) into tem
1619               finally (return (cons '(mlist) tem)))
1620         into temm
1621         finally (return (cons '($matrix) temm)))))
1622
1623 (defun $sublis_and_add (eqns expr &rest  elimin &aux answer)
1624  (setq answer ($sublis eqns expr))
1625  (show elimin)
1626  (setq answer (append answer (mapcar 'bring-to-left-side (cdr  eqns ))))
1627  (setq answer ($factor (apply '$eliminate_nonzero_factors
1628                                        answer elimin))))
1629
1630 ;(defun pcomplexity  (poly)
1631 ;  (cond ((atom poly) 0)
1632 ;       (t (loop for (deg cof) on (cdr poly) by #'cddr
1633 ;                summing (+ deg (pcomplexity cof))))))
1634
1635 ;(defun $remove_linears (eqns &aux (best t)used-eqns pcompw rat-eqns tem)
1636 ;  (setq rat-eqns
1637 ;  (loop for v in (cdr eqns) until (eq v '$case)
1638 ;       collecting (st-rat v) ))
1639 ;  (loop while best
1640 ;       do
1641 ;       (loop for w in rat-eqns
1642 ;             initially (setq best   nil)
1643 ;             with prev-compl = 4000000
1644 ;             when (setq tem (coll-linear w))
1645 ;             do
1646 ;             (show tem)
1647 ;             (cond ((< (setq pcompw (pcomplexity w)) prev-compl)
1648 ;                    (setq prev-compl pcompw ) (setq best (cons (car tem) w))))
1649 ;             finally (cond (best
1650 ;                            (format t "~%Substituting for ~A by using"
1651 ;                                    (get (car best ) 'disrep))
1652 ;                            (sh (cdr best))
1653 ;                            (loop for v in rat-eqns
1654 ;                                  collecting
1655 ;                                  (psublis (list(cons (car best)
1656 ;                                                      (pcoeff (cdr best) 1
1657 ;                                                              (list (car best)))))
1658 ;                                           (pcoeff (cdr best) (car best))
1659 ;                                           (cdr best))
1660 ;                                  into new-eqns
1661 ;                                  finally (setq rat-eqns new-eqns)))))
1662 ;       (push (cdr best) used-eqns) )
1663 ;  (cond (used-eqns
1664 ;        (cons '(mlist) (append (loop for v in rat-eqns
1665 ;                      collecting (header-poly v))
1666 ;                (member '$case eqns :test #'eq)
1667 ;                (loop for v in used-eqns
1668 ;                      collecting (header-poly v)))))))
1669
1670
1671 (defun $remove_Linears (eqns  &aux subs lins)
1672   (loop for v in (cdr eqns)
1673         when (setq lins (cdr  ($coll_linear eqns)))
1674         appending lins into lin-vars
1675         and
1676         collecting v into lin-eqns
1677         finally (cond (lins
1678                        (setq subs
1679                              ($fast_linsolve (cons '(mlist) lin-eqns)
1680                                              (cons '(mlist)
1681                                                    (zl-UNION lin-vars))))
1682                        (format t "~%Eliminating linear variables ~A" lin-vars)
1683                        ))
1684         (cond (subs  (setq eqns ($sublis subs eqns))
1685                      (setq eqns (append (delete 0 ($ratsimp eqns))
1686                                         (loop for v in subs collecting
1687                                               ($numerator
1688                                                 (sub* (second v) (third v)))))))
1689               (t eqns))))
1690
1691
1692 ;
1693 ;(defmacro te (x v)
1694 ;  (cond ((functionp x) `(,x ,v))
1695 ;       (t 5)))
1696
1697 (defmacro find-minimal (in-list ordering &optional   such-that ind   )
1698   (cond (such-that
1699          (cond ((functionp such-that)(setq ind '-ind-)
1700                 (setq such-that `(,such-that -ind-)))
1701                (t
1702          (check-arg ind (not (null ind)) "non nil.  Must specify index")))
1703          ` (loop for ,ind in ,in-list
1704                  with prev-min
1705                  when  ,such-that
1706                  do (cond (prev-min
1707                            (cond ((funcall ,ordering ,ind prev-min)
1708                                   (setq prev-min ,ind))))
1709                           (t (setq prev-min ,ind)))
1710                  finally (return prev-min)))
1711         (t   `(loop for v in ,in-list
1712                     with prev-min
1713                     do (cond (prev-min
1714                               (cond ((funcall , ordering v prev-min)
1715                                      (setq prev-min v))))
1716                              (t (setq prev-min v)))
1717                     finally (return prev-min)))))
1718
1719
1720 ;(defun find-minimal (in-list &optional  such-that ordering  &aux  prev-min)
1721 ;    (cond (such-that
1722 ;        (loop for v in in-list
1723 ;              when (funcall such-that v)
1724 ;              do (cond (prev-min
1725 ;                        (cond ((funcall ordering v prev-min)
1726 ;                               (setq prev-min v))))
1727 ;                       (t (setq prev-min v)))))
1728 ;       (t  (loop for v in in-list
1729 ;                         do (cond (prev-min
1730 ;                           (cond ((funcall ordering v prev-min)
1731 ;                                  (setq prev-min v))))
1732 ;                          (t (setq prev-min v))))))
1733 ;
1734 ;  prev-min)
1735
1736 ;(defun find-good-variable-order (ideal)
1737 ;  (setq vars (mapcar '$list_variables (cdr ideal)))
1738 ; (setq vars  (mapcar 'cdr vars))
1739 ; (setq all-vars (union (apply 'append vars)))
1740 ; (setq f #'(lambda (v) (loop for w in vars
1741 ;                            when (member v :test #'eq) vars
1742 ;                            count1
1743 ; (sort all-vars #'(lambda (u v)
1744 ;
1745 ; (setq vars (sort vars #'(lambda( u v)(< (length u) (length v)))))
1746 ; (setq all-vars (union (apply 'append  vars))))
1747
1748 (defun find-good-variable-order (ideal &aux f all-vars vars)
1749   (declare (special  f vars))
1750   (setq vars (mapcar '$list_variables (cdr ideal)))
1751   (setq vars  (mapcar 'cdr vars))
1752    (setq all-vars (zl-union (apply 'append  vars)))
1753   (setq f #'(lambda (v) (loop for w in vars when (member v vars) count 1 into tem
1754                               finally (return tem))))
1755    (sort all-vars #'(lambda( u v)(< (funcall f u) (funcall f v)))))
1756
1757 (defun may-invertp (poly invertible-poly)
1758   (if ($zerop poly)
1759       nil
1760       (numberp (denom (ratreduce invertible-poly (square-free poly))))))
1761
1762 (defvar $char_set nil)
1763 (defvar $ideal nil)
1764
1765 (defun poly-linearp (poly  var may-invert &aux cof)
1766   (cond ((numberp poly ) nil)
1767         ((eql (pdegree poly var) 1)
1768          (setq cof (pcoeff poly (list var 1 1)))
1769          (cond ((atom cof) cof)
1770                ((may-invertp cof may-invert) cof)
1771                (t nil)))))
1772
1773 ;(defun $te (eqns may-invert &aux ans)
1774 ;  (setq ans (remove-linears (mapcar 'st-rat (cdr eqns) ) (st-rat may-invert)))
1775 ;  (cons '(mlist) (mapcar 'new-disrep ans)))
1776
1777 ;;;this seems the most efficient.  If may-invert is x^2-1 then u*(x+1)+h(w,v)
1778 ;;;will cause the replacement u--> -h(w,v)/(x+1) in the remaining eqns.
1779 ;;;then it calls itself on the remaining equations.
1780 (defun remove-linears (eqns &optional (may-invert 1) &aux sub cof)
1781   (check-arg eqns (affine-polynomialp (car eqns)) "first term not  polynomial")
1782   (loop for v in eqns do
1783        (block sue
1784          (loop for vari in (list-variables v)
1785             when (setq cof (poly-linearp v  vari may-invert))
1786             do (setq sub (cons vari (pdifference (ptimes (list vari 1 1) cof) v))
1787                      denom cof)
1788               (loop for uu in eqns
1789                  when (not (equal uu v))
1790                  collecting (psublis (list sub) denom uu) into new-eqn
1791                  finally (return-from sue (cons (num (ratreduce v cof))
1792                                                 (remove-linears new-eqn may-invert))))))
1793      finally (return eqns)))
1794
1795 (defvar *to-invert* nil)
1796
1797 (defun $ritt_set (ideal-generators &optional (may-invert 1) &aux rat-ideal tem id)
1798   "makes $char_set from an ideal-generators pushing the leading coefs into *to-invert*"
1799   (setq $ideal (setq id (make-ideal)))
1800
1801   (setf (ideal-variable-correspondence id)
1802         (order-variables (find-good-variable-order ideal-generators)))
1803      (ideal-variable-correspondence id)
1804     (setq rat-ideal
1805           (loop for u in (cdr ideal-generators)
1806                 with tem
1807                 when (not (pzerop (setq tem (st-rat u))))
1808                 collecting tem ))
1809     (setq may-invert (st-rat may-invert))
1810     (setq rat-ideal (remove-linears rat-ideal may-invert))
1811     (setq vars (reverse (sort  (list-variables rat-ideal)'pointergp)))
1812     (show vars)
1813     (show (car rat-ideal))
1814     (setq *to-invert* (list may-invert))
1815     (setq rat-ideal
1816           (loop for v in rat-ideal
1817                 collecting (remove-common-factors v may-invert)))
1818     (setq $char_set
1819           (loop for v in vars
1820                 with temm
1821                 when (setq temm (find-minimal
1822                                  rat-ideal
1823                             #'(lambda (u v)(< (p-le u) (p-le v)))
1824                                  (and (listp pol) (eq (p-var pol) v)) pol))
1825
1826                 do (show temm); (sh v)
1827                 and
1828                 collecting  temm))
1829     (loop for v in rat-ideal
1830           do
1831           (setq tem (ritt-reduce v $char_set))
1832           (mshow tem)
1833           (show (length $char_set))
1834           (cond ((pzerop tem) nil)
1835                 (t (setq $char_set (add-to-chain tem $char_set)))))
1836     (setq $char_set (eliminate-factors $char_set *to-invert*))
1837     (setf (ideal-char-set id) $char_set)
1838     (setf (ideal-localization  id) *to-invert*)
1839     ($disrep_ideal id))
1840
1841 (defun eliminate-factors (good factors-to-elim )
1842   (loop for v in factors-to-elim
1843         with mon
1844         do (setq mon (pexpt v 5))
1845         (setq good
1846         (loop for w in good
1847               collecting (num (ratreduce w mon)))))
1848   good)
1849
1850 (defun header-fake (expr)
1851   (cond ((or (numberp expr)(get  (car expr) 'disrep))
1852          (cons '(mrat nil nil) (cons expr 1)))
1853         ((and (or (numberp (car expr))( get (caar expr) 'disrep))
1854               (or (numberp (cdr expr)) (get (cadr expr) 'disrep)))
1855          (cons '(mrat nil nil) expr))
1856         (t (merror "not a rat'l fun or poly (at least no disrep prop)"))))
1857 (defun fake-disrep (pol)
1858    ($totaldisrep (header-fake pol)))
1859
1860 (defun $disrep_ideal (ideal)
1861   (list '(mlist)
1862         (cons '(mlist)(mapcar 'fake-disrep (ideal-char-set ideal)))
1863         (cons '(mlist)(mapcar 'fake-disrep (ideal-localization ideal)))))
1864
1865 ;(defun ideal-disrep item ideal)
1866 ;  (loop with added = t
1867 ;       while added
1868 ;       do (setq added nil)
1869 ;       (loop for v in rat-ideal
1870 ;             do (show (length $char_set))
1871 ;              (setq tem (ritt-reduce v $char_set))
1872 ;             when (not (pzerop tem))
1873 ;             do
1874 ;             (mshow tem)
1875 ;             (setq $char_set (add-to-chain tem $char_set))
1876 ;             (show (length $char_set))
1877 ;             else do (setq rat-ideal (delete v rat-ideal)) (format t "its zero")
1878 ;             finally (setq rat-ideal nil)))
1879 ;  $char_set)
1880
1881 (defun square-free (p)
1882   (let ((facts (psqfr p)))
1883     (if (cddr facts)
1884         (loop for v in (cddr facts) by #'cddr
1885            with answer = (car facts) do
1886              (setq answer (ptimes v answer))
1887            finally (return answer))
1888         (car facts))))
1889
1890 (defun $square_free_numerators (expr)
1891   (cond ((atom expr) expr)
1892         ((affine-polynomialp expr)(square-free expr))
1893         ((rational-functionp expr) (square-free (num expr)))
1894         (($ratp expr) (cons (car expr) (cons (square-free (num (cdr expr)))
1895                                              1)))
1896         ((mbagp expr) (cons (car expr)
1897                             (mapcar '$square_free_numerators (cdr expr))))
1898         (t (let ((tem (new-rat expr)))
1899
1900              (header-poly (cons (square-free (num tem)) 1))))))
1901
1902
1903
1904 (defun add-to-chain (poly chain)
1905   (setq poly (square-free poly))
1906   (cond ((eql poly 0) chain)
1907         (t  (loop for v in chain
1908                        when (eq (p-var poly)(p-var v))
1909                        do (cond ((< (p-le poly) (p-le v))(show 'deleting)
1910                                  (setq chain (delete v chain :test #'equal)))
1911                                 ((>= (p-le poly)(p-le v))
1912                                  (mshow poly v) (merror "not reduced")))
1913                        finally(format t "~%adding .." )
1914                        (sh poly)
1915                        (setq chain
1916                              (sort (cons poly chain)
1917                                          #'(lambda (u v)
1918                                              (pointergp (p-var u) (p-var v)))))
1919                        (return chain)))))
1920
1921
1922 (defun order-variables (list-of-vars &aux vc)
1923   (cond (($listp list-of-vars) (setq list-of-vars (cdr list-of-vars))))
1924   (setq vc (make-variable-correspondence))
1925   (loop for v in list-of-vars
1926         for i from 1
1927         with w
1928         collecting (setq w (gensym)) into gen
1929         do
1930         (setf (get w 'disrep) v)
1931         (setf (symbol-value w) i)
1932         finally (setf (vc-genvar vc) gen)
1933         (setf (vc-varlist vc) (copy-list list-of-vars)))
1934   vc)
1935
1936 (defmacro with-vc (var-corr &rest body)
1937  `(let ((*genvar* (vc-genvar ,var-corr))
1938         (*varlist* (vc-varlist ,var-corr)))
1939     (unwind-protect
1940       (progn ,@ body)
1941       (setf (vc-genvar ,var-corr) *genvar*)
1942       (setf (vc-varlist ,var-corr ) *varlist*))))
1943
1944 ;(defun sh (f)
1945 ;  (cond ($display2d (displa (header-poly f)))
1946 ;       (t (string-grind (header-poly f) :stream t))))
1947 ;(defun shl (l) (mapcar 'sh l))
1948 ;(defun shl (l)
1949 ;  (cond ($display2d (mapcar 'sh l))
1950 ;       (t (loop for v in l
1951 ;               for i from 0
1952 ;               initially (format t "~%[")
1953 ;               do  (sh (header-poly v))
1954 ;               when (< i (1- (length l))) do(format t ",~%")
1955 ;               finally (format t "]")))))
1956
1957
1958 ;(defun add-to-chain (poly chain)
1959 ;  (cond ((eql poly 0) chain)
1960 ;       (t
1961 ;  (loop for v on chain
1962 ;       collecting (car v) into tem
1963 ;       when (eq (p-var poly)(p-var (car v)))
1964 ;       do (cond ((< (p-le poly) (p-le (car v)))
1965 ;                 (return (nconc tem (cdr v)))))
1966 ;
1967 ;       finally (return (cons poly chain))))))
1968
1969
1970
1971
1972 (defun must-ritt-reducep (poly ch-set)
1973   (cond ((atom poly) nil)
1974         (t
1975   (loop for v in ch-set
1976         when (eq (p-var poly) (p-var v))
1977         do (cond ((>= (p-le poly) (p-le v))
1978                   (return t)))))))
1979
1980
1981 (defun ritt-reduce (poly ch-set &aux tem)
1982   ;;assumes the ch-set is sorted by main variables
1983   (loop while (must-ritt-reducep poly ch-set)
1984         do
1985         (loop for v in ch-set
1986               when (numberp poly)
1987               do (return poly)
1988               when (eq (p-var poly) (p-var v))
1989               do (cond ((>= (p-le poly) (p-le v))
1990                         (setq poly (second (setq tem (vdivide poly v))))
1991 ;                       (push-new (third tem) *to-invert*) 'equal
1992                          (cond ((not (numberp (third tem)))
1993                                 (setq *to-invert*
1994                                       (list (nplcm (third tem)
1995                                                     (car *to-invert*))))))
1996                         (show *to-invert*)))
1997               finally (return poly)))
1998   poly)
1999
2000 ;;rational-map  will be type rmap
2001 ;;  funs (f1 f2 ... fn) denom gg
2002 ;;zopen set : (((p1,  pn)gp) (q1,..qn) gq),gg
2003 ;;where zi=qi(x1,x2,  xn)/gp and xi=pi(z1,..,zn)/gq
2004
2005 (defmfun my-testdivide (x y)
2006   (ignore-rat-err (pquotient x y)))
2007
2008 (defun new-testdivide (f g &aux quot)
2009   (setq quot (ratreduce f g))
2010   (setq quot (cond ((equal (denom quot) 1) (num quot))
2011         (t nil)))
2012   (iassert (equal quot (my-testdivide f g)))
2013   quot)
2014
2015 (defun eliminate-highest-power-dividing-homogeneously  (list-fns divisor &optional
2016                                                         (highest-deg 10000)
2017                                                         &aux quot )
2018   (cond ((and (numberp divisor)(Equal (abs divisor) 1))
2019          list-fns)
2020         (t
2021          (loop named sue
2022                with prev-list-quotients = list-fns
2023                for i from 1
2024                do
2025                (loop for v in prev-list-quotients
2026                      do
2027                      (setq quot        (my-testdivide v divisor))
2028                      when (or  (null quot) (> i highest-deg))
2029                      do
2030                      (return-from sue prev-list-quotients)
2031                      else
2032                      collecting quot into new-quotients
2033                      finally (setq prev-list-quotients new-quotients))))))
2034
2035 ;
2036 ;(defun eliminate-common-factors (list-fns &aux facts simple-fn)
2037 ;  (loop for v in list-fns collecting (pcomplexity v)
2038 ;             into tem
2039 ;             when (not (pzerop v))
2040 ;             minimize (car tem) into the-min
2041 ;             finally
2042 ;              (setq simple-fn (nth (setq the-min (find-position-in-list the-min tem)) list-fns)))
2043 ;  (setq facts (npfactor simple-fn))
2044 ;  (loop for (pol deg)  on facts by #'cddr
2045 ;       with quot = list-fns
2046 ;       do (setq quot (eliminate-highest-power-dividing-homogeneously quot pol deg))
2047 ;       finally (return quot)))
2048
2049 (defun sort-remember (list predicate &key key)
2050   (setq list (loop for i from 0 for v in list
2051         collecting (cons i v)))
2052   (setq list (sort list predicate :key (or (and key  #'(lambda (x)  (key (cdr x))))
2053                                 #'cdr)))
2054   (loop for v on list
2055         collecting (caar v) into ordering
2056         do (setf (car v) (cdar v))
2057         finally (return (values list ordering))))
2058
2059 (defun un-sort (list ordering &aux (newlist (make-list (length ordering))))
2060   (loop for v in ordering
2061         for w in list
2062         do (setf (nth v newlist) w))
2063   newlist)
2064
2065 (defun sort-by-ordering (list ordering)
2066   (loop for v in ordering
2067         collecting (nth v list)))
2068
2069 (defun copy-sort (list predicate &key key slow-key remember)
2070   "copies the list and sorts by predicate. If slow-key, it uses only one application of the key per item.
2071 If slow-key or remember it returns a second value : the ordering so that (un-sort  result ordering)
2072 would restore the list"
2073   (cond  (slow-key
2074           (multiple-value-bind (result order)
2075               (sort-remember (mapcar slow-key list) predicate)
2076             (values (sort-by-ordering list order) order)))
2077          (remember (sort-remember list predicate :key key))
2078          (t (sort (copy-list list) predicate :key key))))
2079
2080 ;;this works for rational functions as well should be fixed to use complexity to.
2081 (defun eliminate-common-factors (list-functions &aux num-gcd denom-gcd firs)
2082   "try to speed up by putting simple ones first"
2083   (multiple-value-bind (list-fns order)
2084       (copy-sort list-functions '< :slow-key
2085                  #'(lambda (x) (+ (pcomplexity (function-numerator x))
2086                                   (pcomplexity (function-denominator x)))))
2087     (setq firs (first list-fns))
2088     (cond ((affine-polynomialp firs)
2089            (setq num-gcd   firs)
2090            (setq denom-gcd 1))
2091           ((rational-functionp firs)
2092            (setq num-gcd (num  firs))
2093            (setq denom-gcd (denom firs))))
2094     (loop for v in (cdr list-fns)
2095           when (affine-polynomialp v)
2096           do (setq num-gcd (pgcd num-gcd v))
2097           (setq denom-gcd 1)
2098           else
2099           do (check-arg v 'rational-functionp "rational function")
2100           (setq num-gcd (pgcd num-gcd (num v))))
2101     (un-sort (loop for v in list-fns
2102                    collecting
2103                    (st-rat (cons (pquotient (function-numerator v) num-gcd)
2104                                  (pquotient (function-denominator v) denom-gcd))))
2105              order)))
2106
2107
2108
2109
2110 (defun reduce-rational-map-old (rmap &aux answ (gc (rmap-denom rmap)) (genvar *genvar*))
2111   (let ((genvar (nreverse(sort (list-variables (cons (rmap-denom rmap)(rmap-fns rmap)))
2112                       'pointergp))))
2113   (loop for v in (rmap-fns rmap)
2114         do
2115         (setq gc (pgcd gc v)))
2116   (setq answ (make-rmap))
2117   (setf (rmap-fns answ)
2118         (loop for v in (rmap-fns rmap)
2119               collecting (pquotient v gc)))
2120   (setf (rmap-denom answ) (pquotient (rmap-denom rmap)gc))
2121   answ))
2122
2123
2124
2125 (defun reduce-rational-map (rmap &aux red-fns fns (genvar *genvar*) answ)
2126   (let ((genvar (nreverse(sort (list-variables (cons (rmap-denom rmap)(rmap-fns rmap)))
2127                                'pointergp))))
2128     (setq fns (cons (rmap-denom rmap) (rmap-fns rmap)))
2129     (setq red-fns   (eliminate-common-factors fns))
2130     (setq answ   (zl-copy-structure rmap rmap- fns (cdr red-fns)
2131                                          denom (car red-fns)))
2132     answ))
2133
2134 (defun convert-rmap-to-new (name)
2135          (zl-copy-structure
2136                name rmap-
2137                     fns
2138                     (loop for v in (rmap-fns name)
2139                           collecting (ratreduce v (rmap-denom name)))))
2140
2141
2142
2143 (defun new-rmap-p (rmap)
2144   (and (eq (car rmap) 'rmap)
2145        (loop for v in (rmap-fns rmap)
2146              when (not (rational-functionp v))
2147              do (return nil) finally (return t))))
2148
2149 (defmacro new-rmap (f)
2150   `(cond ((not (new-rmap-p ,f))(format t "~%Converting an rmap")
2151           (setq ,f (convert-rmap-to-new ,f)))))
2152
2153
2154 (defun my-pairlis (l g)
2155   (loop for v in l
2156         for w in g
2157         collecting (cons v w)))
2158
2159 (defremember compose-rmap (f g)
2160   (let (fns-f subs)
2161     (new-rmap f) (new-rmap g)
2162     (setq subs (loop for gg in (rmap-fns g)
2163                   for v in *xxx* collecting (cons v gg)))
2164     (setq fns-f (loop for ff in (rmap-fns f) collecting (simple-rat-sublis subs ff)))
2165     (construct-rmap fns-f)))
2166
2167 ;;should take [(x1+x2)
2168
2169 (defun construct-rmap (list-funs &aux (answ 1) rat-fns)
2170   (cond (($listp list-funs)(setq list-funs (cdr list-funs))))
2171   (setq rat-fns (loop for v in  list-funs
2172                       when (affine-polynomialp v) collecting (cons v 1)
2173                       else when (rational-functionp v) collecting v
2174                       else when (get v 'disrep) collecting (cons (list v 1 1) 1)
2175                       else collecting (new-rat v)))
2176   (loop for w in rat-fns
2177         do (setq answ (nplcm answ (denom w))))
2178   (make-rmap :fns rat-fns :denom answ))
2179
2180
2181
2182 (defun describe-rmap (rmap)
2183   (cond ((eq (car rmap) 'rmap)
2184          (loop for v in (rmap-fns rmap)
2185                collecting  (header-poly  v) into tem
2186                finally (displa (cons '(mlist) tem)) (format t "Common denom is ..")
2187                (displa ($factor (new-disrep (rmap-denom rmap))))))))
2188
2189 (defvar *give-coordinates* t)
2190
2191 (defun describe-zopen (op)
2192   (cond ((zopen-history op)
2193          (format t "~%The opens history is ~A " (zopen-history op))))
2194   (cond (*give-coordinates*
2195            (format t "~%Zopen with Coord and inverse")
2196            (describe-rmap (zopen-coord op))
2197                               (describe-rmap (zopen-inv op))))
2198   (format t "~%inequality is ..")
2199   (displa ($factor (new-disrep  (zopen-inequality op)))))
2200
2201 (defun xxx (i)
2202   (list (nth (1- i) *xxx* ) 1 1))
2203
2204 ;; I is the index of the special slot so this will return
2205 ;;the ith cover of the blowup of the FIRSTK coordinates of
2206 ;; DIM-affine space
2207 (defun ichart (i firstk dim &aux fns qss pss)
2208   (setq fns
2209     (loop for j from 1 to dim
2210           when (and (<= j firstk) (not (eql j i)))
2211           collecting  (xxx j)
2212           else
2213           collecting   (ptimes (xxx i) (xxx j))))
2214   (setq fns (loop for v in fns
2215                   with den = (xxx i)
2216                   collecting (ratreduce v den)))
2217   (setq pss  (construct-rmap fns ))
2218   (setq fns
2219         (loop for j from 1 to dim
2220               when (and (<= j firstk) (not (eql i j)))
2221               collecting  (ptimes (xxx i) (xxx j))
2222               else collecting  (xxx j) ))
2223   (setq qss (construct-rmap fns ))
2224   (make-zopen :coord pss :inv qss  :inequality 1))
2225
2226 (defun des (obj &optional (stream *standard-output*))
2227   (let ((*standard-output* stream))
2228     (cond
2229       ((atom obj)(format t "~s" obj))
2230       ((or (affine-polynomialp obj)(rational-functionp obj))(sh obj))
2231       (t (case (car obj)
2232            (zopen (show (zopen-history obj)))
2233 ;           (describe-zopen obj))
2234            (rmap (format t "~%Rmap with coordinate functions") (describe-rmap obj))
2235            (ldata (format t "~%Ldata with the comoponent equations")(describe-ldata obj))
2236            (pre-ldata-sheaves (format t "~%Pre Ldata sheaves")(describe-pls obj))
2237            (components (describe-components obj))
2238            (s-var (format t "~%S Variety with basic open sets :")(des (cdr obj) stream))
2239
2240            (t (cond ((macsyma-typep obj) (string-grind obj :stream stream))
2241                     (t
2242                      (loop for v in obj
2243                     for i from 0
2244                     do
2245                     (cond ((and (listp v)(member (car v) '(ldata zopen) :test #'eq))
2246                            (format t "~%Number ~D :" i)))
2247                     (des v stream))))))))
2248     obj))
2249
2250 (defun des (obj &optional (stream *standard-output*))
2251   (let ((*standard-output* stream))
2252     (cond
2253       ((atom obj)(format t "~s" obj))
2254       ((or (affine-polynomialp obj)(rational-functionp obj))(sh obj))
2255       (t (case (car obj)
2256            (zopen
2257             (describe-zopen obj))
2258            (rmap (format t "~%Rmap with coordinate functions") (describe-rmap obj))
2259            (ldata (format t "~%Ldata with the comoponent equations")(describe-ldata obj))
2260            (pre-ldata-sheaves (format t "~%Pre Ldata sheaves")(describe-pls obj))
2261            (components (describe-components obj))
2262            (s-var (format t "~%S Variety with basic open sets :")(des (cdr obj) stream))
2263
2264            (t (cond ((macsyma-typep obj) (string-grind obj :stream stream))
2265                     (t
2266                      (loop for v in obj
2267                     for i from 0
2268                     do
2269                     (cond ((and (listp v) (member (car v) '(ldata zopen) :test #'eq))
2270                            (format t "Number ~D :" i)))
2271                     (des v stream))))))))
2272     obj))
2273
2274
2275
2276 (eval-when
2277     (:compile-toplevel :load-toplevel :execute)
2278   (defmacro pls-opens (pls) `(sv-zopens (pls-s-var ,pls))))
2279
2280 (defvar *reorder-eqns* t)
2281
2282 (defun describe-pls (pls &aux dim eqns  ch-set)
2283   (format t "~%It has ~D opens having ~A components."
2284           (length (pls-opens pls)) (loop for u in (pls-data pls)
2285                                      collecting (length u)))
2286   (loop for v in (sv-zopens (pls-s-var pls))
2287         for w in (pls-data pls)
2288         for i from 0
2289         do  (des v)
2290         (setq dim (length (rmap-fns (zopen-coord v))))
2291         (loop for u in w
2292               when   *reorder-eqns*
2293               when
2294               (progn (multiple-value (eqns ch-set) (order-equations (ldata-eqns u)))
2295                      ch-set)
2296               collecting (- dim (length (ldata-eqns u))) into tem
2297               else collecting "?" into tem
2298               when *reorder-eqns*
2299               collecting (make-ldata :eqns eqns :inequality
2300                                      (ldata-inequality u)
2301                                      :usedup (ldata-usedup u))
2302               into newl
2303               finally
2304               (format t "~%with the ~D ldata" (length w))
2305               (cond (*reorder-eqns*         (format t " of dimension ~A "  tem)
2306                                             (setq w newl)))
2307               (format t "on  open number ~D: "i))
2308         (des w)))
2309
2310 (defun make-component-history (pls &key (add-to-open-history t) &aux dim eqns ch-set)
2311   (loop named sue for op in (pls-opens pls)
2312         for lis-dat in (pls-data pls)
2313         do
2314         (setq dim (length (rmap-fns (zopen-coord op))))
2315         collecting
2316         (loop for u in lis-dat
2317               when
2318               (progn (multiple-value (eqns ch-set) (order-equations (ldata-eqns u)))
2319                      ch-set)
2320               collecting (- dim (length (ldata-eqns u))) into tem
2321               else collecting "?" into tem
2322               finally (cond (add-to-open-history
2323                              (push (cons 'dimensions tem) (zopen-history op))))
2324               (return (cons 'dimensions tem)))))
2325
2326 (defun zopen-dim (zop)
2327   (length (rmap-fns (zopen-coord zop))))
2328
2329 ;;tried to simplify and write better.
2330 ;(defun new-ldata-simplifications (ldata &key (open-g 1)
2331 ;                                 &aux fns used-up answ var  *clear-above*)
2332 ;  (prog
2333 ;    sue
2334 ; nil
2335 ;    (setq fns (copy-list (ldata-eqns ldata)))
2336 ; eliminate-linears
2337 ;    (loop for f in fns
2338 ;
2339 ;         when (and        (not (member f used-up)) (setq var (any-linearp f open-g)))
2340 ;         do
2341 ;         (setq fns  (replace-functions fns
2342 ;                                       f var))
2343 ;         (setq used-up  (replace-functions used-up
2344 ;                                           f var))
2345 ;         (push f used-up)
2346 ;         (go eliminate-linears))
2347 ;    (setq fns (union-equal fns))
2348 ; eliminate-invertible-leading
2349 ;    (loop for f in fns
2350 ;         when (setq var (any-invertible-leading-coefficient f open-g))
2351 ;         when (not (loop for v in *clear-above*
2352 ;                         with deg = (pdegree f var)
2353 ;                         when (and (eq (car v) var) (>= deg (second v)))
2354 ;                         do (return nil) finally (return t)))
2355 ;         do
2356 ;         (setq fns  (replace-functions fns
2357 ;                                       f var))
2358 ;         (setq fns (cons f fns))
2359 ;         (go eliminate-linears))
2360 ;
2361 ;    (setq ldata (copy-structure ldata ldata- eqns (delete 0 (union-equal used-up fns))))
2362 ; make-dichotomy
2363 ;    (setq answ (new-make-dichotomy ldata :open-g open-g))
2364 ;    ;;if no dichotomy continue
2365 ;    (cond ((eql (length answ) 1) (setq ldata (car answ)))
2366 ;         (t (return (delete-redundant-ldata answ))))
2367 ; divide-dichotomy
2368 ;    (setq answ (new-divide-dichotomy ldata :open-g open-g))
2369 ;    (cond ((eql (length answ) 1) (setq ldata (car answ)))
2370 ;         (t (return (delete-redundant-ldata answ))))
2371 ;    (return answ)))
2372 ;
2373 ;
2374 ;;;not sure about the *refine-opens* = nil mode working.
2375 ;(defun new-MAKE-DICHOTOMY (ldata &key (open-g 1)&aux all-facs stop-simplify  eqns-modv ld  answ gg dich lin-dich)
2376 ;  "If stop-simplify is true then it only works if have linear dichotomy.  It returns
2377 ;  a list of  ldata "
2378 ;  (setq dich (find-good-dichotomy  ldata))
2379 ;  (setq all-facs *all-factors*)
2380 ;  (setq dich (order-dichotomy dich))
2381 ;  (show dich)
2382 ;  (setq gg (nplcm open-g (ldata-inequality ldata)))
2383 ;  (cond ((null *refine-opens*) (setq open-g gg) ))
2384 ;  (setq lin-dich
2385 ;       (loop for v in dich when (not (any-linearp v gg)) do (return nil)
2386 ;             finally (return (and dich t))))
2387 ;  (show lin-dich)
2388 ;;; I think the gm-prepared business should all be done after.
2389 ;;; this may not be true.  The simplifications from finding a gm-prepared
2390 ;;; and performing the elimination of variables might be necessary.
2391 ;  (cond ((null lin-dich)
2392 ;        (check-for-gm-prepared (ldata-eqns ldata) open-g)
2393 ;        (show *stop-simplify*)))
2394 ;  (setq stop-simplify *stop-simplify*)
2395 ;  (cond ((and ;;(null *refine-opens*)
2396 ;             *stop-simplify* )
2397 ;        (setq answ (ldata-refinement ldata (car *stop-simplify*)
2398 ;                                     (second *stop-simplify*) :inequality open-g))
2399 ;;       (mshow answ)  (format t "**Is the refinement ")
2400 ;        (setq *stop-simplify* nil)
2401 ;        (setq answ (loop for v in answ
2402 ;                                         appending (new-ldata-simplifications v :open-g
2403 ;                                                                              open-g))))
2404 ;  ;;priority 1 Linear-dichotomy
2405 ;  ;;         2 gm-prepared equation
2406 ;  ;;         3 any dichotomy
2407 ;  ;;proceed with dich if dich is linear or if found no gm-prepared
2408 ;       ((or lin-dich (null *stop-simplify*))
2409 ;        (cond (dich
2410 ;               (loop for v in dich
2411 ;                     with so-far = 1
2412 ;                     appending
2413 ;                     (progn
2414 ;                       (setq eqns-modv
2415 ;                             (loop for facs in all-facs
2416 ;                                   when (not (member v facs))
2417 ;                                   collecting
2418 ;                                   (apply 'gen-ptimes (loop for ter in facs by #'cddr
2419 ;                                                            when (not
2420 ;                                                                   (may-invertp
2421 ;                                                                     ter so-far))
2422 ;                                                            collecting ter))))
2423 ;                       (cond ((member nil eqns-modv :test #'eq) (merror "nil should not be here")))
2424 ;                       (cond ((eq v nil) (merror "nil should not be here")))
2425 ;                       (setq ld (make-ldata))
2426 ;                       (setf (ldata-eqns ld) (cons v eqns-modv))
2427 ;                       (setf (ldata-inequality ld)(nplcm gg so-far))
2428 ;                       (setf so-far (nplcm so-far v))
2429 ;                       (cond ((Null *stop-simplify*)
2430 ;                              (new-LDATA-SIMPLIFICATIONS ld :open-g open-g))
2431 ;                             (t (list  ld))))
2432 ;                     into list-of-ld
2433 ;                     finally
2434 ;                     (cond ((null list-of-ld)
2435 ;                            (setq answ (list (make-ldata eqns '(1) inequality 1))))
2436 ;                                                    (t (setq answ list-of-ld)))))
2437 ;              (t (setq answ (list ldata))))))
2438 ;
2439 ;  (cond ((null answ) (setq answ (list ldata))))
2440 ;;;  (cond ((null answ) (setq answ (list (make-ldata :eqns '(1) :inequality 1)))))
2441 ;;  (setq answ (delete-redundant-ldata answ)))
2442 ;  answ)
2443 ;
2444 ;(defun new-divide-dichotomy (ldata &key (open-g 1) &aux (eqns (ldata-eqns ldata)) f new-eqns do-dichotomy
2445 ;       try                  occurs vars highest-vars orig-rep repeat eqns-rep)
2446 ;
2447 ;  "endeavors to turn ldata into an  triangular list of eqns
2448 ; so that each equation has  possibly one more variable occurring than the
2449 ; previous.  It takes a good order for the variables and then takes the first variable
2450 ; to be highest in two succeeding eqns, and does a division to try to correct this.  If
2451 ; the leading variable is not invertible it does a dichotomy.   This dichotomy must be
2452 ; at the level of open sets, since otherwise we will not get component containment."
2453 ;  ;;ordering should take into account open-g
2454 ;  (setq vars  (good-order-variables eqns))
2455 ;  (setq occurs (second vars))
2456 ;  (setq vars (first vars))
2457 ;  (setq highest-vars
2458 ;       (loop for v in occurs
2459 ;             collecting (loop for u in vars
2460 ;                              when (member u v :test #'eq)
2461 ;                              do (return u))))
2462 ;  (show vars highest-vars)
2463 ;  (setq repeat
2464 ;       (loop named rep for v in vars
2465 ;             do
2466 ;             (loop for w on highest-vars
2467 ;                   when (and (eq (car w) v)
2468 ;                             (member v (cdr w) :test #'eq))
2469 ;                   do (return-from rep v))))
2470 ;  (cond
2471 ;    (repeat
2472 ;     (setq eqns-rep  (loop for v in eqns
2473 ;                          for u in highest-vars
2474 ;                          when (eq u repeat)
2475 ;                          collecting v))
2476 ;     (setq orig-rep (copy-list eqns-rep))
2477 ;     (show (length eqns-rep))
2478 ;     (setq eqns-rep (sort-key eqns-rep '< 'pdegree repeat))
2479 ;     (loop for v on eqns-rep while (>= (length v ) 2)
2480 ;           do
2481 ;          (cond ((eq ( pdegree (first v) repeat)
2482 ;                     (pdegree (second v) repeat))
2483 ;                 ;;choose the least complex leading coefficient to divide by
2484 ;                 (setq try
2485 ;                       (sort-key (firstn 2 v) '<
2486 ;                                 #'(lambda (u var) (pcomplexity
2487 ;                                                     (leading-coefficient u var)))
2488 ;                                 repeat)))
2489 ;                (t (setq try (firstn 2 v))))
2490 ;          (setq f (second try))
2491 ;          (multiple-value-bind (rem c-reqd)
2492 ;              (gen-prem  f (first try) repeat)
2493 ;            (show rem)
2494 ;            (cond ((null rem) (merror "empty")))
2495 ;            (cond ((may-invertp c-reqd open-g)
2496 ;                   (return (new-ldata-simplifications
2497 ;                             (copy-structure ldata ldata-
2498 ;                                             eqns
2499 ;                                             (cons rem (delete f(copy-list  eqns))))
2500 ;                             :open-g open-g)))))
2501 ;          finally (return (list ldata))))))
2502 ;
2503 ;
2504 ;
2505
2506 ;resolution(f):=
2507 ;block([g,h],
2508 ;      for i thru 5   do
2509 ;      (g: coeff(f,x),
2510 ;       h: f-x*g,
2511 ;       f: expand(subst(x+g/2,x,h)+g^2/4)),
2512 ;       return(f));
2513 ;(defun resolution (f n var &aux h g gg ggg)
2514 ;  (loop for i below n
2515 ;       with mon = (list var 1 1)
2516 ;       do (setq g  (pcoeff (num f) mon))
2517 ;       (setq h (ratdifference f (setq gg (ratreduce (ptimes g mon) (denom f)))))
2518 ;       (setq ff (gen-rat-sublis
2519 ;                 (list var)
2520 ;                 (list (setq ggg (ratplus (cons mon 1) (ratreduce g (* 2 (denom f))))))
2521 ;                 h))
2522 ;       (setq g (ratreduce g (denom f)))
2523 ;       (setq ff  (ratplus ff (rattimes (cons 1 4)  (rattimes g g nil) t)))
2524 ;       (show (denom h))
2525 ;       (setq f ff)
2526 ;       finally (return f)))
2527
2528
2529 ;(defun new-solve-ldata (ldata &key (open-g 1))
2530 ;  (user:function-let
2531 ;    (ldata-simplifications #'(lambda (&rest l) (list (car l))))
2532 ;
2533 ;   (loop until (equal list-ld prev-ld)
2534 ;        with list-ld = (list ldata)
2535 ;        do
2536 ;     (format t "~%make-dichotomy")
2537 ;     (setq list-ld
2538 ;          (loop for ld in list-ld
2539 ;                appending (make-dichotomy ld :open-g open-g)))
2540 ;     (mshow list-ld)
2541 ;     (format t "~%divide-dichotomy")
2542 ;     (setq list-ld
2543 ;          (loop for ld in list-ld
2544 ;                appending (divide-dichotomy ld :open-g open-g)))
2545 ;     (mshow list-ld))))
2546 ;
2547 ;
2548 ;;;ideas have a simpflag slot
2549 ;;;which contains things like 'all-irreducible 'no-linears 'no-
2550 ;;;a clear above slot which has '(x 3 y 1) to indicate that only one
2551 ;;;note if z+y is a polynomial can't use it for z and y so it has to
2552 ;;;be in the used up group.