share/fourier_elim/fourier_elim.lisp

   1 #| Copyright 2007, 2008, 2021 by Barton Willis
   2
   3   This is free software; you can redistribute it and/or
   4   modify it under the terms of the GNU General Public License,
   5   http://www.gnu.org/copyleft/gpl.html.
   6
   7  This software has NO WARRANTY, not even the implied warranty of
   8  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9 |#
  10
  11
  12 (if (not ($get '$to_poly '$version)) ($load '$to_poly))
  13
  14 (mfuncall '$declare '$one_to_one '$feature)
  15 (mfuncall '$declare '$sinh '$one_to_one)
  16 (mfuncall '$declare '$log  '$one_to_one)
  17 (mfuncall '$declare '$tanh '$one_to_one)
  18 (mfuncall '$declare '$log '$increasing)
  19
  20 ;; The macro opcons is defined in src; opapply is defined in to_poly. So we
  21 ;; no longer define them here.
  22
  23 ;; Maybe I should use csign instead of these functions...
  24 (defun number-sign (x)
  25   (cond ((or (integerp x) (floatp x))
  26          (cond ((< x 0) '$neg)
  27                ((= x 0) '$zero)
  28                ((> x 0) '$pos)))
  29
  30         (($ratnump x) (number-sign ($num x)))
  31         (($bfloatp x) (number-sign (second x)))
  32         ((eq '$minf x) '$neg)
  33         ((member x '($inf $%pi $%e $%phi) :test #'eq) '$pos)
  34         ((member x '($%i '$infinity) :test #'eq) '$complex)
  35         (($constantp x) (constant-expression-sign x))
  36         (t nil)))
  37
  38 (defun constant-expression-sign (x)
  39   (let ((f1) (f2) (f3))
  40
  41     ;; I can do better. A big float running error evaluator would be nice.
  42     ;; (or interval arithmetic).  But for now let's evaluate with 25, 50,
  43     ;; and 75 digits. Provided all three agree, and |f3| > 2^-50, return the common sign.
  44
  45     ;; I need to check for complex...
  46
  47     (setq f1 (bind-fpprec 25 ($bfloat x)))
  48     (setq f2 (bind-fpprec 50 ($bfloat x)))
  49     (setq f3 (bind-fpprec 75 ($bfloat x)))
  50
  51     (if (and ($bfloatp f1) ($bfloatp f2) ($bfloatp f3)
  52              (eq (number-sign f1) (number-sign f2))
  53              (eq (number-sign f2) (number-sign f3))
  54              (> (first (last f3)) -50)) (number-sign f1) nil)))
  55
  56 ;; This isn't exactly the same as CARTESIAN-PRODUCT in src/nset.lisp,
  57 ;; so give it a different name to avoid name collision.
  58 (defun fourier_elim-cartesian-product (b)
  59   (cond ((null b)
  60          nil)
  61         (t
  62          (let ((a) (acc (mapcar #'list (car b))))
  63            (setq b (cdr b))
  64            (dolist (bi b acc)
  65              (setq a nil)
  66              (dolist (bij bi (setq acc a))
  67                (setq a (append a (mapcar #'(lambda (x) `(,@x ,bij)) acc)))))))))
  68
  69 (defun expand-and-over-or (e)
  70   (cond (($mapatom e) e)
  71         ((op-equalp e 'mand)
  72          (setq e (mapcar #'(lambda (s) (if (op-equalp s 'mor) (margs s) (list s))) (margs e)))
  73          (setq e (fourier_elim-cartesian-product e))
  74          (setq e (mapcar #'(lambda (s) (opapply 'mand s)) e))
  75          (opapply 'mor e))
  76         (t (opapply (mop e) (mapcar 'expand-and-over-or (margs e))))))
  77
  78 (defun make-positive-product (l)
  79   (let* ((n (length l)) (m (expt 2 n)) (k) (p) (acc nil))
  80     (dotimes (i m)
  81       (if (evenp (logcount i))
  82           (progn
  83             (setq p nil)
  84             (setq k 0)
  85             (dolist (li l)
  86               (push (apply 'fe> (if (logbitp k i) (list 0 li nil) (list li 0 nil))) p)
  87               (incf k))
  88
  89             (setq p (opapply 'mand p))
  90             (push p acc))))
  91     (opapply 'mor acc)))
  92
  93 (defun in-real-domain (e)
  94   (cond ((eq e '$%i) nil)
  95         (($mapatom e) t)
  96
  97         ((op-equalp e 'mplus 'mtimes '%atan '%cos '%sin '%cosh '%sinh '%asinh)
  98          (opapply 'mand (mapcar 'in-real-domain (margs e))))
  99
 100         ((op-equalp e '%acos '%asin)
 101          (setq e (first (margs e)))
 102          (opapply 'mand (list (fe>= e -1) (fe>= 1 e))))
 103
 104         ((op-equalp e '%log)
 105          (fe> (first (margs e)) 0))
 106
 107         ((op-equalp e '%acosh)
 108          (fe>= (first (margs e)) 1))
 109
 110         ((op-equalp e 'mexpt)
 111          (let ((x (first (margs e))) (n (second (margs e))))
 112            (cond ((member ($compare n 0) '("<" "<=") :test 'equal)
 113                   (take '(mand) (m-neq x 0) (in-real-domain x) (in-real-domain n)))
 114                  (t
 115                   (take '(mand) (m-neq x 0) (in-real-domain x) (in-real-domain n))))))
 116         (t (m= e (take '($conjugate) e)))))
 117
 118 (defun freeof-floats (e)
 119   (if ($mapatom e) (not (or (floatp e) ($bfloatp e)))
 120     (every 'freeof-floats (margs (ratdisrep e)))))
 121
 122 ;; Splitify an expression; this does e -> ((e1, boolean) (e2, boolean), ...). Examples:
 123
 124 ;;   x -->  ((x t)),
 125 ;;   |x| --> ((-x, x < 0) (x, x >= 0)),
 126 ;;   max(a,b) --> ((a, a > b) (b, b >= a)),
 127 ;;   a + |b| --> a + ((-b, b < 0) (b, b >=0)) --> ((a-b, b < 0) (a + b, b >= 0)),
 128
 129 (defun $splitify (e)
 130   (cons '(mlist) (mapcar #'(lambda (s) (push '(mlist) s)) (splitify e))))
 131
 132 (defun splitify (e)
 133   (let ((acc nil) (f))
 134     (cond (($mapatom e) (list (list e t)))
 135
 136           ((op-equalp e '$max)
 137            (splitify (max-to-abs (margs e))))
 138
 139           ((op-equalp e '$min)
 140            (splitify (min-to-abs (margs e))))
 141
 142           ((op-equalp e 'mabs)
 143            (setq e (first (margs e)))
 144            (list (list (neg e) (fe> 0 e)) (list e (fe>= e 0))))
 145
 146           ((op-equalp e 'mexpt)
 147            (mapcar #'(lambda (s) (list (take '(mexpt) (first s) (third e)) (second s))) (splitify (second e))))
 148
 149           ((op-equalp e 'mplus 'mtimes)
 150            (setq f (if (op-equalp e 'mplus) 'add 'mult))
 151            (setq e (mapcar 'splitify (margs e)))
 152            (setq e (fourier_elim-cartesian-product e))
 153            (dolist (ek e acc)
 154              (push
 155               (reduce
 156                #'(lambda (a b)
 157                    (list (funcall f (first a) (first b)) (take '(mand) (second a) (second b)))) ek) acc)))
 158
 159           (t (list (list e t))))))
 160
 161 (defun fe> (a b &optional (expand nil) (use-splitify t))
 162   (let* ((z) (sgn) (z-split) (acc nil))
 163     (setq a ($ratdisrep a))
 164     (setq b ($ratdisrep b))
 165
 166     ;; This uses p / q > 0 == q^2 * (p / q) > 0 == p * q > 0. Skip the
 167     ;; p / q > 0 --> p * q > 0 transformation when z contains a floating
 168     ;; point number.
 169
 170     (setq z (sub a b))
 171     (setq z (if (freeof-floats z) ($factor (mul ($num z) ($denom z))) z))
 172     (setq sgn (compare-using-empty-context a b))
 173
 174     (cond ((equal sgn ">") t)
 175
 176           ((member sgn '("<" "<=" "=") :test 'equal) nil)
 177
 178           ;; Catch four easy cases before we splitify z. Without the checks for the easy cases,
 179           ;; we'd get things like |x| > 1 --> (-1 < x & x < 0) or (x = 0) or (x < 1, x > 0).
 180
 181           ;; First, |a| > b --> (b < 0) or (a > b, b >= 0) or (-a > b, b >= 0).
 182
 183           ((and (op-equalp a 'mabs) ($freeof 'mabs '$min '$max b))
 184            (setq a (first (margs a)))
 185            (opapply 'mor (list
 186                           (opapply 'mand (list (fe> 0 b)))
 187                           (opapply 'mand (list (fe>= b 0) (fe> a b)))
 188                           (opapply 'mand (list (fe>= b 0) (fe> (neg a) b))))))
 189
 190           ;; Second, a > |b| == a > b and  a > -b.
 191
 192           ((and (op-equalp b 'mabs) ($freeof 'mabs '$min '$max a))
 193            (setq b (first (margs b)))
 194            (opapply 'mand (list (fe> a b) (fe> a (neg b)))))
 195
 196           ;; Third, min(a1,a2,..., an) > b --> a1 > b and a2 > b and .. an > b.
 197
 198           ((and (op-equalp a '$min) ($freeof 'mabs '$min '$max b))
 199            (opapply 'mand (mapcar #'(lambda (s) (fe> s b)) (margs a))))
 200
 201           ;; Fourth, a > max(b1,b2, ..., bn) --> a > b1 and a > b2 and ... and a > bn.
 202
 203           ((and (op-equalp b '$max) ($freeof 'mabs '$min '$max a))
 204            (opapply 'mand (mapcar #'(lambda (s) (fe> a s)) (margs b))))
 205
 206           ;; Do z^n > 0 --> z # 0  n even,  z > 0, n odd.
 207           ((and (op-equalp z 'mexpt) (integerp (third z)))
 208            (if (even (third z)) (m-neq (second z) 0) (fe> (second z) 0)))
 209
 210           ;; Do f(a) > f(b), where f is increasing --> a > b.
 211           ((and (not ($mapatom a)) (not ($mapatom b)) (eq (mop a) (mop b)) ($featurep (mop a) '$increasing))
 212            (opapply 'mand (list
 213                            (fe> (first (margs a)) (first (margs b)))
 214                            (in-real-domain a)
 215                            (in-real-domain b))))
 216
 217           ;; Do f(a) > f(b), where f is decreasing --> a < b.
 218           ((and (not ($mapatom a)) (not ($mapatom b)) (eq (mop a) (mop b)) ($featurep (mop a) '$decreasing))
 219            (fe> (first (margs b)) (first (margs a))))
 220
 221           ;; Do a^x > a^y, where a > 1 --> x > y,
 222           ((and (not ($mapatom a)) (not ($mapatom b)) (eq (mop a) (mop b))
 223                 (op-equalp a 'mexpt) (eql (second a) (second b))
 224                 (equal (compare-using-empty-context (second a) 1) ">"))
 225            (fe> (third a) (third b)))
 226
 227           ;; Do a * b > 0 --> (a > 0, b > 0) or (a < 0, b < 0). We only do this when
 228           ;; z has two or more non-constant factors. This check seems spendy--is there
 229           ;; a way to bailout before we get here?
 230
 231           ((and (op-equalp z 'mtimes) expand)
 232            (make-positive-product (margs z)))
 233           ;; Finally, take care of the abs, max, and min cases that the previous
 234           ;; four cases miss.
 235
 236           ((and use-splitify (op-equalp z '$max '$min 'mabs 'mtimes 'mplus))
 237            (setq z-split (splitify z))
 238            (dolist (zk z-split)
 239              (push `((mand) ,(fe> (first zk) 0 expand nil) ,@(rest zk)) acc))
 240            (push '(mor) acc)
 241            (simplifya acc nil))
 242
 243           (t
 244            (opapply 'mgreaterp (list z 0))))))
 245
 246 (defun m= (a b &optional (use-splitify t))
 247   (let* ((z (sub a b)) (nz) (acc) (z-split) (sgn  (compare-using-empty-context a b)))
 248     (setq z (if (freeof-floats z) ($factor z) z))
 249
 250     (cond
 251
 252      ;; If compare says a = b, return true.
 253      ((equal sgn "=") t)
 254
 255      ;; If compare says a # b, return false.
 256      ((member sgn '("<" ">" "#") :test 'equal) nil)
 257
 258      ;; for complex numbers, look at the real and imaginary parts.
 259      ((and (complex-number-p a '$numberp) (complex-number-p b '$numberp))
 260       (take '(mand) (m= ($realpart a) ($realpart b)) (m= ($imagpart a) ($imagpart b))))
 261
 262      ;; z^n = 0 --> false if n <= 0 else z = 0.
 263      ((and (op-equalp z 'mexpt) (mnump (third z)))
 264       (if (member (number-sign (third z)) '($neg $zero) :test #'eq) nil (m= (second z) 0)))
 265
 266      ;; f(a) = f(b), where f is one-to-one --> a = b.
 267      ((and (not ($mapatom a)) (not ($mapatom b)) (eq (mop a) (mop b)) ($featurep (mop a) '$one_to_one))
 268       (opapply 'mand (append (list (in-real-domain a) (in-real-domain b)) (mapcar #'m= (margs a) (margs b)))))
 269
 270      ;; a * b  = 0 --> a = 0 or b = 0; also a / b --> a = 0 and b # 0.
 271      ((op-equalp z 'mtimes)
 272       (setq nz (m-neq ($ratdenom z) 0))
 273       (expand-and-over-or
 274        (take '(mand) nz (opapply 'mor (mapcar #'(lambda (s) (m= s 0)) (margs z))))))
 275
 276      ((and use-splitify (op-equalp z '$max '$min 'mabs 'mtimes 'mplus))
 277       (setq z-split (splitify z))
 278       (dolist (zk z-split)
 279         (push `((mand) ,(m= (first zk) 0 nil) ,@(rest zk)) acc))
 280       (push '(mor) acc)
 281       (expand-and-over-or (simplifya acc nil)))
 282
 283      (t (take '(mequal) z 0)))))
 284
 285 (defun m-neq (a b)
 286   (let ((save-context $context) (new-context (gensym)) (sgn))
 287     (unwind-protect
 288      (progn
 289        (setq sgn (mnqp a b))
 290        (cond ((or (eq sgn t) (eq sgn nil)) sgn)
 291              ((eq $domain '$real) (opcons 'mor (fe> a b t) (fe> b a t)))
 292              (t (take '(mnot) (m= a b)))))
 293      (if ($member new-context $contexts) ($killcontext new-context))
 294      (setq $context save-context))))
 295
 296 (defun fe>= (a b)
 297   (let ((sgn (compare-using-empty-context a b)))
 298     (cond ((member sgn '(">=" ">") :test 'equal) t)
 299           ((equal sgn "<") nil)
 300           (t (opcons 'mor (fe> a b  t) (m= a b))))))
 301
 302 (defun standardize-inequality (e)
 303   (let ((a) (b))
 304     (cond ((op-equalp e 'mand)
 305            (opapply 'mand (mapcar 'standardize-inequality (margs e))))
 306
 307           ((op-equalp e 'mor)
 308            (opapply 'mor (mapcar 'standardize-inequality (margs e))))
 309
 310           ((or (mrelationp e) (op-equalp e '$equal))
 311            (setq a (second e))
 312            (setq b (third e))
 313
 314            (cond ((op-equalp e 'mlessp) (fe> b a t))
 315                  ((op-equalp e 'mleqp) (opcons 'mor (fe> b a  t) (m= a b )))
 316                  ((op-equalp e 'mequal '$equal) (m= a b ))
 317                  ((op-equalp e 'mgreaterp) (fe> a b  t))
 318                  ((op-equalp e 'mgeqp) (fe>= a b))
 319                  ((op-equalp e 'mnotequal) (m-neq a b))
 320                  (t e)))
 321           (t e))))
 322
 323 ;; affine_p(p,vars) returns true iff p is an affine polynomial in vars, that is,
 324 ;; that it is a polynomial in vars (a list) whose total degree in vars is no greater than 1.
 325
 326 ;; Stavros Macrakis wrote $affine_p and the tests (see rtest_fourier_elim.mac) for this function.
 327
 328 (defun $affine_p (p vl)
 329   (setq vl (require-list vl "affine_p"))
 330   (let* (($ratfac nil)
 331          ($ratprint nil)
 332          (rat ($rat p)))
 333     (and (eq (caar rat) 'mrat) ; don't handle bags etc.
 334          (not (memq 'trunc (car rat))) ; don't handle taylor series (even in other vars)
 335          (let* (;; in affine poly, only numer can include vars
 336                 (num ($ratnumer rat))
 337                 ;; (what vars are actually used; cf. $ratfreeof/$showratvars)
 338                 (numvars (caddar (minimize-varlist num)))
 339                 ;; ... and denominator cannot depend on vars at all
 340                 (den ($ratdenom rat))
 341                 (denvars (caddar (minimize-varlist den)))
 342                 (truncnum))
 343            (and
 344             ;; everything in denvars must be freeof vl
 345             (every #'(lambda (term)
 346                        (every #'(lambda (var) (freeof var term)) vl))
 347                    denvars)
 348             ;; everything in numvars must be either in vl or freeof vl
 349             (every #'(lambda (term)
 350                        (or (memalike term vl)
 351                            (every #'(lambda (var) (freeof var term)) vl)))
 352                    numvars)
 353             ;; there must be no terms of degree > 1
 354             (progn
 355               ;; calculate p without terms of degree > 1
 356               (let (($ratwtlvl 1)
 357                     ;; ignore prevailing *ratweights (don't append to new ones)
 358                     (*ratweights (mapcar #'(lambda (x) (cons x 1)) vl)))
 359                 ;; adding ($rat 0) performs the truncation; just ($rat num) does not
 360                 (setq truncnum (add ($rat 0) num)))
 361               ;; subtract out: any terms of degree > 1?
 362               (equal 0 (ratdisrep (sub num truncnum)))))))))
 363
 364 (defun linear-elimination (l v)
 365   (let (($linsolve_params nil) ($backsubst t) ($programmode t)
 366               ($linsolvewarn nil) ($globalsolve nil) (subs) (unsolved) (vars))
 367
 368     (setq l ($elim l v))
 369     (cond (($member 1 ($first l)) '$emptyset)
 370           (t
 371            (setq subs ($linsolve ($second l) v))
 372            (setq vars (mapcar '$lhs (margs subs)))
 373            (setq vars (push '(mlist) vars))
 374      (setq unsolved ($first l))
 375      (setq unsolved (cons '(mlist) (mapcar #'(lambda (q) (take '(mequal) q 0)) (cdr ($first l)))))
 376      (simplifya (list '(mlist) subs unsolved vars) t)))))
 377
 378 (defun $fourier_elim (l vars)
 379
 380   (let ((eq-list nil) (pos-list nil) (other-list nil) (acc) ($listconstvars nil)
 381         ($ratprint nil) ($domain '$real))
 382
 383     ;; Check the arguments
 384
 385     (setq l (if (op-equalp l 'mand 'mlist) (margs l) (list l)))
 386     (require-list-or-set vars "$fourier_elim")
 387
 388     ;; Standardize each inequality and simplify. To simplify, apply 'mand. After that,
 389     ;; expand 'and' over 'or.'
 390
 391     (setq l (opapply 'mand (mapcar #'(lambda (s) (standardize-inequality s)) l)))
 392     (setq l (expand-and-over-or l))
 393
 394     (cond ((eq t l) '$universalset)
 395           ((eq nil l) '$emptyset)
 396           ((op-equalp l 'mor)
 397            (setq l (margs l))
 398            (dolist (li l)
 399              (push ($fourier_elim li vars) acc))
 400
 401            (setq acc (delete '$emptyset acc))
 402            (cond ((null acc) '$emptyset)
 403                  ((member t acc) '$universalset)
 404                  (t (opapply 'mor acc))))
 405
 406           (t
 407            (setq l (if (op-equalp l 'mand) (margs l) (list l)))
 408            (setq l (delete t l))
 409
 410            ;; Push all non-affine expressions into other-list; push all equalities
 411            ;; into eq-list, push all > equalities into pos-list, and push everything else
 412            ;; into other-list.
 413
 414            (dolist (li l)
 415              (cond ((not ($affine_p ($lhs li) ($listofvars li))) (push li other-list))
 416                    ((op-equalp li 'mequal) (push li eq-list))
 417                    ((op-equalp li 'mgreaterp) (push ($lhs li) pos-list))
 418                    (t (push li other-list))))
 419
 420            ;; Using eq-list, eliminate as many variables as possible.
 421
 422            (push '(mlist) eq-list)
 423            (push '(mlist) pos-list)
 424            (push '(mlist) other-list)
 425            (setq eq-list (linear-elimination eq-list vars))
 426
 427            (cond ((eq '$emptyset eq-list) (setq pos-list '$emptyset))
 428                  (t
 429                   ;;;(setq elim-vars ($third eq-list))
 430                   (setq other-list (append other-list (margs ($second eq-list))))
 431                   (setq eq-list ($first eq-list))
 432                   (setq pos-list ($substitute eq-list pos-list))
 433                   (setq other-list ($substitute eq-list other-list))))
 434
 435            (cond ((eq pos-list '$emptyset) pos-list)
 436                  (t
 437                   (setq pos-list (fourier-elim (margs pos-list) (margs vars)))
 438                   (if (eq pos-list '$emptyset) pos-list
 439                     ($append eq-list pos-list other-list))))))))
 440
 441 ;; Without the post-fourier-elim-simp, we get
 442
 443 ;; (%i2) eqs : [0 < x, x<1, 0 < y, y <1, x+y+z < 4, z > 0]$
 444 ;; (%i3) fourier_elim(eqs,[z,y,x]);
 445 ;; (%o3) [0<z,z<-y-x+4,0<y,y<min(1,4-x),0<x,x<1]
 446
 447 ;; Since 0 < x < 1, Maxima should be able to deduce that min(1,4-x) = 1. So (%o3)
 448 ;; should simplify to [0<z,z<-y-x+4,0<y,y<1,0<x,x<1].
 449 ;; To do this simplification, it's necessary to set dosimp to true--why, I don't
 450 ;; know (alternatively, use ($expand e 0 0)).
 451
 452 (defun post-fourier-elim-simp (pos)
 453   (let ((save-context $context) (new-context (gensym)))
 454
 455     ;; (a) Declare a new context (b) put every member of the CL list
 456     ;; pos into the fact database (c) simplify each member of pos
 457     ;; (d) kill the new context and restore the old context.
 458
 459     (unwind-protect
 460         (progn
 461           ($newcontext new-context)
 462           (mapcar #'(lambda (s) (mfuncall '$assume s)) pos)
 463           (let ((dosimp t))
 464             (setq pos (mapcar #'(lambda (s) (simplifya s nil)) pos))
 465             (delete-if #'(lambda (s) (equal t (standardize-inequality s))) pos)))
 466
 467       (if ($member new-context $contexts) ($killcontext new-context))
 468       (setq $context save-context))))
 469
 470 (defun fourier-elim (pos vars)
 471   (let ((vi) (acc nil) (w))
 472     (while (and (not (eq pos '$emptyset)) (first pos) vars)
 473       (setq vi (pop vars))
 474
 475       (setq w (fourier-elim-one-variable pos vi))
 476       (setq acc (append acc (first w)))
 477       (setq pos (second w)))
 478     (if (eq '$emptyset pos) '$emptyset
 479       (progn
 480         (setq pos (delete t (mapcar #'(lambda (s) (fe> s 0)) pos)))
 481         (if (consp (member 'nil pos)) '$emptyset
 482                 (opapply 'mlist (post-fourier-elim-simp (append acc pos))))))))
 483
 484 ;; Do Fourier elimination on a list l of inequalities with respect to the
 485 ;; indeterminate x. Each member of l is a polynomial, and each polynomial p in
 486 ;; in this list  implies the inequality p > 0. The function fourier-elim-one-variable
 487 ;; returns a list of the form
 488 ;;
 489 ;;   ((lb < x , x < ub), other),
 490 ;;
 491 ;; where other list of list of the implied inequalities (require that each lower bound is less than
 492 ;; each upper bound) and members of p that do not determine a bound on v.
 493
 494 ;; Each member of l is a polynomial with real coefficients. Also each member p of l
 495 ;; determines an inequality p > 0.
 496
 497 (defun simplify-fe-args (pos)
 498
 499   ;; (1) Delete all manifestly positive members of pos.
 500   ;; (2) If pos has a zero or negative member set pos to (-1).
 501   ;; (3) Remove the extremal members of pos (apply min).
 502
 503   (cond ((member pos '($inf $minf) :test #'eq) (list pos))
 504         (t
 505          (let ((save-context $context) (new-context (gensym)))
 506            (unwind-protect
 507                (progn
 508                  ($newcontext new-context)
 509                  (setq pos (delete-if #'(lambda (s) (eq '$pos (csign s))) pos))
 510                  (if (some #'(lambda (s) (member (csign s) '($neg $nz $zero) :test #'eq)) pos) (setq pos (list -1)))
 511                  (setq pos (apply-min pos))
 512                  (if (op-equalp pos '$min) (margs pos) (list pos)))
 513              (if ($member new-context $contexts) ($killcontext new-context))
 514              (setq $context save-context))))))
 515
 516 (defun fourier-elim-one-variable (pos x &optional (need-simp t))
 517
 518   (let ((cf) (acc nil) (lb nil) (ub nil) (lb-args) (ub-args) (sgn) (sol) (bounds nil))
 519         (if need-simp (setq pos (simplify-fe-args pos)))
 520
 521     ;; Find the lower and upper bounds for x; push the lower bounds into the list 'lb', and the upper
 522     ;; bounds into holds 'ub' holds the upper bounds. All polynomials that don't determine a bound get
 523     ;; pushed into acc.
 524
 525     ;; Example: say 1 - x > 0. The coefficient of x is negative, and solving 1 - x = 0 for
 526     ;; x gives an *upper* bound for x (that is x < 1) . Simiarly, say 1 + x > 0. The
 527     ;; coefficient of x is positive and solving 1 + x = 0 for x gives a *lower* bound for x
 528     ;; (that is x > -1).
 529
 530     (dolist (li pos)
 531       (setq cf ($ratcoef li x))
 532       (setq sgn (number-sign cf))
 533       (if (member sgn '($neg $pos) :test #'eq)
 534           (progn
 535             (setq sol ($expand (div (sub (mul cf x) li) cf)))
 536             (if (eq sgn '$neg) (push sol ub) (push sol lb)))
 537         (push li acc)))
 538
 539     ;; Delete the non-extremal members of lb and ub
 540
 541     (setq lb (apply-max lb))
 542     (setq ub (apply-min ub))
 543     (setq lb-args (if (op-equalp lb '$max) (margs lb) (list lb)))
 544     (setq ub-args (if (op-equalp ub '$min) (margs ub) (list ub)))
 545
 546     (setq lb-args (delete '$minf lb-args))
 547     (setq ub-args (delete '$inf ub-args))
 548
 549     ;; Require that each lower bound be smaller than each upper bound. We could use
 550     ;; fe> instead of sub here. But I think it's not needed, and fe> is more spendy.
 551
 552     (setq acc (append acc (mapcar #'(lambda (s) (sub (second s) (first s)))
 553                                   (fourier_elim-cartesian-product (list lb-args ub-args)))))
 554
 555     ;; Return ((lb < x x < ub) acc).
 556
 557     (if (not (eq '$inf ub)) (push (opcons 'mlessp x ub) bounds))
 558     (if (not (eq '$minf lb)) (push (opcons 'mlessp lb x) bounds))
 559
 560     (list
 561      ;;(list (opcons 'mlessp lb x) (opcons 'mlessp x ub))
 562      bounds
 563      (if (some #'(lambda (s) (or (eql 0 s) (eq nil s))) acc) '$emptyset acc))))
 564
 565 ;; Apply max and min without looking at the current context; if something goes wrong,
 566 ;; cleanup the mess.
 567
 568 (defun apply-max (l)
 569   (let ((save-context $context) (new-context (gensym)))
 570     (unwind-protect
 571         (progn
 572           ($newcontext new-context)
 573           (setq l (mapcar '$expand l))
 574           (opapply '$max l))
 575       (if ($member new-context $contexts) ($killcontext new-context))
 576       (setq $context save-context))))
 577
 578 (defun apply-min (l)
 579   (let ((save-context $context) (new-context (gensym)))
 580     (unwind-protect
 581       (progn
 582         ($newcontext new-context)
 583         (setq l (mapcar '$expand l))
 584         (opapply '$min l))
 585       (if ($member new-context $contexts) ($killcontext new-context))
 586       (setq $context save-context))))
 587
 588 (defun compare-using-empty-context (a b)
 589   (let ((save-context $context) (new-context (gensym)))
 590     (unwind-protect
 591         (progn
 592           ($newcontext new-context)
 593           ($compare a b))
 594       (if ($member new-context $contexts) ($killcontext new-context))
 595       (setq $context save-context))))
 596