share/affine/sheafc.lisp

   1 ;;; -*- Mode:Lisp; Package:CL-MAXIMA; Syntax:COMMON-LISP; Base:10 -*-
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;                                                                    ;;;;;
   4 ;;;     Copyright (c) 1984 by William Schelter,University of Texas     ;;;;;
   5 ;;;     All rights reserved                                            ;;;;;
   6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   7
   8 (in-package :maxima)
   9
  10 ;(defun grobner-subset (id1 id2 &optional ( may-invert 1) &aux answ unit-ideal)
  11 ;  (declare   (values  answ unit-ideal))
  12 ;  (cond ((not (eql 1 may-invert))
  13 ;        (setq *poly-simplifications* (grobner-remember (cons (pdifference
  14 ;                                                               (ptimes (st-rat '$zzzzz)
  15 ;                                                                       may-invert)
  16 ;                                                               1)
  17 ;                                                             id2))))
  18 ;       (t (setq *poly-simplifications* (grobner-remember id2))))
  19 ;  (setq answ
  20 ;       (cond (($zerop (polysimp 1))(setq unit-ideal 'unit-ideal))
  21 ;             (t
  22 ;              (loop for v in id1 when (not ($zerop (polysimp v)))
  23 ;                    do (return nil)
  24 ;                    finally (return t)) )))
  25 ;  (values answ unit-ideal))
  26
  27 (defun grobner-subset (id1 id2 &optional ( may-invert 1) &aux answ unit-ideal)
  28 ;  (declare   (values  answ unit-ideal))
  29   (setq unit-ideal (unit-idealp id2 may-invert))
  30   (setq answ (loop for v in id1 when (not ($zerop (polysimp v)))
  31                                             do (return nil)
  32                              finally (return t)) )
  33   (values answ unit-ideal))
  34
  35 ;(defun unit-idealp (list-eqns &optional (may-invert 1) &aux )
  36 ;  (cond ((not (eql 1 may-invert))
  37 ;        (setq *poly-simplifications* (grobner-remember (cons (pdifference
  38 ;                                                               (ptimes (st-rat '$zzzzz)
  39 ;                                                                       may-invert)
  40 ;                                                               1)
  41 ;                                                             list-eqns))))
  42 ;       (t (setq *poly-simplifications* (grobner-remember list-eqns))))
  43 ;  (cond (($zerop (polysimp 1)) 'unit-ideal)))
  44
  45 (defvar *invert-separately* t)
  46
  47 (defun fast-grobner-basis (list-eqns &optional (may-invert 1) (string-for-inverse-name "")
  48                             &aux eqns facts newvars inverses
  49                     (pcomplexity-of-inverses 80))
  50 ;  (declare    (values   *poly-simplifications* newvars))
  51    (cond
  52     ((loop for v in list-eqns
  53            when  (may-invertp v may-invert)
  54            do (setq *poly-simplifications* (list 1 (rzero)))  (return 'done)))
  55     ((not (numberp may-invert))
  56      (cond ((and *invert-separately* (> (pcomplexity may-invert)
  57                                         pcomplexity-of-inverses))
  58             (setq facts (npfactor may-invert)))
  59            (t (setq facts (list (square-free may-invert) 1))))
  60      (setq newvars (unused-variables (truncate (length facts) 2) string-for-inverse-name))
  61      (setq inverses
  62            (loop for (pol deg) on facts by #'cddr
  63                  for new in newvars
  64                  when (< (pcomplexity pol) pcomplexity-of-inverses)
  65                  collecting
  66                  (pdifference (ptimes (list new 1 1)
  67                                       pol) 1)
  68                  else do
  69                  (format t "~3%*******  That inverse was a little to complex to use: .. ****~3%")
  70                  (des pol)))
  71      (setq eqns (sort (copy-list (append inverses  list-eqns)) 'alphalessp))
  72      (setq *poly-simplifications* (grobner-remember eqns)))
  73     (t (setq *poly-simplifications* (grobner-remember (sort (copy-list list-eqns) 'alphalessp)))))
  74    (values   *poly-simplifications* newvars))
  75
  76
  77 (defvar *invalid-contractions* nil)
  78 (defun contract-ideal-localization (list-eqns may-invert &aux simps newvars)
  79     (catch 'took-too-long
  80       (multiple-value-setq (simps newvars) (fast-grobner-basis list-eqns may-invert 'zzzzzzzz)))
  81     (cond ((null simps) (format t "~%******Invalid contraction")
  82            (push (list list-eqns may-invert) *invalid-contractions*) list-eqns)
  83           (t
  84     (loop for v in (poly-relations-from-simplifications)
  85           when (loop for var in newvars
  86                      when (not (zerop (pdegree v var)))
  87                      do (return nil) finally (return t))
  88           collecting v))))
  89
  90 (defvar *too-long-for-grobner-basis* nil)
  91 (defun unit-idealp (list-eqns &optional (may-invert 1) &aux answer)
  92
  93   "Has a side effect of setting the *poly-simplifications* to be the correct ones"
  94   (format t "~%Entering unit-ideal")
  95   (if-verbose (show (gen-pcomplexity list-eqns) (pcomplexity may-invert)))
  96   (cond
  97     ;;this is handled by grobner-remember now
  98 ;    ((member (list list-eqns may-invert) *too-long-for-grobner-basis*)
  99 ;        (setq *poly-simplifications* nil))
 100         (t
 101   (setq answer (catch 'took-too-long  (fast-grobner-basis list-eqns may-invert "ZZ")))
 102   (cond ((eq answer 'took-too-long) (push (list list-eqns may-invert) *too-long-for-grobner-basis*)
 103                                           (setq *poly-simplifications* nil))
 104         (($zerop (polysimp 1)) 'unit-ideal)))))
 105
 106 (defun ptruncate (poly var deg)
 107   (cond ((<= (pdegree poly var) deg)
 108          poly)
 109         (t (ptruncate1 poly var deg))))
 110 (defun ptruncate1 (poly var degree &aux answ tem)
 111   (cond ((atom poly) poly)
 112         ((pointergp (p-var poly) var)
 113          (loop for (deg cof) on (cdr poly) by #'cddr
 114                when (not (pzerop (setq tem
 115                                        (ptruncate1 cof var degree))))
 116                collecting deg into pol
 117                and
 118                collecting tem into pol
 119                finally (return (cond (pol (cons (p-var poly) pol))
 120                                      (t 0)))))
 121         ((eq (p-var poly) var)
 122          (loop for v on (cdr poly) by #'cddr
 123                when (<= (car v) degree)
 124                do (return (setq answ v)))
 125          (cond (answ (cons var answ))
 126                (t 0)))
 127         (t poly)))
 128
 129 (defun sort-key (data pred key &rest rest-key-arg)
 130   (declare (special rest-key-arg key pred))
 131   (sort (copy-list data) #'(lambda (u v)
 132                  (funcall pred (apply key u rest-key-arg)
 133                           (apply key v rest-key-arg)))))
 134
 135 (defun leading-coefficient (poly var)
 136   (pcoeff poly (list var (pdegree poly var) 1)))
 137
 138
 139 ;;condition on domain of op1-->op2
 140 ;; is that the  part of the denom of the map op1-->op2  which is relatively
 141 ;;prime to the inequality on op1, should be invertible on op2.
 142 ;;condition on the range is that the image of a point in op1
 143 ;;should satisfy the inequality of op2
 144 (defremember open-subsetp (op1 op2 )
 145   (let ( map12 map21 gg gf domain-ok)
 146     (setq map21  (find-ring-map op2 op1))
 147     (setq gf (numerator (apply-rmap map12 (rmap-denom map21))))
 148     (setq map12 (find-ring-map op1 op2))
 149     (setq gg (numerator (apply-rmap map12 (zopen-inequality op1))))
 150     (cond ((may-invertp gf (ptimes gg (zopen-inequality op2))) (setq domain-ok t))
 151           (t nil))
 152     ;;check the image lands in
 153     (cond ((and domain-ok
 154                 (may-invertp (numerator (apply-rmap map21 (zopen-inequality op2)))
 155                              (zopen-inequality op1))) t)
 156           (t nil))))
 157
 158 ;
 159 ;  (setq gg (nplcm (numerator  (apply-rmap map (zopen-inequality op2) )) (rmap-denom map)))
 160 ;  (cond ((may-invertp gg (zopen-inequality op1)) t)
 161 ;       (t nil)))
 162 ;;works but makes to many refinements
 163 (defun open-refinement-for-list (zopen list-eqns
 164                                  &aux ( gg  (zopen-inequality zopen)) tem answ )
 165   (format t "Should be using best open cover see file")
 166   (cond ((null list-eqns)(list zopen))
 167         (t
 168          (loop for pol on list-eqns
 169                when (and (not (any-linearp (car pol)gg))
 170                          (setq tem (gm-prepared (car pol) :inequal gg)))
 171                do (setq answ (open-refinement zopen (car pol)  (length tem)))
 172                (return (loop for vv in answ
 173                              appending (open-refinement-for-list vv (cdr pol))))
 174                finally (return (list zopen))))))
 175
 176
 177 ;;;
 178 (defun open-refinement-for-list (zopen list-eqns
 179                                  &aux ( gg  (zopen-inequality zopen)) tem answ  refs)
 180   (cond ((null list-eqns)(list zopen))
 181         (t
 182          (loop for pol in list-eqns
 183                when (and (not (any-linearp  pol gg))
 184                          (setq tem (gm-prepared  pol :inequal gg)))
 185                do
 186                (setq answ (best-open-cover list-eqns (zopen-inequality zopen)))
 187                (iassert (unit-idealp answ (zopen-inequality zopen)))
 188                (loop for gg in answ
 189                      when (not (or
 190                                  (unit-idealp list-eqns gg)
 191                                  (linear-ldatap (make-ldata :eqns list-eqns):open-g gg)))
 192                      do (fsignal "This refinement is not sufficient."))
 193                ;;what here????
 194
 195                (format t "~%Refining for the equations ~/maxima::tilde-q-fsh/ with inequality ~/maxima::tilde-q-fsh/ ~
 196                     ~%Using inequalities ~/maxima::tilde-q-fsh/ " list-eqns gg answ)
 197
 198                (setq refs
 199                      (loop for gg in answ
 200                            collecting (zl-copy-structure zopen zopen-
 201                                                       inequality gg)))
 202                ;;could stick something recursive in here to make linear variables in each
 203                ;;eqn disjoint
 204 ;             (return (loop for vv in refs
 205 ;                           appending (open-refinement-for-list vv list-eqns)))
 206                (return refs)
 207                finally (return (list zopen))))))
 208
 209
 210
 211
 212 ;(defun eliminate-larger (a-list &key ( key 'identity) test test-not (destructive nil)
 213 ;                        &aux maybe-not)
 214 ;  (cond ((null destructive)(setq a-list (copy-list a-list))))
 215 ;  (cond (test-not (setq test test-not)
 216 ;                 (setq maybe-not 'null))
 217 ;       (t (setq maybe-not 'identity)))
 218 ;  (loop for v in a-list
 219 ;             for i from 0
 220 ;             appending (loop for  w in a-list
 221 ;                             for j from 0
 222 ;                             when (and (not (eql i j))
 223 ;                                       (funcall maybe-not (funcall test  (funcall key v)
 224 ;                                                                (funcall key w))))
 225 ;                             do (setq a-list (delete w a-list))))
 226 ;  a-list)
 227
 228 ;;above does not seem to work
 229 (defun eliminate-larger (a-list &key ( key 'identity) test test-not verify
 230                          (destructive nil)
 231                          &aux maybe-not answ orig)
 232   "Eliminates equal duplicates in A-LIST and then if test is given it
 233   eliminates items w such that (funcall test v w) is true for zl-SOME v not equal to w,
 234   and if test-not is specified it eliminates w such that (null (funcall test v w)) holds"
 235   (cond ((null destructive)(setq a-list (copy-list a-list))))
 236   (cond (test-not (setq test test-not)
 237                   (setq maybe-not 'null))
 238         (t (setq maybe-not 'identity)))
 239   ;;in case two things are equal we don't want to delete both!!
 240   (setq orig a-list)
 241   (setq a-list  (union-equal a-list))
 242   (setq answ (copy-list a-list))
 243   (loop for v in a-list
 244      for i from 0
 245      do (loop for w in a-list
 246            for j from 0
 247            when (and (not (eql i j))
 248                      (funcall maybe-not (funcall test  (funcall key v)
 249                                                  (funcall key w))))
 250            do
 251              (setq answ (delete w answ :test #'equal))))
 252   (cond (verify
 253          (loop for v in orig
 254                 do (loop for w in answ
 255                           when
 256                           (funcall maybe-not (funcall test  (funcall key w)
 257                                                       (funcall key v)))
 258
 259                           do (return nil)
 260                           finally (merror "~A is not larger than anything in the answer ~a" v answ)))))
 261   answ)
 262
 263
 264 (defun eliminate-smaller (a-list &key ( key 'identity) test test-not verify (destructive nil)
 265                           &aux maybe-not answ orig)
 266   "Eliminates equal duplicates in A-LIST and then if test is given it
 267   eliminates items w such that (funcall test v w) is true for zl-SOME v not equal to w,
 268   and if test-not is specified it eliminates w such that (null (funcall test v w)) holds"
 269   (cond ((null destructive)(setq a-list (copy-list a-list))))
 270   (cond (test-not (setq test test-not)
 271                   (setq maybe-not 'null))
 272         (t (setq maybe-not 'identity)))
 273   ;;in case two things are equal we don't want to delete both!!
 274   (setq orig a-list)
 275   (setq a-list  (union-equal a-list))
 276   (setq answ (copy-list a-list))
 277   (loop for v in a-list
 278      for i from 0
 279      do (loop for  w in a-list
 280            for j from 0
 281            when (and (not (eql i j))
 282                      (funcall maybe-not (funcall test  (funcall key v)
 283                                                  (funcall key w))))
 284            do
 285              (setq answ (delete v answ :test #'equal))))
 286   (cond (verify
 287          (loop for v in orig
 288             do (loop for w in answ
 289                   when
 290                     (funcall maybe-not (funcall test  (funcall key w)
 291                                                 (funcall key v)))
 292
 293                   do (return nil)
 294                   finally (merror "~A is not larger than anything in the answer ~a" v answ)))))
 295   answ)
 296
 297 (defun add-pls-zopen-history (pls)
 298   (loop for v in (pls-opens pls)
 299         for i from 0
 300         collecting (add-zopen-history v i) into opens
 301         finally (setf (pls-opens pls ) opens)))
 302
 303 (defun find-position-in-list (item list &key (test #'equal))
 304   (loop for v in list
 305         for i from 0
 306         when (funcall test v item)
 307         do (return i)))
 308
 309 (defun normalize-pls (pls  open ld-number &aux nth-open all-opens-data opens comp
 310                       lis-data MAPL tem  refs codimension-component)
 311   "makes the list-eqns become first coords in open and in the other opens, it
 312   should also look after refining if necessary: ie. if one of eqns is gm-prepared
 313   it must refine"
 314 ;  (declare (values normal-pls opens-not-to-blowup component-codimension))
 315   (setq opens (pls-opens pls))
 316   (setq lis-data  (pls-data pls))
 317   (cond ((numberp open)(setq nth-open open))
 318        (t (setq nth-open(find-position-in-list open opens))))
 319   (setq comp (match-components pls nth-open :ldata-number ld-number))
 320   (setq codimension-component (length (ldata-eqns (nth ld-number
 321                                                        (nth nth-open lis-data)))))
 322   (setq comp
 323         (loop for v in comp
 324         when (eq (cdr (assoc nth-open v :test #'equal)) ld-number)
 325         do (return v)))
 326
 327   (loop for op  in opens
 328         for lis-ld in lis-data
 329         for i from 0
 330         when (setq ld-number  (cdr (assoc i comp :test #'equal)))
 331         do
 332         (setq refs (open-refinement-for-list op (ldata-eqns (nth ld-number lis-ld))))
 333         (loop for op1 in refs
 334               do (push (list op1  lis-ld (ldata-eqns (nth ld-number lis-ld))) all-opens-data))
 335         else
 336         do  (push (list op  lis-ld nil) all-opens-data))
 337 ;  (setq all-opens-data  (eliminate-larger all-opens-data :key 'car
 338 ;                                         :test 'open-subsetp ))
 339   (setq all-opens-data (nreverse  all-opens-data))
 340   (loop for v in all-opens-data
 341         for i from 0
 342         when (third v)
 343         collecting (setq tem (normalize-zopen (car v ) (third v))) into all-opens
 344         and
 345         do (setq MAPL (find-ring-map (car v) tem))
 346         and
 347         collecting (loop for ld in (second v)
 348                          collecting (apply-rmap MAPL ld))
 349         into data
 350         else
 351         collecting i into opens-not-to-blowup
 352         and
 353         collecting (car v) into all-opens
 354         and collecting (second v) into data
 355         finally (return (values
 356                           (make-pre-ldata-sheaves
 357                                      :s-var      (make-s-var :zopens  all-opens)
 358                                      :data  data)
 359                           opens-not-to-blowup
 360                           codimension-component))))
 361
 362 (defun a-factor (fact-list1 fact-list2 &aux tem)
 363   (loop for (pol deg) on fact-list1 by #'cddr
 364         when (not  (and (setq tem (get-alt pol fact-list2))
 365                         (>= tem deg)))
 366         do (return nil)
 367         finally (return t)))
 368
 369 (defun multiply-out-factorization (facts)
 370   (loop for (pol deg) on facts by #'cddr
 371      with answ = 1
 372      do (setq answ (ptimes answ (pexpt pol deg)))
 373      finally (return answ)))
 374
 375 (defun monomialp (pol)
 376   (cond((atom pol) t)
 377        (t (and (eql (length pol) 3) (monomialp (third pol))))))
 378 (defun one-variablep (f)
 379   (and (eql (length f) 3) (numberp (third f))))
 380
 381 (defvar *leave-unit-ldata* nil)
 382 (defun ldata-simplify-homogeneous (ldata &optional (may-invert 1) &aux all-others used  eqns facts answ simpler (added t) eqn tem)
 383   (process-sleep 15)
 384   (setq facts (loop for v in (ldata-eqns ldata)
 385
 386                     do   (process-sleep 15)
 387                     when (not ($zerop v))
 388                     collecting (setq tem  (non-constant-factors v may-invert))
 389                     ))
 390                     ;;when (null tem) do (mshow v) (break 'here)))
 391   (loop for v in facts
 392         when (null v)
 393         do
 394         (return (setq answ 'unit)))
 395   (cond
 396     ((null answ)  (setq facts (eliminate-larger facts :test 'a-factor :verify t))
 397      (loop for v in facts
 398            when (loop for  (pol deg) on v by #'cddr
 399                       when (not (one-variablep pol))
 400                       do (return nil)
 401                       finally (return t))
 402
 403            collecting (multiply-out-factorization v) into tem
 404            else collecting (multiply-out-factorization v) into others
 405            finally (setq used tem) (grobner-basis tem) (setq all-others others))
 406      (loop named sue while added
 407            for ii from 0
 408            do
 409            (setq added nil)
 410            (loop for w in all-others
 411                  do (setq eqn (numerator (polysimp w)))
 412                  when (not ($zerop eqn))
 413                  do
 414                  (cond ((numberp eqn)   (setq answ 'unit) (return-from sue 'done))
 415                        ((monomialp eqn)
 416                         (add-to-poly-simplifications eqn) (push eqn used)
 417                         (setq added t) (return 'again)))
 418                  and
 419                  collecting eqn into tem
 420                  when (may-invertp eqn may-invert) do (return-from sue (setq answ 'unit))
 421                  finally (setq simpler tem))
 422            finally (setq eqns (append  (delete 0 simpler :test #'equal) used)))))
 423   (setq answ  (cond ((eq answ 'unit) (format t "~%It turned out to be trivial")
 424                      (cond (*leave-unit-ldata* (list (make-ldata :eqns '(1))) )))
 425                     (t (format t "~%There are ~A equations of complexity ~A of which ~A ~
 426                                    are monomials"
 427                                (length eqns) (gen-pcomplexity eqns) (length used))
 428                      (list (zl-copy-structure ldata ldata- eqns eqns))))))
 429
 430
 431 ;;;try to do without factoring.
 432 ;(defun ldata-simplify-homogeneous (ldata &optional (may-invert 1)
 433 ;                                  &aux all-others used  eqns facts answ)
 434 ;  (process-sleep 15)
 435 ;  (setq facts (loop for v in (ldata-eqns ldata)
 436 ;                   do   (process-sleep 15)
 437 ;                   collecting (non-constant-factors v may-invert)))
 438 ;  (loop for v in facts
 439 ;       when (null v)
 440 ;       do (return (setq answ 'unit)))
 441 ;  (cond ((null answ)  (setq facts (eliminate-larger facts :test 'a-factor))
 442 ;        (loop for v in facts
 443 ;              when (loop for  (pol deg) on v by #'cddr
 444 ;                         when (not (one-variablep pol))
 445 ;                         do (return nil)
 446 ;                         finally (return t))
 447 ;
 448 ;              collecting (multiply-out-factorization v) into tem
 449 ;              else collecting (multiply-out-factorization v) into others
 450 ;              finally (setq used tem) (grobner-basis tem) (setq all-others others))
 451 ;        (loop for w in all-others
 452 ;              collecting (numerator (polysimp w)) into tem
 453 ;              finally (setq eqns (append  (delete 0 tem) used)))))
 454 ;  (cond ((eq answ 'unit) (format t "~%It turned out to be trivial")nil )
 455 ;       (t
 456 ;        (list (copy-structure ldata ldata- eqns eqns)))))
 457
 458
 459
 460
 461 ;;if the ldata on one open are trivial it will give empty list corresponding to that open.
 462 (defun simplify-svar-homogeneous (pls &aux  all-data tem)
 463   (setq pls (copy-list-structure pls))
 464   (setq all-data
 465         (loop
 466            for op in (pls-opens pls)
 467            for w in (pls-data pls)
 468            collecting (cons op
 469                             (loop for ld in w
 470                                when (setq tem (ldata-simplify-homogeneous ld (zopen-inequality op)))
 471                                appending tem))))
 472   (construct-pre-ldata-sheaves  :s-var (make-s-var :zopens (mapcar #'car all-data))
 473                                 :data (mapcar #'cdr all-data)))
 474
 475 (defun intersection-inequality-in-op1 (op1 op2 &aux hh )
 476   (let  ((map1 (find-ring-map op1 op2))
 477          (map2 (find-ring-map op2 op1)))
 478     (setq hh (apply-rmap-to-square-free-factors map2 (zopen-inequality op2)))
 479     (setq hh (sftimes
 480                hh (apply-rmap-to-square-free-factors map2 (rmap-denom map1))))
 481     (setq hh (sftimes  hh (zopen-inequality op1)))
 482     (setq hh (sftimes hh (rmap-denom map2)))))
 483
 484 (defun refine-pls (pls list-eqns &aux refs changed new-pls)
 485   (loop for op in (pls-opens pls)
 486         for lis in list-eqns
 487         for lis-dat in (pls-data pls)
 488         appending
 489         (progn
 490           (setq refs (open-refinement-for-list op lis ))
 491 ;         (loop for ope in refs
 492 ;              (cond ((and lis original-open
 493 ;                       (null (unit-idealp lis (zopen-inequality ope)))
 494 ;                       (null (linear-ldatap
 495 ;                               (make-ldata eqns lis)
 496 ;                               :open-g  (zopen-inequality ope))))
 497 ;                     (if tried (fsignal "this won't work"))
 498 ;                     (setq hh (intersection-inequality-in-op1 ope original-open))
 499 ;                     (setq tem1 (try-harder-and-intersect
 500 ;                                  op
 501 ;                             (make-ldata eqns eqns1) lis-dat hh))
 502 ;                     (setq lis (ldata-eqns tem1)) (setq tried t)
 503 ;                     (format t "~%*****Had to intersect ~VQ with original ldata to obtain ~VQ" eqns1 'fsh tem1 'fsh)
 504
 505                (cond ((> (length refs ) 1)(setq changed t)))
 506                refs)  into all-ops
 507         appending
 508         (make-list (length refs) :initial-element lis-dat  ) into all-dat
 509         appending
 510         (make-list (length refs) :initial-element lis  ) into the-eqns
 511         finally
 512         (setq new-pls(construct-pre-ldata-sheaves :opens all-ops :data all-dat))
 513         (cond (changed (add-pls-zopen-history new-pls)))
 514         (return
 515           (values new-pls
 516                   the-eqns ))))
 517
 518 (defun all-linearp ( eqns &optional (inequal 1))
 519   (loop for v in eqns when (not (any-linearp v inequal))
 520         do (return nil) finally (return t)))
 521
 522 ;(defun prepare-for-blowup (pls open-number eqns &aux tr-data ld from-open ref-pls norm-op
 523 ;                          transl-eqns)
 524 ;  "the eqns will be translated around to other opens and the coords and ldata will be
 525 ; normalized on the components it meets.   A pls and opens-not-to-blowup are returned"
 526 ;  (setq from-open (nth open-number(pls-opens pls)))
 527 ;  (setq ld  (make-ldata eqns eqns))
 528 ;  (loop for op in (pls-opens pls)
 529 ;       for op-number from 0
 530 ;       collecting
 531 ;       (ldata-eqns
 532 ;         (translate-reduced-component open-number op-number ld pls)) into lis-eqns
 533 ;       finally (multiple-value ( ref-pls transl-eqns)
 534 ;                 (refine-pls pls lis-eqns)))
 535 ;  (loop for op in (pls-opens ref-pls)
 536 ;       for eqns in transl-eqns
 537 ;       for i from 0
 538 ;       for lis-dat in (pls-data ref-pls)
 539 ;       when (and eqns lis-dat (all-linearp eqns (zopen-inequality op)))
 540 ;       do
 541 ;       (multiple-value
 542 ;          (norm-op tr-data)
 543 ;             (normalize-zopen op eqns :data lis-dat))
 544 ;       and collecting (cons norm-op tr-data) into all-data
 545 ;       else
 546 ;       collecting (cons op lis-dat) into all-data
 547 ;       and collecting i into opens-not-to-blowup
 548 ;       finally
 549 ;       (return(values (construct-pre-ldata-sheaves :opens (mapcar 'car all-data)
 550 ;                                      :data (mapcar 'cdr all-data))
 551 ;                      opens-not-to-blowup))))
 552 ;;Better: the data of what to blow up may come from the associated reduced sheaf,
 553 ;;so to specify the data we should specify an open and the equations, and the
 554 ;;pls to blowup.
 555
 556 (defun zopen-special-subset (op1 op2)
 557    (and (equal (zopen-coord op1) (zopen-coord op2))
 558         (may-invertp (zopen-inequality op2) (zopen-inequality op1))))
 559 (defun zopen-equal (op1 op2)
 560   (and (equal (zopen-coord op1) (zopen-coord op2))
 561        (equal (zopen-inequality op1) (zopen-inequality op2))))
 562
 563
 564 (defun try-harder-and-intersect (op trans-ld lis-dat hh &aux answ  *refine-opens*)
 565
 566   (setq answ (simplify-ldata (make-ldata :eqns (append (ldata-eqns trans-ld)
 567                                                       (ldata-eqns (car lis-dat))))
 568                              :open-g hh ))
 569   (fsignal 'hi)
 570   (format t "~%*** Having to try harder ***.")
 571   (cond ((equal(length answ) 1) (setq answ (car answ)))
 572         ((null answ) nil)
 573         (t (fsignal "Two components in the translate")))
 574   (cond ((linear-ldatap answ :open-g (zopen-inequality op)) answ)
 575         (t (fsignal "Real trouble : the ldata are not linear after translation"))))
 576
 577 (defun prepare-for-blowup (pls open-for-equations eqns
 578                            &aux tr-data ld from-open ref-pls norm-op orig-eqns
 579                            transl-eqns int-inequal list-eqns tem)
 580   "the eqns will be translated around to other opens and the coords and ldata will be
 581  normalized on the components it meets.   A pls and opens-not-to-blowup are returned"
 582   (setq from-open open-for-equations)
 583   (setq orig-eqns eqns)
 584   (setq pls (copy-list-structure pls))
 585   (cond ((ldatap eqns)(setq ld eqns))
 586         (t   (setq ld  (make-ldata :eqns eqns))))
 587   (loop for op in (pls-opens pls)
 588         for op-number from 0
 589         for lis-dat in (pls-data pls)
 590         do (setq tem nil)
 591         when lis-dat
 592         collecting
 593         (progn
 594           (multiple-value-setq (tem int-inequal)
 595                                (translate-reduced-component-and-reduce
 596                                  from-open op ld
 597                                  :homogeneous-ldata-on-to-open (car lis-dat)) )
 598           (ldata-eqns tem))
 599         into lis-eqns
 600         else collecting nil into lis-eqns
 601         collecting int-inequal into lis-int-inequal
 602         when (and tem (not  (= (length (ldata-eqns tem)) (length (ldata-eqns ld)))))
 603         do (merror "wrong number of  equations")
 604         do
 605         (format t "~2%On open ~A the inequality is ~/maxima::tilde-q-fsh/ (before refinement).~
 606                    ~%The translated equations are ~/maxima::tilde-q-fsh/. "
 607                 op-number (zopen-inequality op) (ldata-eqns tem))
 608 ;       when tem do (fsignal 'hi)
 609
 610         finally (multiple-value ( ref-pls transl-eqns)
 611                   (refine-pls pls (setq list-eqns lis-eqns)))
 612
 613         (setq ref-pls (copy-list-structure ref-pls)))
 614   (loop for op in (pls-opens ref-pls)
 615         for eqns in transl-eqns
 616         for i from 0
 617         for lis-dat in (pls-data ref-pls)
 618         when (and eqns lis-dat (all-linearp eqns (zopen-inequality op)))
 619         do
 620         (multiple-value
 621           (norm-op tr-data)
 622           (normalize-zopen op eqns :data lis-dat))
 623         (format t "~%On this open ~/maxima::tilde-q-fsh/ will correspond to ~/maxima::tilde-q-fsh/" eqns
 624                 (loop for i from 1 to (length eqns) collecting (xxx i)))
 625         (check-normalization open-for-equations norm-op orig-eqns)
 626 ;       (break 'prepare1)
 627         (format t "~%normalizing open ~A . " i)
 628         and collecting (cons norm-op tr-data) into all-data
 629         else
 630         do
 631         (format t "%The open ~A does not meet" i)
 632                                                 ;(break 'prepare)
 633         and
 634         collecting (cons op lis-dat) into all-data
 635         and collecting i into opens-not-to-blowup
 636         finally
 637         (cond ((eql (length opens-not-to-blowup) (length all-data))
 638                (merror "Not going to blowup anything")))
 639         (return(values (construct-pre-ldata-sheaves :opens (mapcar 'car all-data)
 640                                                     :data (mapcar 'cdr all-data))
 641                        opens-not-to-blowup ))))
 642
 643
 644
 645 (defun check-normalization (original-open  norm-op original-eqns &aux answ)
 646   (cond ((not (ldatap  original-eqns))(setq  original-eqns
 647                                         (make-ldata :eqns    original-eqns))))
 648   (setq answ (translate-component-and-reduce original-open norm-op
 649                                   original-eqns))
 650   (setq answ (sort (copy-list (ldata-eqns answ)) 'alphalessp))
 651   (loop for v in answ
 652         for w in *xxx*
 653         when (not (and (eq (first v) w) (eql (second v) 1)))
 654         do (merror "bad normalization")))
 655
 656
 657 (defun hhh  (pls reduced-pls op-num dat-num)
 658   (blowup-pls-and-reduced-pls pls reduced-pls (nth op-num (pls-opens reduced-pls))
 659                               (ldata-eqns (pls-ldata reduced-pls op-num :ldata-number dat-num))))
 660 (defun blowup-pls-and-reduced-pls (pls reduced-pls opens-for-eqns eqns
 661                                    &aux new-pls new-red-pls)
 662   (setq new-pls (blowup-pls pls opens-for-eqns eqns :simplify-homogeneous t))
 663   (setq new-red-pls (blowup-pls reduced-pls opens-for-eqns eqns
 664                                 :throw-out-components-in-exceptional-divisor t))
 665   (list new-pls (simplify-svar-ldata new-red-pls)))
 666
 667 (defun blowup-pls (pls open-for-eqns eqns  &key simplify-homogeneous
 668                    throw-out-components-in-exceptional-divisor &aux codim)
 669   (setq pls (copy-list-structure pls))
 670   (multiple-value-bind
 671     ( new-pls opens-not-to-blowup)
 672       (prepare-for-blowup pls open-for-eqns eqns)
 673     (show opens-not-to-blowup)
 674     (cond ((ldatap eqns)(setq codim (length (ldata-eqns eqns))))
 675           (t (setq codim (length eqns))))
 676
 677     (format t "~%****      Blowing a subscheme of codimension ~A  *** "codim )
 678     (des eqns)
 679     (blowup-sheaf new-pls codim :opens-not-to-blowup opens-not-to-blowup
 680                   :throw-out-components-in-exceptional-divisor
 681                   throw-out-components-in-exceptional-divisor
 682                   :simplify-homogeneous simplify-homogeneous)))
 683
 684 (defun blowup-pls (pls open-for-eqns eqns  &key simplify-homogeneous
 685                       add-exceptional-divisor-ldata &aux codim)
 686   (setq pls (copy-list-structure pls))
 687   (multiple-value-bind
 688     ( new-pls opens-not-to-blowup)
 689       (prepare-for-blowup pls open-for-eqns eqns)
 690     (show opens-not-to-blowup)
 691     (cond ((ldatap eqns)(setq codim (length (ldata-eqns eqns))))
 692           (t (setq codim (length eqns))))
 693     (format t "~%****      Blowing a subscheme of codimension ~A  *** "codim )
 694     (des eqns)
 695     (blowup-sheaf new-pls codim :opens-not-to-blowup opens-not-to-blowup
 696                    :add-exceptional-divisor-ldata  add-exceptional-divisor-ldata
 697                   :simplify-homogeneous simplify-homogeneous)))
 698
 699
 700 (defun blowup-keep-exceptional (pls  &key red-pls-list open-eqns-list &aux  open eqns
 701                                 ld-num op-num red-pls)
 702   (cond (red-pls-list (setq red-pls (first red-pls-list)
 703                             op-num (second red-pls-list)
 704                             ld-num (third red-pls-list ))
 705                       (setq open (nth op-num (pls-opens red-pls)))
 706                       (setq eqns (nth ld-num (nth op-num (pls-data red-pls)))))
 707         (t (setq open (first open-eqns-list))
 708            (setq eqns (second open-eqns-list))))
 709   (let ((*leave-unit-ldata* t))
 710     (blowup-pls pls open  eqns  :simplify-homogeneous t
 711                 :add-exceptional-divisor-ldata t)))
 712 ;(defun blowup-pls-by-component-of-reduced (pls open-number red-pls
 713 ;                                          red-open-number ld-number &aux eqns)
 714 ;  (setq eqns (ldata-eqns
 715 ;              (apply-rmap (find-ring-map
 716 ;                            (nth red-open-number (pls-opens red-pls))
 717 ;                            (nth open-number (pls-opens pls)))
 718 ;                          (nth ld-number (pls-ldata red-pls red-open-number :ldata-number
 719 ;                                                    ld-number)))))
 720 ;  (format t "~%Blowing up ") (des eqns) (format t "on open ~A" open-number)
 721 ;  (blowup-pls pls open-number eqns))
 722
 723 (defun blowup-pls-by-component-of-reduced (pls  red-pls
 724                                            red-open-number ld-number &aux open eqns)
 725   (setq eqns  (nth ld-number      (nth red-open-number (pls-data red-pls))))
 726   (setq open (nth red-open-number (pls-opens red-pls) ))
 727   (format t "~%Blowing up ") (des eqns) (format t " with coordinates ")
 728   (des open)
 729   (blowup-pls pls  open
 730              eqns))
 731
 732 ;checks if the ldata eqns are in triangular linear form
 733 ;that is if there's a linear variable in eqn1 not occurring in remaining eqns
 734 ;etc.
 735
 736
 737 (defun linear-ldatap (ldata  &key (open-g 1) &aux (eqns (ldata-eqns ldata))
 738                       var-occurs lin-vars tem)
 739   (block sue
 740 ;    (show eqns)
 741   (setq lin-vars (loop for v in eqns
 742         collecting (setq tem (all-linear-variables v open-g))
 743         when (null tem)do (return-from sue nil)))
 744 ; (show lin-vars)
 745   (setq var-occurs (mapcar 'list-variables eqns))
 746   (loop for v in lin-vars
 747         for w on var-occurs
 748         collecting
 749         (loop named kay for vv in v
 750               do
 751               ;;check vv doesn't occur in remaining equations
 752               ;;and if not put it in collection
 753               (loop for ww in (cdr w)
 754                     when (member vv ww :test #'eq) do (return nil)
 755                     finally (return-from kay vv))
 756               finally (return-from sue nil)))))
 757
 758 (defun reduce-linear-ldata (ldata variables &key (open-g) &aux f eqns vars)
 759   (setq eqns (reverse (ldata-eqns ldata)))
 760   (setq vars (reverse variables))
 761   (loop while fns
 762         for v in vars
 763         with fns = eqns
 764         with used-up
 765         do
 766         (setq f (car fns))
 767         (setq fns (replace-functions fns   f v :invertible-g open-g))
 768
 769         (push f used-up)
 770         finally (return used-up)))
 771
 772 ;
 773 ;(defun try-harder-delete-redundant-components (lis-dat &optional invertible-g)
 774 ;  (setq codims  (loop for v in lis-dat
 775 ;        collecting (length (ldata-eqns v))))
 776 ;  (loop for cod1 in codims
 777 ;       for ld1
 778 ;       when (and  (> max (apply 'max codims))
 779 ;                  (not (linear-ldatap ld1) :open-g invertible-g))
 780 ;       do
 781 ;       (setq simp-ld (ldata-simplifications ld1 invertible-g))
 782 ;       (loop for cod2 in codims
 783 ;             for ld2 in lis-dat
 784 ;             when ( > cod2 cod1)
 785 ;             do (grobner-subset '(1) (ldata-eqns ld2) invertible-g)
 786 ;             (loop for s-ld in simp-ld
 787 ;                   do
 788 ;                   (loop for fn in (ldata-eqns s-ld)
 789 ;                         when (null (any-linearp fn invertible-g))
 790 ;                         do (setq deg-1-vars (degree-one-variables fn))
 791 ;                         (loop for v in deg-1-vars
 792 ;                               do (check-containment ld1 ld2 invert))
 793 ;
 794
 795 (defmacro setq-eqns ( &rest alt-list)
 796   (loop for v on alt-list by #'cddr
 797         collecting `(cond ((ldatap ,(second v))(setf  ,(car v ) (ldata-eqns ,(second v))))
 798                           (t (setf ,(car v) ,(second v))))
 799                             into tem
 800         finally (return (cons 'progn tem))))
 801
 802 ;;tries to check that (var ld2) is subset (var ld1) by checking where dich-f is not zero
 803 ;;and where it is zero  .  If simplify is t it will try simplifying ld1 further
 804 ;;into components on the open where dich-f is not zero, and then try to verify
 805 ;;that (var ld2) is contained in one of the components so obtained.
 806 ;;this handles the example below where the first ideal has zero set contained
 807 ;;in the zero set of the second one but one still has to simplify the second one..
 808 ;ld2:[2*X1*X4*X6+X2*X4^2-X1^2*X3,X7,X8,X5,X6^2+X2*X3]
 809 ;ld1:[X5,X8,X4*X7-X4*X6+X1*X3,X7^2-X6^2-X2*X3]
 810 ;(SETQ LLD (RERAT (QUOTE ((LDATA ((X6 1 (X4 1 (X1 1 2)) 0 (X4 2 (X2 1 1) 0 (X3 1 (X1 2 -1)))) (X7 1 1) (X8 1 1) (X5 1 1) (X6 2 1 0 (X3 1 (X2 1 1)))) 1 0) (LDATA ((X5 1 1) (X8 1 1) (X7 1 (X4 1 1) 0 (X6 1 (X4 1 -1) 0 (X3 1 (X1 1 1)))) (X7 2 1 0 (X6 2 -1 0 (X3 1 (X2 1 -1))))) 1 0)))))
 811
 812
 813
 814
 815 (defun check-containment (  ld1 ld2 dich-f &key (open-g 1) simplify
 816                           &aux eqn1 eqn2   lis-ld)
 817   (block  sue
 818     (setq-eqns eqn1 ld1 eqn2 ld2)
 819 ;  (setq cont (grobner-subset eqn1 eqn2  (nplcm open-g dich-f)))
 820 ;  (setq cont2 (grobner-subset (cons dich-f eqn1 ) (cons dich-f eqn2) open-g))
 821 ; (show cont2)
 822
 823     (cond (simplify (setq lis-ld (ldata-simplifications ld1 :open-g (nplcm open-g dich-f))))
 824           (t (setq lis-ld (list ld1))))
 825     ;;run through the splitting of lis-ld into components
 826     (loop for ld in lis-ld
 827           when
 828           (check-containment1 ld eqn2 dich-f :open-g open-g)
 829           do (return-from sue t))))
 830
 831
 832 (defun check-containment1 (  ld1 ld2 dich-f &key (open-g 1)
 833                           &aux eqn1 eqn2 cont cont2 )
 834   (setq-eqns eqn1 ld1 eqn2 ld2)
 835   (setq cont (grobner-subset eqn1 eqn2  (nplcm open-g dich-f)))
 836   (cond (cont
 837          (setq cont2 (grobner-subset (cons dich-f eqn1 ) (cons dich-f eqn2) open-g))
 838          cont2))
 839   (and cont cont2))
 840
 841 (defmacro sim-hom (x)`(setq  ,(intern (string-append "SIMP-" (string x)))
 842           (simplify-svar-homogeneous ,x)))
 843
 844 (defmacro sim-ld (x)`(setq  ,(intern (string-append "RED-" (string-trim "SIMP-"
 845                                                                         (string x))))
 846           (simplify-svar-ldata ,x)))
 847
 848 (defun sub-scheme (pls op-number ld-number)
 849   (construct-pre-ldata-sheaves :opens (list (nth op-number
 850                                                                (pls-opens pls)))
 851                                      :data (list (list (nth ld-number
 852                                                             (nth op-number
 853                                                                  (pls-data pls)))))))
 854
 855 (defun open-sub-scheme (pls &rest op-numbers)
 856   (loop for i in op-numbers
 857         collecting (nth i (pls-opens pls)) into ops
 858         collecting (nth i(pls-data pls) ) into dat
 859         finally
 860         (return (construct-pre-ldata-sheaves :opens ops
 861                                      :data dat))))
 862
 863
 864 (defun intersect-ldata (ld1 ld2 &key ( open-g 1) &aux eqns *refine-opens*)
 865   (setq eqns (append (ldata-eqns ld1) (ldata-eqns ld2)))
 866   (ldata-simplifications (make-ldata :eqns eqns ) :open-g open-g))
 867
 868 (defmacro save-rational (name &optional (file-name "r20:ps:<atp.schelter>kill.tem" aa))
 869   (cond ((null aa)
 870          (setq file-name (merge-pathnames file-name))
 871          (setq file-name (send file-name :new-name (string name)))))
 872
 873   `(with-open-file  (st ',file-name :direction :output)
 874        (let (*print-gensym*)
 875      (format st ";;; -*- mode: lisp; package: cl-maxima; syntax: common-lisp -*-
 876      ~%(setq ~A (rerat '" ',name) (format st "~s" ,name)
 877      (format st "))")
 878      ',file-name)))
 879
 880 (defun save-parts  ( obj file-name &aux  me)
 881   (with-open-file (st  file-name :direction :output)
 882     (let (*print-gensym*)
 883     (format st ";;; -*- mode: lisp; package: cl-maxima; syntax: common-lisp; Mode: LISP -*- ~%")
 884     (loop for v in obj
 885           for i from 1
 886           do
 887           (setq me `(setq ,(intern (format nil "PART~A" i)) (rerat ',v)))
 888           (format st "~%;;next part ~%~S" me))
 889     (send st :pathname))))
 890
 891 (defun save-blew-up ( &aux  me)
 892   (with-open-file (st "r20:ps:<atp.schelter>blew-up.tem" :direction :output)
 893     (let (*print-gensym*)
 894     (format st ";;; -*- mode: lisp; package: cl-maxima; syntax: common-lisp; Mode: LISP -*- ~%")
 895     (loop for v in *blew-up*
 896           for i from 1
 897           do
 898           (setq me            `(setq ,(intern (format nil "PART~A" i)) (rerat ',v)))
 899           (format st "~%;;next part ~%~S" me))
 900     (send st :pathname))))
 901
 902
 903
 904
 905 ;(defun triangular (fns &key (open-g 1) variables-solved-for variables-in-gm-coefficients)
 906 ;  (cond ((null fns) 'ok)
 907 ;       ((eql (length fns ) 1)
 908 ;        (cond ((loop for var in (degree-one-variables (car fns))
 909 ;                     when (may-invertp (pcoeff (car fns) (list var 1 1)) open-g)
 910 ;                     do (return var)))
 911 ;              (t (setq tem (nsubst nil 'ok (gm-all-prepared (car fns) open-g)))
 912 ;                 (loop named sue
 913 ;                   for possible in tem
 914 ;                       do
 915 ;                       (loop for var in possible
 916 ;                             when (not (member var variables-in-gm-coefficients :test #'eq))
 917 ;                             do (return-from sue 'ok-gm))))))
 918 ;       (t (loop for f in fns
 919 ;                    collecting (setq var-1 (degree-one-variables f)) into vars-1
 920 ;                    do
 921 ;                    (loop for va in var-1
 922 ;                          when (and (poly-linearp f va open-g)
 923 ;                                    (setq tem (triangular (delete f (copy-list fns))
 924 ;                                                :variables-solved-for                                                                    variables-solved-for
 925 ;                                                :variables-in-gm-coefficients
 926 ;                                                variables-in-gm-coefficients)))
 927 ;                          do (return (cons va tem)))
 928
 929
 930 ;;;will guarantee that the system has the form
 931 ;;  f1 (u1,   us,x1)
 932 ;;  f2 (u1,   us,x1,x2)
 933 ;;  ..
 934 ;;  fn (u1,   us,x1,x2,   xn)
 935 ;; and  xi is linear in fi and none of the variables in the  coefficient of xi in fi
 936 ;; occur among the x1,x2,x3,.. xi-1.  This will guarantee irreducibility, since if
 937 ;; we did xn>xn-1>..>x1>us>...>u1 then
 938 ;; the system is irreducible.  We show fi is irreducible modulo the f1,.. fi-1.
 939 ;; To see this: If it could be factored it would mean the coefficient of xi, say
 940 ;; ci and the constant di had a common factor.   The coefficient ci is reduced
 941 ;; withrespect to the preceding system, and so
 942
 943 (defun delete-nth ( n a-list)
 944   (loop for i from 0
 945         for v in a-list
 946         when
 947         (not (eql n i))
 948         collecting v))
 949 ;;will return a list of variables ANSWER occurring with degree 1 in various equations
 950 ;;such that the equations may be so reordered such that the
 951 ;;ith variable in Answer
 952 ;; u1+x1,
 953 ;; x1+u1*x2
 954 ;; Is such a system which is not irreducible.
 955
 956 (defun linear-triangularp (fns &key variables-solved-for variables-in-prev-eqns
 957                           varlist-of-fns &aux tem )
 958   (cond ((null varlist-of-fns)
 959          (setq varlist-of-fns (loop for f in fns collecting (list-variables f)))))
 960   (cond ((null fns)variables-solved-for)
 961         (t
 962          (loop named sue
 963                for f in fns
 964                for i from 0
 965                for varl in varlist-of-fns
 966                do (loop for va in varl
 967                         when (and
 968                                (not (member va variables-solved-for :test #'eq))
 969                                (not (member va variables-in-prev-eqns :test #'eq))
 970                                (eql (pdegree f va) 1)
 971                                (setq tem  (linear-triangularp
 972                                                  (delete-nth i fns)
 973                                                  :varlist-of-fns
 974                                                  (delete-nth i varlist-of-fns)
 975                                                  :variables-solved-for
 976                                                  (cons va variables-solved-for)
 977                                                  :variables-in-prev-eqns
 978                                                  (append varl variables-in-prev-eqns))))
 979                         do (return-from sue tem ))))))
 980
 981 ;;returns a list of variables each of which occurs with degree 1 in exactly
 982 ;;one function and does not occur in any other fn.
 983 (defun linear-solvedp (fns &key variables-solved-for variables-in-prev-eqns
 984                           order-functions
 985                           varlist-of-fns &aux tem )
 986   (cond ((null varlist-of-fns)
 987          (setq varlist-of-fns (loop for f in fns collecting (list-variables f)))))
 988   (cond ((null fns)variables-solved-for)
 989         (t
 990          (loop named sue
 991                for f in fns
 992                for i from 0
 993                for varl in varlist-of-fns
 994                do (loop for va in varl
 995                         when (and
 996                                (not (member va variables-solved-for :test #'eq))
 997                                (not (member va variables-in-prev-eqns :test #'eq))
 998                                (not (loop for lis in varlist-of-fns
 999                                           for j from 0
1000                                           when (and  (not (eql i j))
1001                                                      (member va lis :test #'eq))
1002                                           do (return t)))
1003                                (eql (pdegree f va) 1)
1004                                (setq tem
1005                                  (linear-solvedp
1006                                                  (delete-nth i fns)
1007                                                  :order-functions order-functions
1008                                                  :varlist-of-fns
1009                                                  (delete-nth i varlist-of-fns)
1010                                                  :variables-solved-for
1011                                                  (cons va variables-solved-for)
1012                                                  :variables-in-prev-eqns
1013                                                  (append varl variables-in-prev-eqns))))
1014                         do (return-from sue tem))))))
1015
1016 (defun linear-solvedp (fns &key variables-solved-for variables-in-prev-eqns order-functions varlist-of-fns &aux tem ord-fns)
1017   (cond ((null varlist-of-fns)
1018          (setq varlist-of-fns (loop for f in fns collecting (list-variables f)))))
1019   (cond ((null fns)(values variables-solved-for (subst nil t order-functions)))
1020         (t
1021          (loop named sue
1022                 for f in fns
1023                 for i from 0
1024                 for varl in varlist-of-fns
1025                 do (loop for va in varl
1026                           when (and
1027                                 (not (member va variables-solved-for :test #'eq))
1028                                 (not (member va variables-in-prev-eqns :test #'eq))
1029                                 (not (loop for lis in varlist-of-fns
1030                                             for j from 0
1031                                             when (and  (not (eql i j))
1032                                                        (member va lis :test #'eq))
1033                                             do (return t)))
1034                                 (eql (pdegree f va) 1)
1035                                 (progn (multiple-value-setq
1036                                            (tem ord-fns)
1037                                          (linear-solvedp
1038                                           (delete-nth i fns)
1039                                           :varlist-of-fns
1040                                           (delete-nth i varlist-of-fns)
1041                                           :variables-solved-for
1042                                           (cons va variables-solved-for)
1043                                           :order-functions
1044                                           (cond (order-functions
1045                                                  (cons f order-functions)))
1046                                           :variables-in-prev-eqns
1047                                           (append varl variables-in-prev-eqns)))))
1048                           do (return-from sue (values  tem ord-fns)))))))
1049
1050
1051 (defun linear-solvedp (fns &rest ignore &aux bad good all-degs varl)
1052   (setq varl (list-variables fns))
1053   (let ((genvar (nreverse (sort varl 'pointergp))))
1054     (setq all-degs (loop for fn in fns collecting (pdegreevector fn)))
1055     (loop named kay
1056        for degs in all-degs
1057        for fn-number from 0
1058        do (loop named sue for i in degs
1059              for j from 0
1060              when (and (eql i 1) (not (member j bad :test #'eq)))
1061              do
1062                (loop for vv in all-degs
1063                   when (not (eql degs vv))
1064                   do (cond ((not (zerop (nth j vv) ))
1065                             (push j bad) (return 'done)))
1066                   finally (push j good) (return-from sue))
1067              finally (return-from kay (setq good nil))))
1068
1069     (loop for j in (nreverse good)
1070        collecting (nth j genvar))))
1071
1072 ;(defun plain-simplify-svar-ldata (pls &aux *refine-opens* *stop-simplify* tem)
1073 ;  (loop for op in (pls-opens pls)
1074 ;       for lis in (pls-data pls)
1075 ;       do
1076 ;       (loop for ld in lis
1077 ;             appending (ldata-simplifications ld :open-g (zopen-inequality op))
1078 ;             into all-ld
1079 ;             finally
1080 ;             (setq tem (delete-redundant-ldata all-ld :gg (zopen-inequality op)
1081 ;                                             :ignore-ldata-inequalities t)))
1082 ;       collecting tem into simp-ld
1083 ;       finally
1084 ;
1085 ;       (return
1086 ;         (construct-pre-ldata-sheaves :data simp-ld :opens (pls-opens pls)))))
1087
1088
1089
1090 (defun plain-simplify-svar-ldata (pls &aux   list-of-dat simp-one-pls one-pls
1091                                   ineq)
1092   (loop for op in (pls-opens pls)
1093         for i from 0
1094         for lis in (pls-data pls)
1095         do
1096         (setq one-pls (open-sub-scheme pls i))
1097         (setq simp-one-pls (simplify-svar-ldata one-pls))
1098         (cond ((eql (length (pls-opens simp-one-pls)) 1)
1099                (setq list-of-dat (car  (pls-data simp-one-pls))))
1100               (t (loop for ope in (pls-opens simp-one-pls)
1101                        for lis-dat in (pls-data simp-one-pls)
1102                        do (setq ineq (num (ratreduce (zopen-inequality ope)
1103                                                      (zopen-inequality op))))
1104                        appending
1105                        (loop for ld in lis-dat
1106                              collecting (make-ldata :eqns
1107                                                     (contract-ideal-localization
1108                                                       (ldata-eqns ld)
1109                                                       ineq)))
1110                        into all-dat
1111                        finally
1112                        (setq list-of-dat
1113                              (delete-redundant-ldata  all-dat :gg (zopen-inequality
1114                                                                           op))))))
1115         collecting list-of-dat into all-data
1116         finally (return
1117                   (construct-pre-ldata-sheaves
1118                     :opens (pls-opens pls)
1119                     :data all-data))))
1120
1121
1122 (defun svar-simp (ld &key (open-g 1))
1123   (simplify-svar-ldata (construct-pre-ldata-sheaves :opens (list (make-normal-zopen nil 8 open-g)) :data (list (list ld)))))
1124
1125 (defun jacobian (eqns &key ( variables (list-variables eqns))
1126                  ( codim (length eqns)) &aux det mat row-tables col-tables)
1127   (setq mat (matrix-rows (jacobian-matrix eqns variables)))
1128 ;  (displa (cons '($matrix) (loop for v in mat collecting (cons'(mlist)
1129 ;                                                            (mapcar 'new-disrep v)))))
1130 ;  (show variables)
1131   (setq row-tables (list-tableaux codim (length eqns)))
1132   (setq col-tables (list-tableaux  codim (length variables)))
1133   (loop named sue
1134         for ro-tabl in row-tables
1135         appending
1136         (loop for tabl in col-tables
1137               do (setq *sparse-matrix* (convert-to-sparse-matrix
1138                                          mat :columns-to-use tabl
1139                                          :rows-to-use ro-tabl
1140
1141                                          :re-use-sparse-matrix  *sparse-matrix*))
1142               (setq det (sp-determinant *sparse-matrix*))
1143               when (numberp det)
1144               do      (cond ((not ($zerop det)) (sp-show-matrix *sparse-matrix*)
1145                              (show tabl ro-tabl)
1146                              (return-from sue '(1))))
1147               else collecting det into jacob
1148               finally (return jacob))))
1149
1150
1151 (defun rank-generic-reduce-jacobian (eqns &optional (variables(list-variables eqns)) &aux mat det )
1152 ;  (declare (values must-invert rank))
1153  (setq mat (matrix-rows (jacobian-matrix eqns variables)))
1154  (setq *sparse-matrix* (convert-to-sparse-matrix mat :re-use-sparse-matrix *sparse-matrix*))
1155  (setq det (sp-determinant *sparse-matrix*))
1156  (values det (sp-number-of-pivots *sparse-matrix*)))
1157 (defun reduce-jacobian (eqns &optional (variables (list-variables eqns)) &aux mat det)
1158 ;  (declare (values must-invert rank))
1159  (setq mat (matrix-rows (jacobian-matrix eqns variables)))
1160  (setq *sparse-matrix* (convert-to-sparse-matrix mat :re-use-sparse-matrix *sparse-matrix*))
1161  (setq det (sp-determinant *sparse-matrix*))
1162  (values det (sp-number-of-pivots *sparse-matrix*)))
1163 ;;(non-singular-complete-intersection-p  (st-rat #$[x^2+w+y^3,u+y^2+x^2,w^2+y^2+y]$) )
1164 ;;the above has a 4 singular points given by the following equations:
1165 ;;the jacobian together with the equations is a mess to try and solve
1166 ;[2*Y+1,
1167 ;4*X^2+4*U+1,
1168 ;8*W-8*U-3,
1169 ;64*U^2+48*U-7]
1170 (defun non-singular-complete-intersection-p (eqns &key (open-g 1)
1171                                              (use-simplification t)
1172                                              &aux new
1173                                              singular-locus eqns-and-jacob)
1174 ;  (declare (values ( singular-locus non-singularp)))
1175   (cond ((ldatap eqns) (setq eqns (ldata-eqns eqns))))
1176   (setq eqns-and-jacob  (mapcar 'square-free (append eqns (jacobian eqns))))
1177   (cond (use-simplification
1178          (setq new (ldata-simplifications
1179                      (make-ldata :eqns eqns-and-jacob) :open-g open-g))
1180          (cond ((null new) (values '(1)  t))
1181                (t  (values   new nil))))
1182         (t
1183          (cond ((unit-idealp eqns-and-jacob open-g)
1184                 (values   '(1) t))
1185                (t (setq singular-locus (poly-relations-from-simplifications ))
1186                   (values (list (make-ldata :eqns singular-locus)) nil))))))
1187
1188 (defun non-singular-complete-intersection-p (eqns &key (open-g 1) codim
1189                                              (use-simplification t)
1190                                              (variables (list-variables eqns))
1191                                              &aux new
1192                                              singular-locus eqns-and-jacob)
1193 ;  (declare (values ( singular-locus non-singularp)))
1194   (cond ((ldatap eqns) (setq eqns (ldata-eqns eqns))))
1195   (cond ((null codim) (setq codim (length eqns))))
1196   (setq eqns-and-jacob  (mapcar 'square-free (append eqns (jacobian eqns :variables variables
1197                                                                     :codim codim))))
1198   (cond (use-simplification
1199          (setq new (ldata-simplifications
1200                      (make-ldata :eqns eqns-and-jacob) :open-g open-g))
1201          (cond ((null new) (values '(1)  t))
1202                (t  (values   new nil))))
1203         (t
1204          (cond ((unit-idealp eqns-and-jacob open-g)
1205                 (values   '(1) t))
1206                (t (setq singular-locus (poly-relations-from-simplifications ))
1207                   (values (list (make-ldata :eqns singular-locus)) nil))))))
1208 (defun simplify-dichotomy (ldata fn &key (open-g 1) &aux answ)
1209   (setq answ (append (ldata-simplifications (make-ldata :eqns (cons fn (ldata-eqns ldata) ))
1210                                             :open-g open-g)
1211                      (loop for v in
1212                           (ldata-simplifications ldata :open-g (sftimes open-g fn))
1213                         collecting
1214                           (make-ldata :eqns (contract-ideal-localization (ldata-eqns v) fn)))))
1215   (delete-redundant-ldata answ :open-g open-g :ignore-ldata-inequalities t))
1216
1217 (defun non-trivial-open-sub-scheme (pls &aux non-trivial)
1218   (setq non-trivial (loop for lis in (pls-data pls)
1219                           for i from 0
1220                           when (and lis (not (member 1 (ldata-eqns (car lis)) :test #'equal)))
1221                           collecting i))
1222   (show non-trivial) (break t)
1223   (apply #'open-sub-scheme pls non-trivial))
1224
1225 (defun show-divisors (pls &aux tem)
1226   (loop for lis in  (pls-data pls)
1227         for i from 0
1228         collecting (setq tem (loop for ld in (cdr lis)
1229                          when (not (numberp (setq tem (car (ldata-eqns ld)))))
1230                          collecting tem))
1231         do (format t "~%On open ~D the divisors are ~/maxima::tilde-q-fsh/" i tem)))
1232
1233 (defmacro alter-ldata (ldat &rest keys)
1234   (loop for (key repl)  on keys by #'cddr
1235         with u = '.ldat.
1236         collecting `(setf (,(intern (format nil "LDATA-~A" key)) .ldat.) ,repl) into body
1237         finally (cond (body (return `(let ((.ldat. ,ldat)) ,@body .ldat.)))
1238                       (t (return `(progn ,ldat))))))
1239
1240 (defun jacobian-dichotomy (ldata &key (open-g 1) &aux ld1 ld2 gg answer)
1241   (multiple-value-bind (c-reqd rank) (rank-generic-reduce-jacobian (ldata-eqns ldata))
1242     (cond ((may-invertp c-reqd (setq gg (sftimes open-g (ldata-inequality ldata))))
1243            (cond ((eql (length (ldata-eqns ldata)) rank) (format t "~%Non singular ldata"))
1244                  (t (format t "~%singular ldata")))
1245            (setq answer (list ldata)))
1246           ((eql rank (length (ldata-eqns ldata)))
1247            (setq ld1 (copy-list-structure ldata))
1248            (setq ld2 (copy-list ldata))
1249            (alter-ldata ld1 eqns (append (ldata-eqns ldata) (list c-reqd))
1250                         inequality gg)
1251            (alter-ldata ld2 inequality (sftimes gg c-reqd))
1252            (format t "~%Breaking into dichotomy on:")
1253            (sh c-reqd)
1254            (format t "%original ldata followed by consequents:" )
1255            (des ldata)(des ld1) (des ld2 )
1256            (setq answer
1257                  (append (ldata-simplifications ld1 :open-g open-g :recursive-p t)
1258                          (ldata-simplifications ld2 :open-g open-g :recursive-p t))))
1259           (t (setq answer (list ldata))))
1260     (delete-redundant-ldata answer :open-g open-g :ignore-ldata-inequalities t)))
1261
1262 (defun $list_components (lis-ldat)
1263   (cond ((eq (car lis-ldat) 'pre-ldata-sheaves)
1264          (setq lis-ldat (apply 'append (pls-data lis-ldat))))
1265         ((ldatap lis-ldat) (setq lis-ldat (list lis-ldat))))
1266   (cons '(mlist) (loop for v in lis-ldat
1267                     collecting (cons '(mlist) (mapcar 'new-disrep (ldata-eqns v))))))
1268
1269
1270 (defun $simple_solve (eqns &key (open-g 1) &aux ld vari solu  )
1271   (setq ld (make-ldata :eqns (st-rat eqns)))
1272   (multiple-value-bind (to-invert rank) (rank-generic-reduce-jacobian (ldata-eqns ld))
1273     (setq to-invert (sftimes 1 to-invert))
1274     (cond((eql  (length (ldata-eqns ld)) rank)
1275           (cond ((may-invertp to-invert open-g))
1276                 (t(format t "having to invert ~/maxima::tilde-q-fsh/" to-invert)
1277                   (setq open-g (sftimes open-g to-invert))))))
1278     (setq vari (linear-ldatap ld :open-g (setq open-g(st-rat open-g))))
1279     (cond ((null vari) (fsignal "equations not linear-ldatap")))
1280     (setq solu (solve-simple-system (ldata-eqns ld) vari :invertible open-g))
1281     (cons '(mlist) (loop for v in vari
1282                       for w in solu
1283                       collecting (list '(mequal) (get v 'disrep)  (new-disrep (st-rat  w)))))))