Author: Samuel LeliÃ¨vre
Views : 151
Description: Using QEPCAD in a CoCalc Sage worksheet -- examples from the SageMath reference manual: http://doc.sagemath.org/html/en/reference/interfaces/sage/interfaces/qepcad.html
Compute Environment: Ubuntu 18.04 (Deprecated)
# 2016-04-04
# Testing QEPCAD on CoCalc, using the examples from

var('a b c d x y z')

(a, b, c, d, x, y, z)
ellipse = 3*x^2 + 2*x*y + y^2 - x + y - 7

F = qf.exists(y, ellipse == 0); F

(E y)[3 x^2 + 2 x y + y^2 - x + y - 7 = 0]
qepcad(F)

8 x^2 - 8 x - 29 <= 0
qepcad(qf.exists(x, ellipse == 0))

8 y^2 + 16 y - 85 <= 0
circle = x^2 + y^2 - 3

F = qf.exactly_k(3, y, circle * ellipse == 0); F

(X3 y)[(3 x^2 + 2 x y + y^2 - x + y - 7) (x^2 + y^2 - 3) = 0]
qepcad(F)

x^2 - 3 <= 0 /\ 8 x^2 - 8 x - 29 <= 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 ]
F = qf.exactly_k(3, y, circle * ellipse == 0); F

(X3 y)[(3 x^2 + 2 x y + y^2 - x + y - 7) (x^2 + y^2 - 3) = 0]
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 ]
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
F = qf.and_(ellipse < 0, circle < 0); F

[3 x^2 + 2 x y + y^2 - x + y - 7 < 0 /\ x^2 + y^2 - 3 < 0]
qepcad(F)

y^2 + 2 x y + y + 3 x^2 - x - 7 < 0 /\ y^2 + x^2 - 3 < 0
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 ]
F = qf.exists(y, qf.and_(circle == 0, x + y > 0)); F

(E y)[x^2 + y^2 - 3 = 0 /\ x + y > 0]
qepcad(F)

x^2 - 3 <= 0 /\ [ x > 0 \/ 2 x^2 - 3 < 0 ]
qepcad(F, solution='extended')

x^2 - 3 <= 0 /\ x > _root_1 2 x^2 - 3
p = qepcad(ellipse == 0, solution='any-point'); p

{'y': 0.9685019685029527?, 'x': -1.468501968502953?}
pellipse = QQ['x,y'](ellipse)

pellipse(**p) == 0

True
pts = qepcad(ellipse != 0, solution='cell-points'); pts

[{'y': 0, 'x': 4}, {'y': 1, 'x': 2.468501968502953?}, {'y': -9, 'x': 2.468501968502953?}, {'y': 9, 'x': 1/2}, {'y': -1, 'x': 1/2}, {'y': -5, 'x': 1/2}, {'y': 3, 'x': -1.468501968502953?}, {'y': -1, 'x': -1.468501968502953?}, {'y': 0, 'x': -3}]
[pellipse(**p) != 0 for p in pts]

[True, True, True, True, True, True, True, True, True]
F = qf.exactly_k(3, y, circle * ellipse == 0); F

(X3 y)[(3 x^2 + 2 x y + y^2 - x + y - 7) (x^2 + y^2 - 3) = 0]
pts = qepcad(F, solution='all-points'); pts

[{'x': 1.732050807568878?}, {'x': 1.731054913462534?}, {'x': 0.678911384208004?}, {'x': -0.9417727377417167?}, {'x': -1.468193559928821?}, {'x': -1.468501968502953?}]
pt = pts[0]

pcombo = QQ['x,y'](circle * ellipse)

intersections = pcombo(y=polygen(AA, 'y'), **pt); intersections

y^4 + 4.464101615137755?*y^3 + 0.2679491924311227?*y^2
intersections.roots()

[(-4.403249005600958?, 1), (-0.06085260953679653?, 1), (0, 2)]
[len(pcombo(y=polygen(AA, 'y'), **p).roots()) for p in pts]

[3, 3, 3, 3, 3, 3]
qe = qepcad(qf.forall(x, x^2 + b*x + c > 0), interact=True); qe

QEPCAD object in phase 'Before Normalization'
qe.d_cell()

Error GETCID: This command is not active here.
qe.go()

