Ticket #3431: trac772-qepcad-interface.patch

a/sage/interfaces/all.py

@@ -21 +21 @@
-from mupad import mupad, mupad_console, Mupad  # NOT functional yet
-from mwrank import mwrank, Mwrank, mwrank_console
-from octave import octave, octave_console, octave_version, Octave
+from qepcad import qepcad, qepcad_console, qepcad_version, qepcad_formula
-from qsieve import qsieve
-from singular import singular, singular_console, singular_version, is_SingularElement, Singular
-from sage0 import sage0 as sage0, sage0_console, sage0_version, Sage

/dev/null

@@ +1 @@
+r"""
+Interface to QEPCAD B
+
+The basic function of QEPCAD is to construct cylindrical algebraic
+decompositions (CADs) of $\RR^k$, given a list of polynomials.  Using
+this CAD, it is possible to perform quantifier elimination and formula
+simplification.
+
+A CAD for a set $A$ of $k$-variate polynomials decomposes $\RR^j$ into
+disjoint cells, for each $j$ in $0 \leq j \leq k$.  The sign of each
+polynomial in $A$ is constant in each cell of $\RR^k$, and for each
+cell in $\RR^j$ ($j > 1$), the projection of that cell into
+$\RR^{j-1}$ is a cell of $\RR^{j-1}$.  (This property makes the
+decomposition ``cylindrical''.)
+
+Given a formula $\exists x. P(a,b,x) = 0$ (for a polynomial $P$), and
+a cylindrical algebraic decomposition for $P$, we can eliminate the
+quantifier (find an equivalent formula in the two variables $a$, $b$
+without the quantifier $\exists$) as follows.  For each cell $C$ in
+$\RR^2$, find the cells of $\RR^3$ which project to $C$.  (This
+collection is called the \emph{stack} over $C$.)  Mark $C$ as true if
+some member of the stack has sign $= 0$; otherwise, mark $C$ as false.
+Then, construct a polynomial formula in $a$, $b$ which specifies
+exactly the true cells (this is always possible).  The same technique
+works if the body of the quantifier is any boolean combination of
+polynomial equalities and inequalities.
+
+Formula simplification is a similar technique.  Given a formula which
+describes a simple set of $\RR^k$ in a complicated way as a boolean
+combination of polynomial equalities and inequalities, QEPCAD can
+construct a CAD for the polynomials and recover a simple equivalent
+formula.
+
+Note that while the following documentation is tutorial in nature, and
+is written for people who may not be familiar with QEPCAD, it is
+documentation for the \sage interface rather than for QEPCAD.  As
+such, it does not cover several issues that are very important to use
+QEPCAD efficiently, such as variable ordering, the efficient use of
+the alternate quantifiers and \code{_root_} expressions, the
+\code{measure-zero-error} command, etc.  For more information on
+QEPCAD, see the online documentation at
+\url{http://www.cs.usna.edu/~qepcad/B/QEPCAD.html} and Chris Brown's
+tutorial handout and slides from
+\url{http://www.cs.usna.edu/~wcbrown/research/ISSAC04/Tutorial.html}.
+(Several of the examples in this documentation came from these
+sources.)
+
+The examples below require that the optional qepcad package is installed.
+
+QEPCAD can be run in a fully automatic fashion, or interactively.
+We first demonstrate the automatic use of QEPCAD.
+
+Since \sage has no built-in support for quantifiers, this interface
+provides \code{qepcad_formula} which helps construct quantified
+formulas in the syntax QEPCAD requires.
+
+sage: var('a,b,c,d,x,y,z')
+(a, b, c, d, x, y, z)
+sage: qf = qepcad_formula
+
+We start with a simple example.  Consider an arbitrarily-selected
+ellipse:
+
+sage: ellipse = 3*x^2 + 2*x*y + y^2 - x + y - 7
+
+What is the projection onto the $x$ axis of this ellipse?  First we 
+construct a formula asking this question.
+
+sage: F = qf.exists(y, ellipse == 0); F
+(E y)[y^2 + 2 x y + y + 3 x^2 - x - 7 = 0]
+
+Then we run qepcad to get the answer.
+
+sage: qepcad(F)
+8 x^2 - 8 x - 29 <= 0
+
+How about the projection onto the $y$ axis?
+
+sage: qepcad(qf.exists(x, ellipse == 0))
+8 y^2 + 16 y - 85 <= 0
+
+QEPCAD deals with more quantifiers than just ``exists'', of course.
+Besides the standard ``forall'', there are also ``for infinitely
+many'', ``for all but finitely many'', ``for a connected subset'', and
+``for exactly $k$''.  The \function{qepcad} documentation has examples
+of all of these; here we'll just give one example.  
+
+First we construct a circle:
+
+sage: circle = x^2 + y^2 - 3
+
+For what values $k$ does a vertical line $x=k$ intersect the combined
+figure of the circle and ellipse exactly three times?
+
+sage: F = qf.exactly_k(3, y, circle * ellipse == 0); F
+(X3 y)[(y^2 + x^2 - 3) (y^2 + 2 x y + y + 3 x^2 - x - 7) = 0]
+sage: qepcad(F)
+8 x^2 - 8 x - 29 <= 0 /\ x^2 - 3 <= 0 /\ 8 x^4 - 26 x^2 - 4 x + 13 >= 0 /\ [ 8 x^4 - 26 x^2 - 4 x + 13 = 0 \/ x^2 - 3 = 0 \/ 8 x^2 - 8 x - 29 = 0 ]
+
+Here we see that the solutions are among the eight ($4 + 2 + 2$) roots
+of the three polynomials inside the brackets, but not all of these
+roots are solutions; the polynomial inequalities outside the brackets
+are needed to select those roots that are solutions.
+
+QEPCAD also supports an extended formula language, where
+$\mbox{_root_}k \quad P(\bar x, y)$ refers to a particular zero of
+$P(\bar x, y)$ (viewed as a polynomial in $y$).  If there are $n$ roots,
+then _root_$1$ refers to the least root and _root_$n$ refers to the greatest; 
+also, _root_$-n$ refers to the least root and _root_$-1$ refers to the 
+greatest.
+
+This extended language is available both on input and output; see the
+QEPCAD documentation for more information on how to use this syntax on
+input.  We can request output that's intended to be easy to interpret
+geometrically; then QEPCAD will use the extended language to produce
+a solution formula without the selection polynomials.
+
+sage: qepcad(F, solution='geometric')
+x = _root_1 8 x^2 - 8 x - 29
+\/
+8 x^4 - 26 x^2 - 4 x + 13 = 0
+\/
+x = _root_-1 x^2 - 3
+
+We then see that the 6 solutions correspond to the vertical tangent on
+the left side of the ellipse, the four intersections between the
+ellipse and the circle, and the vertical tangent on the right side of
+the circle.
+
+Let's do some basic formula simplification and visualization.
+We'll look at the region which is inside both the ellipse and the circle:
+
+sage: F = qf.and_(ellipse < 0, circle < 0); F
+[y^2 + 2 x y + y + 3 x^2 - x - 7 < 0 /\ y^2 + x^2 - 3 < 0]
+sage: qepcad(F)
+y^2 + 2 x y + y + 3 x^2 - x - 7 < 0 /\ y^2 + x^2 - 3 < 0
+
+We get back the same formula we put in.  This is not surprising (we
+started with a pretty simple formula, after all), but it's not very
+enlightening either.  Again, if we ask for a 'geometric' output, then we see
+an output that lets us understand something about the shape of the solution 
+set.
+
+sage: qepcad(F, solution='geometric')
+[
+  [
+    x = _root_-2 8 x^4 - 26 x^2 - 4 x + 13
+    \/
+    x = _root_2 8 x^4 - 26 x^2 - 4 x + 13
+    \/
+    8 x^4 - 26 x^2 - 4 x + 13 < 0
+  ]
+  /\
+  y^2 + 2 x y + y + 3 x^2 - x - 7 < 0
+  /\
+  y^2 + x^2 - 3 < 0
+]
+\/
+[
+  x > _root_2 8 x^4 - 26 x^2 - 4 x + 13
+  /\
+  x < _root_-2 8 x^4 - 26 x^2 - 4 x + 13
+  /\
+  y^2 + x^2 - 3 < 0
+]
+
+There's another reason to prefer output using _root_ expressions; not
+only does it sometimes give added insight into the geometric
+structure, it also can be more efficient to construct.  Consider this
+formula for the projection of a particular semicircle onto the $x$
+axis:
+
+sage: F = qf.exists(y, qf.and_(circle == 0, x + y > 0)); F
+(E y)[y^2 + x^2 - 3 = 0 /\ y + x > 0]
+sage: qepcad(F)
+x^2 - 3 <= 0 /\ [ x > 0 \/ 2 x^2 - 3 < 0 ]
+
+Here, the formula $x > 0$ had to be introduced in order to get a
+solution formula; the original CAD of F did not include the 
+polynomial $x$.  To avoid having QEPCAD do the extra work to come up 
+with a solution formula, we can tell it to use the extended language; 
+it's always possible to construct a solution formula in the extended 
+language without introducing new polynomials.
+
+sage: qepcad(F, solution='extended')
+x^2 - 3 <= 0 /\ x > _root_1 2 x^2 - 3
+
+Up to this point, all the output we've seen has basically been in the
+form of strings; there is no support (yet) for parsing these outputs
+back into \sage polynomials (partly because \sage does not yet have
+support for symbolic conjunctions and disjunctions).  The function 
+\function{qepcad} supports three more output types that give numbers which
+can be manipulated in \sage: any-point, all-points, and cell-points.
+
+These output types give dictionaries mapping variable names to values.
+With any-point, \function{qepcad} either produces a single dictionary
+specifying a point where the formula is true, or raises an exception
+if the formula is false everywhere.  With all-points,
+\function{qepcad} either produces a list of dictionaries for all
+points where the formula is true, or raises an exception if the
+formula is true on infinitely many points.  With cell-points,
+\function{qepcad} produces a list of dictionaries with one point for
+each cell where the formula is true.  (This means you will have at least
+one point in each connected component of the solution, although you will
+often have many more points than that.)
+
+Let's revisit some of the above examples and get some points to play
+with.  We'll start by finding a point on our ellipse.
+
+sage: p = qepcad(ellipse == 0, solution='any-point'); p
+{'y': [0.96850196850295267 .. 0.96850196850295279], 'x': [-1.4685019685029528 .. -1.4685019685029525]}
+
+(Note that despite appearances, these are really exact numbers.)
+
+We can verify that this point is a solution.  To do so, we create
+a copy of ellipse as a polynomial over QQ (instead of a symbolic 
+expresision).
+
+sage: pellipse = QQ['x,y'](ellipse)
+sage: pellipse(**p) == 0
+True
+
+For cell-points, let's look at points \emph{not} on the ellipse.
+
+sage: pts = qepcad(ellipse != 0, solution='cell-points'); pts
+[{'y': 0, 'x': 4}, {'y': 1, 'x': [2.4685019685029523 .. 2.4685019685029528]}, {'y': -9, 'x': [2.4685019685029523 .. 2.4685019685029528]}, {'y': 9, 'x': 1/2}, {'y': -1, 'x': 1/2}, {'y': -5, 'x': 1/2}, {'y': 3, 'x': [-1.4685019685029528 .. -1.4685019685029525]}, {'y': -1, 'x': [-1.4685019685029528 .. -1.4685019685029525]}, {'y': 0, 'x': -3}]
+
+For the points here which are in full-dimensional cells, QEPCAD has the
+freedom to choose rational sample points, and it does so.
+
+And, of course, all these points really are not on the ellipse.
+
+sage: [pellipse(**p) != 0 for p in pts]
+[True, True, True, True, True, True, True, True, True]
+
+Finally, for all-points, let's look again at finding vertical lines that 
+intersect the union of the circle and the ellipse exactly three times.
+
+sage: F = qf.exactly_k(3, y, circle * ellipse == 0); F
+(X3 y)[(y^2 + x^2 - 3) (y^2 + 2 x y + y + 3 x^2 - x - 7) = 0]
+sage: pts = qepcad(F, solution='all-points'); pts
+[{'x': [1.7320508075688771 .. 1.7320508075688775]}, {'x': [1.7310549134625334 .. 1.7310549134625338]}, {'x': [0.67891138420800389 .. 0.67891138420800401]}, {'x': [-0.94177273774171678 .. -0.94177273774171665]}, {'x': [-1.4681935599288210 .. -1.4681935599288207]}, {'x': [-1.4685019685029528 .. -1.4685019685029525]}]
+
+Since $y$ is bound by the quantifier, the solutions only refer to $x$.
+
+We can substitute one of these solutions into the original equation:
+
+sage: pt = pts[0]
+sage: pcombo = QQ['x,y'](circle * ellipse)
+sage: intersections = pcombo(y=polygen(AA, 'y'), **pt); intersections
+y^4 + [4.4641016151377543 .. 4.4641016151377553]*y^3 + [0.26794919243112269 .. 0.26794919243112276]*y^2
+
+and verify that it does have three roots:
+
+sage: intersections.roots()
+[([-4.4032490056009586 .. -4.4032490056009576], 1), ([-0.060852609536796533 .. -0.060852609536796525], 1), ([0.00000000000000000 .. 0.00000000000000000], 2)]
+
+Let's check all six solutions.
+
+sage: [len(pcombo(y=polygen(AA, 'y'), **p).roots()) for p in pts]
+[3, 3, 3, 3, 3, 3]
+
+We said earlier that we can run QEPCAD either automatically or 
+interactively.  Now that we've discussed the automatic modes, let's turn 
+to interactive uses.
+
+If the \function{qepcad} function is passed \code{interact=True}, then
+instead of returning a result, it returns an object of class
+\class{Qepcad} representing a running instance of QEPCAD that you can
+interact with.  For example:
+
+sage: qe = qepcad(qf.forall(x, x^2 + b*x + c > 0), interact=True); qe
+QEPCAD object in phase 'Before Normalization'
+
+This object is a fairly thin wrapper over QEPCAD; most QEPCAD commands
+are available as methods on the \class{Qepcad} object.  Given a
+\class{Qepcad} object \var{qe}, you can type \code{qe.[tab]} to see
+the available QEPCAD commands; to see the documentation for an
+individual QEPCAD command, for example \code{d_setting}, you can type
+\code{qe.d_setting?}.  (In QEPCAD, this command is called
+\code{d-setting}; we systematically replace hyphens with underscores
+for this interface.)
+
+The execution of QEPCAD is divided into four phases; most commands are
+not available during all phases.  We saw above that QEPCAD starts out
+in phase `Before Normalization'; we see that the \code{d_cell} command
+is not available in this phase:
+
+sage: qe.d_cell()
+Error GETCID: This command is not active here.
+
+We will focus here on the fourth (and last) phase, `Before Solution',
+because this interface has special support for some operations in this
+phase.  Consult the QEPCAD documentation for information on the other
+phases.
+
+We can tell QEPCAD to finish off the current phase and move to the
+next with its \code{go} command.  (There is also the \code{step}
+command, which partially completes a phase for phases that have
+multiple steps, and the \code{finish} command, which runs QEPCAD to
+completion.)
+
+sage: qe.go()
+QEPCAD object has moved to phase 'Before Projection (x)'
+sage: qe.go()
+QEPCAD object has moved to phase 'Before Choice'
+sage: qe.go()
+QEPCAD object has moved to phase 'Before Solution'
+
+Note that the \class{Qepcad} object returns the new phase whenever the
+phase changes, as a convenience for interactive use; except that when
+the new phase is `EXITED', the solution formula printed by QEPCAD is
+returned instead.
+
+sage: qe.go()
+4 c - b^2 > 0
+sage: qe
+QEPCAD object in phase 'EXITED'
+
+Let's pick a nice, simple example, return to phase 4, and explore the
+resulting \var{qe} object.
+
+sage: qe = qepcad(circle == 0, interact=True); qe
+QEPCAD object in phase 'Before Normalization'
+sage: qe.go(); qe.go(); qe.go()
+QEPCAD object has moved to phase 'Before Projection (y)'
+QEPCAD object has moved to phase 'Before Choice'
+QEPCAD object has moved to phase 'Before Solution'
+
+We said before that QEPCAD creates ``cylindrical algebraic
+decompositions''; since we have a bivariate polynomial, we get
+decompositions of $\RR^0$, $\RR^1$, and $\RR^2$.  In this case, where
+our example is a circle of radius $\sqrt{3}$ centered on the origin,
+these decompositions are as follows:
+
+The decomposition of $\RR^0$ is trivial (of course).  The
+decomposition of $\RR^1$ has five cells: $x < -\sqrt{3}$, $x =
+-\sqrt{3}$, $-\sqrt{3} < x < \sqrt{3}$, $x = \sqrt{3}$, and $x >
+\sqrt{3}$.  These five cells comprise the \emph{stack} over the single
+cell in the trivial decomposition of $\RR^0$.
+
+These five cells give rise to five stacks in $\RR^2$.  The first and
+fifth stack have just one cell apiece.  The second and fourth stacks
+have three cells: $y < 0$, $y = 0$, and $y > 0$.  The third stack has
+five cells: below the circle, the lower semicircle, the interior of
+the circle, the upper semicircle, and above the circle.
+
+QEPCAD (and this QEPCAD interface) number the cells in a stack
+starting with 1.  Each cell has an \emph{index}, which is a tuple of
+integers describing the path to the cell in the tree of all cells.
+For example, the cell ``below the circle'' has index (3,1) (the first
+cell in the stack over the third cell of $\RR^1$) and the interior of
+the circle has index (3,3).
+
+We can view these cells with the QEPCAD command \code{d_cell}.  For
+instance, let's look at the cell for the upper semicircle:
+
+sage: qe.d_cell(3, 4)
+---------- Information about the cell (3,4) ----------
+Level                       : 2
+Dimension                   : 1
+Number of children          : 0
+Truth value                 : T    by trial evaluation.
+Degrees after substitution  : Not known yet or No polynomial.
+Multiplicities              : ((1,1))
+Signs of Projection Factors
+Level 1  : (-)
+Level 2  : (0)
+----------   Sample point  ---------- 
+The sample point is in a PRIMITIVE representation.
+<BLANKLINE>
+alpha = the unique root of x^2 - 3 between 0 and 4
+      = 1.7320508076-
+<BLANKLINE>
+Coordinate 1 = 0
+             = 0.0000000000
+Coordinate 2 = alpha
+             = 1.7320508076-
+----------------------------------------------------
+
+We see that, the level of this cell is 2, meaning that it is part of
+the decomposition of $\RR^2$.  The dimension is 1, meaning that the
+cell is homeomorphic to a line (rather than a plane or a point).  The
+sample point gives the coordinates of one point in the cell, both
+symbolically and numerically.
+
+For programmatic access to cells, we've defined a \sage wrapper class
+\class{QepcadCell}.  These cells can be created with the \method{cell}
+method; for example:
+
+sage: c = qe.cell(3, 4); c
+QEPCAD cell (3, 4)
+
+A \class{QepcadCell} has accessor methods for the important state held
+within a cell.  For instance:
+
+sage: c.level()
+2
+sage: c.index()
+(3, 4)
+sage: qe.cell(3).number_of_children()
+5
+sage: len(qe.cell(3))
+5
+
+One particularly useful thing we can get from a cell is its sample point,
+as \sage algebraic real numbers.
+
+sage: c.sample_point()
+(0, [1.7320508075688771 .. 1.7320508075688775])
+sage: c.sample_point_dict()
+{'y': [1.7320508075688771 .. 1.7320508075688775], 'x': 0}
+
+We've seen that we can get cells using the \method{cell} method.
+There are several QEPCAD commands that print lists of cells; we can
+also get cells using the \method{make_cells} method, passing it the
+output of one of these commands.
+
+sage: qe.make_cells(qe.d_true_cells())
+[QEPCAD cell (4, 2), QEPCAD cell (3, 4), QEPCAD cell (3, 2), QEPCAD cell (2, 2)]
+
+Also, the cells in the stack over a given cell can be accessed using
+array subscripting or iteration.  (Remember that cells in a stack are
+numbered starting with one; we preserve this convention in the
+array-subscripting syntax.)
+
+sage: c = qe.cell(3)
+sage: c[1]
+QEPCAD cell (3, 1)
+sage: [c2 for c2 in c]
+[QEPCAD cell (3, 1), QEPCAD cell (3, 2), QEPCAD cell (3, 3), QEPCAD cell (3, 4), QEPCAD cell (3, 5)]
+
+We can do one more thing with a cell: we can set its truth value.
+Once the truth values of the cells have been set, we can get QEPCAD to
+produce a formula which is true in exactly the cells we have selected.
+This is useful if QEPCAD's quantifier language is insufficient to
+express your problem.
+
+For example, consider again our combined figure of the circle and the
+ellipse.  Suppose you want to find all vertical lines that intersect
+the circle twice, and also intersect the ellipse twice.  The vertical
+lines that intersect the circle twice can be found by simplifying:
+
+sage: F = qf.exactly_k(2, y, circle == 0); F
+(X2 y)[y^2 + x^2 - 3 = 0]
+
+and the vertical lines that intersect the ellipse twice are expressed by:
+
+sage: G = qf.exactly_k(2, y, ellipse == 0); G
+(X2 y)[y^2 + 2 x y + y + 3 x^2 - x - 7 = 0]
+
+and the lines that intersect both figures would be:
+
+sage: qf.and_(F, G)
+Traceback (most recent call last):
+...
+ValueError: QEPCAD formulas must be in prenex (quantifiers outermost) form
+
+...except that QEPCAD does not support formulas like this; in QEPCAD
+input, all logical connectives must be inside all quantifiers.
+
+Instead, we can get QEPCAD to construct a CAD for our combined figure
+and set the truth values ourselves.  (The exact formula we use doesn't
+matter, since we're going to replace the truth values in the cells; we
+just need to use a formula that uses both polynomials.)
+
+sage: qe = qepcad(qf.and_(ellipse == 0, circle == 0), interact=True)
+sage: qe.go(); qe.go(); qe.go()
+QEPCAD object has moved to phase 'Before Projection (y)'
+QEPCAD object has moved to phase 'Before Choice'
+QEPCAD object has moved to phase 'Before Solution'
+
+Now we want to find all cells $c$ in the decomposition of $\RR^1$ such
+that the stack over $c$ contains exactly two cells on the ellipse, and
+also contains exactly two cells on the circle.
+
+Our input polynomials are ``level-2 projection factors'', we see:
+
+sage: qe.d_proj_factors()
+P_1,1  = fac(J_1,1) = fac(dis(A_2,1))
+       = 8 x^2 - 8 x - 29
+P_1,2  = fac(J_1,2) = fac(dis(A_2,2))
+       = x^2 - 3
+P_1,3  = fac(J_1,3) = fac(res(A_2,1|A_2,2))
+       = 8 x^4 - 26 x^2 - 4 x + 13
+A_2,1  = input
+       = y^2 + 2 x y + y + 3 x^2 - x - 7
+A_2,2  = input
+       = y^2 + x^2 - 3
+
+so we can test whether a cell is on the ellipse by checking that the
+sign of the corresponding projection factor is 0 in our cell.
+For instance, the cell (12,2) is on the ellipse:
+
+sage: qe.cell(12,2).signs()[1][0]
+0
+
+So we can update the truth values as desired like this:
+
+sage: for c in qe.cell():
+...       count_ellipse = 0
+...       count_circle = 0
+...       for c2 in c:
+...           count_ellipse += (c2.signs()[1][0] == 0)
+...           count_circle += (c2.signs()[1][1] == 0)
+...       c.set_truth(count_ellipse == 2 and count_circle == 2)
+
+and then we can get our desired solution formula.  (The 'G' stands for
+'geometric', and gives solutions using the same rules as
+\code{solution='geometric'} described above.)
+
+sage: qe.solution_extension('G')
+8 x^2 - 8 x - 29 < 0
+/\
+x^2 - 3 < 0
+
+AUTHORS:
+    -- Carl Witty (2008-03): initial version
+"""
+
+##########################################################################
+#
+#       Copyright (C) 2008 Carl Witty <Carl.Witty@gmail.com>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#
+#                  http://www.gnu.org/licenses/
+#
+##########################################################################
+
+from __future__ import with_statement
+
+import sage.misc.misc
+import pexpect
+import re
+import sys
+
+from sage.misc.flatten import flatten
+from sage.misc.sage_eval import sage_eval
+from sage.misc.preparser import implicit_mul
+
+from expect import Expect, ExpectFunction, AsciiArtString
+
+def _qepcad_cmd(memcells=None):
+    r"""
+    Construct a QEPCAD command line.  Optionally set the
+    number of memory cells to use.
+
+    EXAMPLES:
+        sage: from sage.interfaces.qepcad import _qepcad_cmd
+        sage: _qepcad_cmd()
+        'env qe=.../local/ qepcad '
+        sage: _qepcad_cmd(memcells=8000000)
+        'env qe=.../local/ qepcad +N8000000'
+    """
+    if memcells is not None:
+        memcells_arg = '+N%s' % memcells
+    else:
+        memcells_arg = ''
+    return "env qe=%s qepcad %s"%(sage.misc.misc.SAGE_LOCAL, memcells_arg)
+
+_command_info_cache = None
+
+def _update_command_info():
+    r"""
+    Read the file \file{qepcad.help} to find the list of commands
+    supported by QEPCAD.  Used for tab-completion and documentation.
+
+    EXAMPLES:
+        sage: from sage.interfaces.qepcad import _update_command_info, _command_info_cache
+        sage: _update_command_info() # optional
+        sage: _command_info_cache['approx_precision'] # optional
+        ('46', 'abcde', 'm', 'approx-precision N\n\nApproximate algeraic numbers to N decimal places.\n', None)
+    """
+    global _command_info_cache
+    if _command_info_cache is not None:
+        return
+
+    cache = {}
+
+    with open('%s/bin/qepcad.help'%sage.misc.misc.SAGE_LOCAL) as help:
+        assert(help.readline().strip() == '@')
+
+        while True:
+            cmd_line = help.readline()
+            while len(cmd_line.strip()) == 0:
+                cmd_line = help.readline()
+            cmd_line = cmd_line.strip()
+            if cmd_line == '@@@':
+                break
+
+            (cmd, id, phases, kind) = cmd_line.split()
+            assert(help.readline().strip() == '@')
+
+            help_text = ''
+            help_line = help.readline()
+            while help_line.strip() != '@':
+                help_text += help_line
+                help_line = help.readline()
+
+            # I went through qepcad.help and picked out commands that
+            # I thought might be worth a little extra tweaking...
+            special = None
+            
+            # These commands have been tweaked.
+            if cmd in ['d-all-cells-in-subtree', 'd-cell', 'd-pcad',
+                       'd-pscad', 'd-stack', 'manual-choose-cell']:
+                special = 'cell'
+            if cmd in ['ipfzt', 'rational-sample', 'triv-convert', 'use-db',
+                       'use-selected-cells-cond', 'verbose']:
+                special = 'yn'
+
+            # The tweaking for these commands has not been implemented yet.
+            if cmd in ['selected-cells-cond']:
+                special = 'formula'
+            if cmd in ['ch-pivot', 'rem-pf', 'rem-pj']:
+                special = 'i,j'
+            if cmd in ['d-2d-cad', 'set-truth-value', 'solution-extension']:
+                special = 'interactive'
+            if cmd in ['p-2d-cad', 'trace-alg', 'trace-data']:
+                special = 'optional'
+
+            cmd = cmd.replace('-', '_')
+
+            cache[cmd] = (id, phases, kind, help_text, special)
+
+    _command_info_cache = cache
+            
+_rewrote_qepcadrc = False
+def _rewrite_qepcadrc():
+    r"""
+    Write the file \file{default.qepcadrc} to specify a path to Singular
+    (which QEPCAD uses for Gr\"obner basis computations).
+
+    EXAMPLES:
+        sage: from sage.interfaces.qepcad import _rewrite_qepcadrc
+        sage: _rewrite_qepcadrc()
+        sage: from sage.misc.misc import SAGE_LOCAL
+        sage: open('%s/default.qepcadrc'%SAGE_LOCAL).readlines()[-1]
+        'SINGULAR .../local//bin'
+    """
+    global _rewrote_qepcadrc
+    if _rewrote_qepcadrc: return
+
+    SL = sage.misc.misc.SAGE_LOCAL
+    fn = '%s/default.qepcadrc'%SL
+    text = \
+"""# THIS FILE IS AUTOMATICALLY GENERATED -- DO NOT EDIT
+
+#####################################################
+# QEPCAD rc file.
+# This file allows for some customization of QEPCADB.
+# Right now, the ability to give a path to Singular,
+# so that it gets used for some computer algebra
+# computations is the only feature.
+#####################################################
+SINGULAR %s/bin""" % SL
+
+    open(fn, 'w').write(text)
+
+# QEPCAD does not have a typical "computer algebra system" interaction
+# model.  Instead, you run QEPCAD once for each problem you wish to solve,
+# then interact with it while you solve that problem.
+
+# For this reason, most of the methods on the Expect class are
+# inapplicable.  To avoid confusing the user with useless commands
+# in tab completion, we split the control into two classes: 
+# Qepcad_expect is a subclass of Expect, and is never seen by the user;
+# Qepcad is a wrapper for Qepcad_expect, and is what the user interacts with.
+
+class Qepcad_expect(Expect):
+    r"""
+    The low-level wrapper for QEPCAD.
+    """
+    def __init__(self, memcells=None,
+                 maxread=100000,
+                 logfile=None,
+                 server=None):
+        r"""
+        Initialize a low-level wrapper for QEPCAD.  You can specify
+        the number of memory cells that QEPCAD allocates on startup
+        (which controls the maximum problem size QEPCAD can handle),
+        and specify a logfile.  You can also specify a server, and the
+        interface will run QEPCAD on that server, using ssh.  (UNTESTED)
+
+        EXAMPLES:
+            sage: from sage.interfaces.qepcad import Qepcad_expect
+            sage: Qepcad_expect(memcells=100000, logfile=sys.stdout)
+            Qepcad
+        """
+        _rewrite_qepcadrc()
+        Expect.__init__(self,
+                        name="QEPCAD",
+                        # yuck: when QEPCAD first starts,
+                        # it doesn't give prompts
+                        prompt="\nEnter an .*:\r",
+                        command=_qepcad_cmd(memcells),
+                        maxread=maxread,
+                        server=server,
+                        restart_on_ctrlc=False,
+                        verbose_start=False,
+                        logfile=logfile)
+
+class Qepcad:
+    r"""
+    The wrapper for QEPCAD.
+    """
+
+    def __init__(self, formula,
+                 vars=None, logfile=None, verbose=False,
+                 memcells=None, server=None):
+        r"""
+        Construct a QEPCAD wrapper object.  Requires a formula, which
+        may be a \class{qformula} as returned by the methods of
+        \code{qepcad_formula}, a symbolic equality or inequality, a
+        polynomial $p$ (meaning $p = 0$), or a string, which is passed
+        straight to QEPCAD.
+
+        \var{vars} specifies the variables to use; this gives the variable
+        ordering, which may be very important.  If \var{formula} is
+        given as a string, then \var{vars} is required; otherwise,
+        if \var{vars} is omitted, then a default ordering is used
+        (alphabetical ordering for the free variables).
+
+        A logfile can be specified with \var{logfile}.
+        If \code{verbose=True} is given, then the logfile is automatically
+        set to \code{sys.stdout}, so all QEPCAD interaction is echoed to
+        the terminal.
+
+        You can set the amount of memory that QEPCAD allocates with
+        \var{memcells}, and you can use \var{server} to run QEPCAD on
+        another machine using ssh.  (UNTESTED)
+        
+        Usually you will not call this directly, but use \function{qepcad}
+        to do so.  Check the \function{qepcad} documentation for more
+        information.
+
+        EXAMPLES:
+            sage: from sage.interfaces.qepcad import Qepcad
+            sage: Qepcad(x^2 - 1 == 0) # optional
+            QEPCAD object in phase 'Before Normalization'
+        """
+        self._cell_cache = {}
+
+        if verbose:
+            logfile=sys.stdout
+
+        varlist = None
+        if vars is not None:
+            if isinstance(vars, str):
+                varlist = vars.strip('()').split(',')
+            else:
+                varlist = [str(v) for v in vars]
+
+        if isinstance(formula, str):
+            if varlist is None:
+                raise ValueError, "vars must be specified if formula is a string"
+            
+            free_vars = len(varlist) - formula.count('(')
+        else:
+            formula = qepcad_formula.formula(formula)
+            fvars = formula.vars
+            fqvars = formula.qvars
+            if varlist is None:
+                varlist = sorted(fvars) + fqvars
+            else:
+                # Do some error-checking here.  Parse the given vars
+                # and ensure they match up with the variables in the formula.
+                if frozenset(varlist) != (fvars | frozenset(fqvars)):
+                    raise ValueError, "specified vars don't match vars in formula"
+                if len(fqvars) and varlist[-len(fqvars):] != fqvars:
+                    raise ValueError, "specified vars don't match quantified vars"
+            free_vars = len(fvars)
+            formula = repr(formula)
+
+        self._varlist = varlist
+        self._free_vars = free_vars
+
+        varlist = [v.replace('_', '') for v in varlist]
+        if len(frozenset(varlist)) != len(varlist):
+            raise ValueError, "variables collide after stripping underscores"
+        formula = formula.replace('_', '')
+
+        qex = Qepcad_expect(logfile=logfile)
+        qex._send('[ input from Sage ]')
+        qex._send('(' + ','.join(varlist) + ')')
+        qex._send(str(free_vars))
+        # I hope this prompt is distinctive enough...
+        # if not, we could list all the cases separately
+        qex._change_prompt(['\r\n([^\r\n]*) >\r\n', pexpect.EOF])
+        qex.eval(formula + '.')
+        self._qex = qex
+
+    def __repr__(self):
+        r"""
+        Return a string representation of this \class{Qepcad} object.
+
+        EXAMPLES:
+            sage: qepcad(x - 1 == 0, interact=True) # optional
+            QEPCAD object in phase 'Before Normalization'
+        """
+        return "QEPCAD object in phase '%s'"%self.phase()
+
+    def assume(self, assume):
+        r"""
+        The following documentation is from \file{qepcad.help}:
+
+        Add an assumption to the problem.  These will not be
+        included in the solution formula.
+
+        For example, with input  (E x)[ a x^2 + b x + c = 0],
+        if we issue the command
+
+            assume [ a /= 0 ]
+
+        we'll get the solution formula b^2 - 4 a c >= 0.  Without
+        the assumption we'd get something like [a = 0 /\ b /= 0] \/ 
+        [a /= 0 /\ 4 a c - b^2 <= 0] \/ [a = 0 /\ b = 0 /\ c = 0].
+
+        EXAMPLES:
+            sage: var('a,b,c,x')
+            (a, b, c, x)
+            sage: qf = qepcad_formula
+            sage: qe = qepcad(qf.exists(x, a*x^2 + b*x + c == 0), interact=True) # optional
+            sage: qe.assume(a != 0) # optional
+            sage: qe.finish() # optional
+            4 a c - b^2 <= 0
+        """
+        if not isinstance(assume, str):
+            assume = qepcad_formula.formula(assume)
+            if len(assume.qvars):
+                raise ValueError, "assumptions cannot be quantified"
+            if not assume.vars.issubset(frozenset(self._varlist[:self._free_vars])):
+                raise ValueError, "assumption contains variables not present in formula"
+            assume = repr(assume)
+        assume = assume.replace('_', '')
+        result = self._eval_line("assume [%s]" % assume)
+        if len(result):
+            return AsciiArtString(result)
+
+    def solution_extension(self, kind):
+        r"""
+        The following documentation is modified from \file{qepcad.help}:
+
+        solution-extension x
+ 
+        Use an alternative solution formula construction method.  The
+        parameter x is allowed to be T,E, or G.  If x is T, then a
+        formula in the usual language of Tarski formulas is produced.
+        If x is E, a formula in the language of Extended Tarski formulas
+        is produced.  If x is G, then a geometry-based formula is
+        produced.
+
+        EXAMPLES:
+            sage: var('x,y')
+            (x, y)
+            sage: qf = qepcad_formula
+            sage: qe = qepcad(qf.and_(x^2 + y^2 - 3 == 0, x + y > 0), interact=True) # optional
+            sage: qe.go(); qe.go(); qe.go() # optional
+            QEPCAD object has moved to phase 'Before Projection (y)'
+            QEPCAD object has moved to phase 'Before Choice'
+            QEPCAD object has moved to phase 'Before Solution'
+            sage: qe.solution_extension('E') # optional
+            x > _root_1 2 x^2 - 3 /\ y^2 + x^2 - 3 = 0 /\ [ 2 x^2 - 3 > 0 \/ y = _root_-1 y^2 + x^2 - 3 ]
+            sage: qe.solution_extension('G') # optional
+            [
+              [
+                2 x^2 - 3 < 0
+                \/
+                x = _root_-1 2 x^2 - 3
+              ]
+              /\
+              y = _root_-1 y^2 + x^2 - 3
+            ]
+            \/
+            [
+              x^2 - 3 <= 0
+              /\
+              x > _root_-1 2 x^2 - 3
+              /\
+              y^2 + x^2 - 3 = 0
+            ]            
+            sage: qe.solution_extension('T') # optional
+            y + x > 0 /\ y^2 + x^2 - 3 = 0
+        """
+        if kind == 'I':
+            raise ValueError, "Interactive solution construction not handled by Sage interface"
+        result = self._eval_line('solution-extension %s'%kind)
+        tagline = 'An equivalent quantifier-free formula:'
+        loc = result.find(tagline)
+        if loc >= 0:
+            result = result[loc + len(tagline):]
+        result = result.strip()
+        if len(result):
+            return AsciiArtString(result)
+
+    def set_truth_value(self, index, nv):
+        r"""
+        Given a cell index (or a cell) and an integer, set the truth value
+        of the cell to that integer.  Valid integers are 0 (false), 
+        1 (true), and 2 (undetermined).
+
+        EXAMPLES:
+            sage: qe = qepcad(x == 1, interact=True) # optional
+            sage: qe.go(); qe.go(); qe.go() # optional
+            QEPCAD object has moved to phase 'At the end of projection phase'
+            QEPCAD object has moved to phase 'Before Choice'
+            QEPCAD object has moved to phase 'Before Solution'
+            sage: qe.set_truth_value(1, 1) # optional
+        """
+        index_str = _format_cell_index([index])
+        self._eval_line('set-truth-value\n%s\n%s' % (index_str, nv))
+
+    def phase(self):
+        r"""
+        Return the current phase of the QEPCAD program.
+
+        EXAMPLES:
+            sage: qe = qepcad(x > 2/3, interact=True) # optional
+            sage: qe.phase() # optional
+            'Before Normalization'
+            sage: qe.go() # optional
+            QEPCAD object has moved to phase 'At the end of projection phase'
+            sage: qe.phase() # optional
+            'At the end of projection phase'
+            sage: qe.go() # optional
+            QEPCAD object has moved to phase 'Before Choice'
+            sage: qe.phase() # optional
+            'Before Choice'
+            sage: qe.go() # optional
+            QEPCAD object has moved to phase 'Before Solution'
+            sage: qe.phase() # optional
+            'Before Solution'
+            sage: qe.go() # optional
+            3 x - 2 > 0
+            sage: qe.phase() # optional
+            'EXITED'
+        """
+        match = self._qex.expect().match
+        if match == pexpect.EOF:
+            return 'EXITED'
+        else:
+            return match.group(1)
+
+    def _parse_answer_stats(self):
+        r"""
+        Parse the final string printed by QEPCAD, which should be
+        the simplified quantifier-free formula followed by some
+        statistics.  Return a pair of the formula and the statistics
+        (both as strings).
+
+        EXAMPLES:
+            sage: qe = qepcad(x^2 > 2, interact=True) # optional
+            sage: qe.finish() # optional
+            x^2 - 2 > 0
+            sage: (ans, stats) = qe._parse_answer_stats() # optional
+            sage: ans # optional
+            'x^2 - 2 > 0'
+            sage: stats # random, optional
+            '-----------------------------------------------------------------------------\r\n0 Garbage collections, 0 Cells and 0 Arrays reclaimed, in 0 milliseconds.\r\n492514 Cells in AVAIL, 500000 Cells in SPACE.\r\n\r\nSystem time: 16 milliseconds.\r\nSystem time after the initialization: 4 milliseconds.\r\n-----------------------------------------------------------------------------\r\n'
+        """
+        if self.phase() != 'EXITED':
+            raise ValueError, "QEPCAD is not finished yet"
+        final = self._qex.expect().before
+        match = re.search('\nAn equivalent quantifier-free formula:(.*)\n=+  The End  =+\r\n\r\n(.*)$', final, re.DOTALL)
+
+        if match:
+            return (match.group(1).strip(), match.group(2))
+        else:
+            return (final, '')
+
+    def answer(self):
+        r"""
+        For a QEPCAD instance which is finished, return the
+        simplified quantifier-free formula that it printed just before
+        exiting.
+
+        EXAMPLES:
+            sage: qe = qepcad(x^3 - x == 0, interact=True) # optional
+            sage: qe.finish() # optional
+            x - 1 <= 0 /\ x + 1 >= 0 /\ [ x = 0 \/ x - 1 = 0 \/ x + 1 = 0 ]
+            sage: qe.answer() # optional
+            x - 1 <= 0 /\ x + 1 >= 0 /\ [ x = 0 \/ x - 1 = 0 \/ x + 1 = 0 ]
+        """
+        return AsciiArtString(self._parse_answer_stats()[0])
+
+    def final_stats(self):
+        r"""
+        For a QEPCAD instance which is finished, return the
+        statistics that it printed just before exiting.
+
+        EXAMPLES:
+            sage: qe = qepcad(x == 0, interact=True) # optional
+            sage: qe.finish() # optional
+            x = 0
+            sage: qe.final_stats() # random, optional
+            -----------------------------------------------------------------------------
+            0 Garbage collections, 0 Cells and 0 Arrays reclaimed, in 0 milliseconds.
+            492840 Cells in AVAIL, 500000 Cells in SPACE.
+            System time: 8 milliseconds.
+            System time after the initialization: 4 milliseconds.
+            -----------------------------------------------------------------------------
+        """
+        return AsciiArtString(self._parse_answer_stats()[1])
+
+    def trait_names(self):
+        r"""
+        Returns a list of the QEPCAD commands which are available as
+        extra methods on a \class{Qepcad} object.
+
+        EXAMPLES:
+            sage: qe = qepcad(x^2 < 0, interact=True) # optional
+            sage: len(qe.trait_names()) # optional
+            95
+            sage: 'd_cell' in qe.trait_names() # optional
+            True
+        """
+        _update_command_info()
+        return _command_info_cache.keys()
+
+    def cell(self, *index):
+        r"""
+        Given a cell index, returns a \class{QepcadCell} wrapper for that
+        cell.  Uses a cache for efficiency.
+
+        EXAMPLES:
+            sage: qe = qepcad(x + 3 == 42, interact=True) # optional
+            sage: qe.go(); qe.go(); qe.go() # optional
+            QEPCAD object has moved to phase 'At the end of projection phase'
+            QEPCAD object has moved to phase 'Before Choice'
+            QEPCAD object has moved to phase 'Before Solution'
+            sage: qe.cell(2) # optional
+