QEPCAD object has moved to phase 'Before Projection (x)'
qe.go()

QEPCAD object has moved to phase 'Before Choice'
qe.go()

QEPCAD object has moved to phase 'Before Solution'
qe.go()

4 c - b^2 > 0
qe

qe = qepcad(circle == 0, interact=True); qe

QEPCAD object in phase 'Before Normalization'
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'
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. alpha = the unique root of x^2 - 3 between 0 and 4 = 1.7320508076- Coordinate 1 = 0 = 0.0000000000 Coordinate 2 = alpha = 1.7320508076- ----------------------------------------------------
c = qe.cell(3, 4); c

c.level()

2
c.index()

(3, 4)
qe.cell(3).number_of_children()

5
len(qe.cell(3))

5
c.sample_point()

(0, 1.732050807568878?)
c.sample_point_dict()

{'y': 1.732050807568878?, 'x': 0}
qe.make_cells(qe.d_true_cells())

c = qe.cell(3)

c[1]

[c2 for c2 in c]

F = qf.exactly_k(2, y, circle == 0); F

(X2 y)[x^2 + y^2 - 3 = 0]
G = qf.exactly_k(2, y, ellipse == 0); G

(X2 y)[3 x^2 + 2 x y + y^2 - x + y - 7 = 0]
qf.and_(F, G)

Error in lines 1-1 Traceback (most recent call last): File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 904, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/interfaces/qepcad.py", line 1962, in and_ formula_strs, vars = self._combine_formulas(formulas) File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/interfaces/qepcad.py", line 1858, in _combine_formulas raise ValueError("QEPCAD formulas must be in prenex" ValueError: QEPCAD formulas must be in prenex (quantifiers outermost) form
qe = qepcad(qf.and_(ellipse == 0, circle == 0), interact=True)

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'
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
qe.cell(12,2).signs()[1][0]

0
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)

qe.solution_extension('G')

8 x^2 - 8 x - 29 < 0 /\ x^2 - 3 < 0
open('%s/default.qepcadrc'%SAGE_LOCAL).readlines()[-1]

'SINGULAR /projects/sage/sage-6.10/local/bin\n'
from sage.interfaces.qepcad import _qepcad_atoms

var('a,b,c,d,x,y,z')

(a, b, c, d, x, y, z)
qf = qepcad_formula

ellipse = 3*x^2 + 2*x*y + y^2 - x + y - 7

circle = x^2 + y^2 - 3

F = qf.exactly_k(3, y, circle * ellipse == 0)

_qepcad_atoms(qepcad(F))

set(['x^2 - 3 = 0', '8 x^2 - 8 x - 29 <= 0', '8 x^4 - 26 x^2 - 4 x + 13 = 0', '8 x^2 - 8 x - 29 = 0', 'x^2 - 3 <= 0', '8 x^4 - 26 x^2 - 4 x + 13 >= 0'])
_qepcad_atoms(qepcad(F, solution='geometric'))

set(['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'])
F = qf.and_(ellipse < 0, circle < 0)

_qepcad_atoms(qepcad(F))

set(['y^2 + 2 x y + y + 3 x^2 - x - 7 < 0', 'y^2 + x^2 - 3 < 0'])
_qepcad_atoms(qepcad(F, solution='geometric'))

set(['8 x^4 - 26 x^2 - 4 x + 13 < 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 + 2 x y + y + 3 x^2 - x - 7 < 0', 'x < _root_-2 8 x^4 - 26 x^2 - 4 x + 13', 'y^2 + x^2 - 3 < 0', 'x = _root_-2 8 x^4 - 26 x^2 - 4 x + 13'])
F = qf.exists(y, qf.and_(circle == 0, x + y > 0))

_qepcad_atoms(qepcad(F))

set(['x > 0', '2 x^2 - 3 < 0', 'x^2 - 3 <= 0'])
_qepcad_atoms(qepcad(F, solution='extended'))

set(['x > _root_1 2 x^2 - 3', 'x^2 - 3 <= 0'])
qe = qepcad(qf.and_(ellipse == 0, circle == 0), interact=True)

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'
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)

_qepcad_atoms(qe.solution_extension('G'))

set(['x^2 - 3 < 0', '8 x^2 - 8 x - 29 < 0'])