CAD Implementation - GSoC 2015

412 views
Skip to first unread message

Luv Agarwal

unread,
Mar 1, 2015, 7:03:42 PM3/1/15
to sy...@googlegroups.com

Hi,

I am Luv Agarwal, a sophomore pursuing b-tech in computer science. From the past 20 days I have been hopping around CAD (cylindrical algebraic decomposition) and Quantifier Elimination basics, algorithms, improvements and other related stuff to successfully and efficiently implement it. Other than these I read about other uses of CAD like Inequalities solving, collision detection which I suppose can be implemented in the course of Quantifier elimination implementation.

It will be really good if SymPy can perform CAD as it can then be used for quantifier elimination, theorem proving in algebraic geometry, evaluating multiple integrals and in assumptions engine. (Source: https://mattpap.github.io/masters-thesis/html/src/conclusions.html#cylindrical-algebraic-decomposition)

I Plan to implement partial cylindrical algebraic decomposition and quantifier elimination (By - George E. Collins and Hoon Hong). I will check if this partial CAD can be used in other CAD applications. If not then I will make sure that I think of an implementation which can perform both partial as well as full CAD.

My current status of quantifier elimination:

I have read the following:
  1. Collin's original paper Quantifier Elimination for Real Closed Fields by Cylindrical Algebraic Decomposition which contains the core algorithms.
  2. Few papers which talk about the usage of CAD.
  3. Hong's implementation of quantifier elimination from partial cylindrical algebraic decomposition that introduces partial cad which is an enhancement of the complete CAD to perform quantifier elimination.
  4. A small but powerful improvement By Hong in a CAD projection operator.

After all this heavy reading I have a good understanding about the working of these algorithms and what will be needed to implement it. (I have implemented a small  component (projection phase) of CAD for testing purposes). Though I haven't thought about it's  implementation details in SymPy, but I have thought in general. I will soon post those details together with my GSoC proposal.

Also, can someone please update me why FOL (GsoC-2014 by Soumya Dipti Biswas) is not yet merged as quantifier elimination demands it.

Algorithm Functionality:

Input: A quantified formula F.
Output: A quantifier free formula F' (in free variables of F) equivalent to F.

Example:
    >>> F = ForAll(x, x*y=0)
    >>> Reduce(F)
    >>> y = 0

It would be a great help if I can get some community feedbacks.

Thanks.

Aaron Meurer

unread,
Mar 1, 2015, 7:35:57 PM3/1/15
to sy...@googlegroups.com
There were some issues with the class names conflicting with existing
names. This is the pull request
https://github.com/sympy/sympy/pull/7608. Finishing up this pull
request would be a good place to start. I haven't heard from Soumya in
a while, so I don't know if he plans to do it.

I would have to go through it again to remember what the exact issues
were. I would start by reading the log on the pull request and the
logs at https://gitter.im/sympy/sympy/SDB-GSoC.

>
> Algorithm Functionality:
>
> Input: A quantified formula F.
> Output: A quantifier free formula F' (in free variables of F) equivalent to
> F.
>
> Example:
> >>> F = ForAll(x, x*y=0)
> >>> Reduce(F)
> >>> y = 0
>
> It would be a great help if I can get some community feedbacks.

Do you know of any existing computer algebra systems that implement
CAD? It's a good idea to look at them, to see how they structure their
API, how they implement things (if they are open source), what
limitations they have, etc.

Aaron Meurer

>
> Thanks.
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/440e0d91-c457-457e-9c1f-2e22832e5e23%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Luv Agarwal

unread,
Mar 1, 2015, 8:11:12 PM3/1/15
to sy...@googlegroups.com
Maple and Mathematica implements it. Also there is an implementation of this partial CAD written in C called as qepcad (code hosted at github) of which sage also provides an interface (http://www.sagemath.org/doc/reference/interfaces/sage/interfaces/qepcad.html).

Thanks

Joachim Durchholz

unread,
Mar 2, 2015, 4:04:08 AM3/2/15
to sy...@googlegroups.com
Am 02.03.2015 um 01:35 schrieb Aaron Meurer:
> There were some issues with the class names conflicting with existing
> names.

I'm making good progress with removing the class registry and should be
finished within two weeks.
So conflicts within the class registry aren't going to be a problem
anymore. I'm not aware of any other showstopper-ish problems with class
name conflicts; I'd still consider them to be quite undesirable because
they can confuse readers, users, and maintainers.

Aaron Meurer

unread,
Mar 2, 2015, 11:52:38 AM3/2/15
to sy...@googlegroups.com
OK, so the technical problem will be going away. That pull request
really should be reusing classes from the core and the rest of the
logic module, though.

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/54F42784.1070103%40durchholz.org.

Joachim Durchholz

unread,
Mar 2, 2015, 12:23:24 PM3/2/15
to sy...@googlegroups.com
Am 02.03.2015 um 17:52 schrieb Aaron Meurer:
> That pull request
> really should be reusing classes from the core and the rest of the
> logic module, though.

Reusing classes as far as possible is really important.
Every wheel reinvented is another stone added to maintenance burden.
(Sorry for the mixed metaphor.)

Luv Agarwal

unread,
Mar 6, 2015, 9:28:17 PM3/6/15
to sy...@googlegroups.com
What should be the level of abstraction of implementation details in the gsoc application?. Is it good to include heavy implementation details (that I have thought of yet)?. Should I include the details of as many algorithms as I can?. If not then where should I discuss those details with community?.

Aaron Meurer

unread,
Mar 7, 2015, 4:23:22 PM3/7/15
to sy...@googlegroups.com
A suggestion is to show a "fake" session where you show what the
module will do (like a doctest for the final function showing inputs
and outputs). This shows that you at least are thinking about the API,
and that you understand the problem.

For something as complicated as CAD, I would also make sure to include
many details on the algorithm, to show that you understand it. But do
keep the implementation in mind. It's more common for people to make
GSoC proposals too theoretic than it is to make them not theoretic
enough.

Also see https://asmeurersympy.wordpress.com/2011/04/27/advice-for-future-prospective-gsoc-students/.

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/7c331d28-3dfc-4cb4-8f64-1c7d8d86f3ee%40googlegroups.com.

Luv Agarwal

unread,
Mar 16, 2015, 5:02:04 PM3/16/15
to sy...@googlegroups.com
I have made an initial draft at [this](https://github.com/sympy/sympy/wiki/GSoC-2015-Application-Luv-Agarwal:-Cylindrical-Algebraic-Decomposition). I will be adding more implementing details asap.

It would be great if I can get feedbacks.

Thanks

Aaron Meurer

unread,
Mar 16, 2015, 5:08:11 PM3/16/15
to sy...@googlegroups.com
This looks like a good start.

I would add in the pseudocode (maybe in a new section if it is too
much) in each part, what functions will be needed to do that step, and
if those functions are implemented or not (or if something not
implemented will need to be implemented). These should be specific
functions inside of SymPy. That way, it is very clear what is already
implemented and what needs to be implemented. For example, we already
have a root isolation algorithm in SymPy. Is it sufficient for what
you are need here?

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/4c37d75f-dc7c-4919-b08f-6c83f49f30a4%40googlegroups.com.

Aaron Meurer

unread,
Mar 16, 2015, 5:10:23 PM3/16/15
to sy...@googlegroups.com
Does the variable order to CAD matter? In Python, {x, y} is a set, so
{x, y} == {y, x}. Should cylindrically_decompose(x**2 + y**2 < 1, {x,
y}) really be cylindrically_decompose(x**2 + y**2 < 1, [x, y])?

Aaron Meurer

Luv Agarwal

unread,
Mar 20, 2015, 9:40:07 PM3/20/15
to sy...@googlegroups.com
I have added the implementation details of few functions in the draft.

How  can I implement refine_root (and dup_isolate_real_roots_sqf) for polynomials over algebraic numbers?

Kalevi Suominen

unread,
Mar 21, 2015, 2:45:34 AM3/21/15
to sy...@googlegroups.com


On Saturday, March 21, 2015 at 3:40:07 AM UTC+2, Luv Agarwal wrote:
I have added the implementation details of few functions in the draft.

How  can I implement refine_root (and dup_isolate_real_roots_sqf) for polynomials over algebraic numbers?

One would need a method for determining the sign of an algebraic number. Those in AlgebraicField (is_positive,
is_negative, ...) are not suitable for this purpose.
Message has been deleted

Luv Agarwal

unread,
Mar 22, 2015, 5:33:43 AM3/22/15
to sy...@googlegroups.com


On Sunday, March 22, 2015 at 2:49:37 PM UTC+5:30, Luv Agarwal wrote:
Thanks, I looked at them.
How can I implement dup_isolate_real_roots_sqf for algebraic domains?

Kalevi Suominen

unread,
Mar 22, 2015, 7:16:40 AM3/22/15
to sy...@googlegroups.com


On Sunday, March 22, 2015 at 11:33:43 AM UTC+2, Luv Agarwal wrote:
Thanks, I looked at them.
How can I implement dup_isolate_real_roots_sqf for algebraic domains?

It seems that the algorithm in rootisolation should work for real algebraic numbers as well.
The critical points are the inequalities and log() in  dup_root_upper_bound  and  is_negative()  in
dup_sign_variations.
Making them work for algebraic numbers is non-trivial. Otherwise the code for root isolation
could be used as such except for relaxing the checks for QQ or ZZ.

Luv Agarwal

unread,
Mar 23, 2015, 8:44:03 PM3/23/15
to sy...@googlegroups.com
Thanks a lot Kalevi. I have taken a look at all these and understand the requirements for root isolation in Algebraic domain.

I have updated my proposal.
Please have a look especially at the Timeline part. Is there anything that I should add or remove?.

Thanks
Luv Agarwal

Kalevi Suominen

unread,
Mar 24, 2015, 8:26:59 AM3/24/15
to sy...@googlegroups.com


On Tuesday, March 24, 2015 at 2:44:03 AM UTC+2, Luv Agarwal wrote

I have updated my proposal.
Please have a look especially at the Timeline part. Is there anything that I should add or remove?.

Thanks
Luv Agarwal

Thank you, Luv. I feel there is something I should add about root isolation. It may also affect the
timetable.

Implementing  log()  for algebraic numbers should not be difficult. One can pass to numerical values.
They are approximate, of course, but then  log  itself is only used for an estimate.

The comparisons of algebraic numbers are harder to implement. As noted before, the methods
presented in AlgebraicField cannot be used because they do not conform to the standard order of
real numbers. Instead, one could proceed as follows.
 
If  a  and  b  are algebraic numbers, the inequality  a > b  is equivalent to  a - b > 0. So it is enough to
find the sign of  a - b. In other words, to find a field  K  of algebraic numbers containing  a  and  b , and
to implement  K.sign(x).

An algebraic number field  K = Q(theta)  is generated by a single algebraic number  theta , called the
primitive element of the field. If the degree of the minimal polynomial of  theta  is  n , each element of
K  has a unique representation as a polynomial  a = p(theta)  with rational coefficients and of degree < n.

To compute the sign of  a , one looks for an interval  (u, v)  containing  theta  and such that the polynomial 
p(t)  has no roots in the interval. Then the sign of  a  is the sign of  p(t)  at any  t  in the interval.

One can start with an interval isolating the root  theta  of the minimal polynomial and then refine it until
it contains no roots of  p. Such an interval exists if  p != 0 , since then also  p(theta) != 0. One can do the
refinement by simply bisecting the interval, or else following the method of  dup_isolate_real_roots_sqf
based on representing the root by a continued fraction.

The classical method of determining the number of roots in an interval is by means of Sturm's theorem.
The method used in  dup_isolate_real_roots_sqf  is simpler in principle. It is based on using a
Moebius mapping to transform the interval to the positive real axis. If the coefficients of the transformed
polynomial have the same sign (as reported by dup_sign_variations), it has no positive roots.
(This is special case of Descartes' rule of signs which is trivial to verify.) Then the original polynomial
has no roots in the refined isolating interval.
An extra bonus is that the coefficients of the Moebius mapping are readily obtained from the convergents
of the continued fractions.

In summary: Implementation of inequalities in algebraic number fields is a substantial task. You may want
to reconsider your timeline on that part.

Best wishes,
Kalevi

Aaron Meurer

unread,
Mar 24, 2015, 5:28:11 PM3/24/15
to sy...@googlegroups.com
The current algebraic extension code in the polys (like
Poly(extension=[sqrt(x)])) is quite weak. Will this project require
expanding those abilities?

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/3bd8afe9-f2b9-49cf-a2a2-9381e6f0eb35%40googlegroups.com.

Kalevi Suominen

unread,
Mar 24, 2015, 5:54:23 PM3/24/15
to sy...@googlegroups.com


On Tuesday, March 24, 2015 at 11:28:11 PM UTC+2, Aaron Meurer wrote:
The current algebraic extension code in the polys (like
Poly(extension=[sqrt(x)])) is quite weak. Will this project require
expanding those abilities?

Aaron,
I have not delved deeper into the matter, but it looks like the current methods in
numberfields would essentially suffice, at least if they function as expected.
The central one is  primitive_element()  which enables the
construction of an extension containing any (finite) family of algebraic numbers.
On the other hand, the methods in AlgebraicField require revision.

Kalevi 

Aaron Meurer

unread,
Mar 24, 2015, 5:58:34 PM3/24/15
to sy...@googlegroups.com
There was also a GSoC project on this a couple of years ago
(https://github.com/sympy/sympy/wiki/GSoC-2013-Application-Katja-Sophie-Hotz:-Faster-Algorithms-for-Polynomials-over-Algebraic-Number-Fields),
and I am not fully caught up on what still needs to be implemented.
IIRC, though, there were some things that were technically
implemented, but were too slow to be usable in practice (particularly
multiple extensions, if I am not mistaken).

Aaron Meurer
> https://groups.google.com/d/msgid/sympy/c9fbe82d-800b-4fe8-8c84-54b646a07cd4%40googlegroups.com.

Luv Agarwal

unread,
Mar 24, 2015, 7:18:41 PM3/24/15
to sy...@googlegroups.com


On Wednesday, March 25, 2015 at 3:28:34 AM UTC+5:30, Aaron Meurer wrote:
There was also a GSoC project on this a couple of years ago
(https://github.com/sympy/sympy/wiki/GSoC-2013-Application-Katja-Sophie-Hotz:-Faster-Algorithms-for-Polynomials-over-Algebraic-Number-Fields),
and I am not fully caught up on what still needs to be implemented.
IIRC, though, there were some things that were technically
implemented, but were too slow to be usable in practice (particularly
multiple extensions, if I am not mistaken).
In my project this is how I deal with multiple extensions. Whenever multiple algebraic numbers need to be substituted in a multivariate polynomial then instead of substituting them directly, a single algebraic number (say theta) is calculated using simple (2 numbers at a time) such that all these algebraic numbers can be expressed as polynomials in theta and then substituted in multivariate polynomial all in terms of theta. So the new domain of substituted polynomial is a simple extension and the papers say that this conversion is intensive but Partial CAD (which is part of my project) provides some tricks to sometimes escape these calculations.

Thanks

Luv Agarwal

unread,
Mar 24, 2015, 7:54:24 PM3/24/15
to sy...@googlegroups.com


On Wednesday, March 25, 2015 at 3:24:23 AM UTC+5:30, Kalevi Suominen wrote:


On Tuesday, March 24, 2015 at 11:28:11 PM UTC+2, Aaron Meurer wrote:
The current algebraic extension code in the polys (like
Poly(extension=[sqrt(x)])) is quite weak. Will this project require
expanding those abilities?

Aaron,
I have not delved deeper into the matter, but it looks like the current methods in
numberfields would essentially suffice, at least if they function as expected.
The central one is  primitive_element()  which enables the
Is primitive_element doing something similar to simple?. If yes, will there be difference in efficiency of both of these as they take input in different form?.

Thanks

Kalevi Suominen

unread,
Mar 25, 2015, 8:41:44 AM3/25/15
to sy...@googlegroups.com


On Wednesday, March 25, 2015 at 1:54:24 AM UTC+2, Luv Agarwal wrote:

Is primitive_element doing something similar to simple?. If yes, will there be difference in efficiency of both of these as they take input in different form?.

primitive_element and simple do more or less the same thing. Both construct an extension containing given
algebraic elements. Their differences are

    1.  primitive_element takes any number of elements, simple only two.

    2.  simple computes isolating intervals, primitive_element does not care.

In more detail, primitive_element  uses  sqf_norm  which is just steps 1 and 2 of simple combined.
(1. norm is computed as a resultant, 2. t  is chosen to make norm square free.)

As to efficiency, I don't think there would be essential difference, since the algorithms are basically
the same. Setting  ex=True (experimental?)  in primitive_element leads to a parallel variant using
Groebner bases. It might be more efficient for a large number of elements.

Kalevi

Luv Agarwal

unread,
Mar 25, 2015, 6:27:13 PM3/25/15
to sy...@googlegroups.com


On Tuesday, March 24, 2015 at 5:56:59 PM UTC+5:30, Kalevi Suominen wrote:


On Tuesday, March 24, 2015 at 2:44:03 AM UTC+2, Luv Agarwal wrote
I have updated my proposal.
Please have a look especially at the Timeline part. Is there anything that I should add or remove?.

Thanks
Luv Agarwal

Thank you, Luv. I feel there is something I should add about root isolation. It may also affect the
timetable.

Implementing  log()  for algebraic numbers should not be difficult. One can pass to numerical values.
They are approximate, of course, but then  log  itself is only used for an estimate.

The comparisons of algebraic numbers are harder to implement. As noted before, the methods
presented in AlgebraicField cannot be used because they do not conform to the standard order of
real numbers. Instead, one could proceed as follows.
 
If  a  and  b  are algebraic numbers, the inequality  a > b  is equivalent to  a - b > 0. So it is enough to
find the sign of  a - b. In other words, to find a field  K  of algebraic numbers containing  a  and  b , and
to implement  K.sign(x).

An algebraic number field  K = Q(theta)  is generated by a single algebraic number  theta , called the
primitive element of the field. If the degree of the minimal polynomial of  theta  is  n , each element of
K  has a unique representation as a polynomial  a = p(theta)  with rational coefficients and of degree < n.
This can be done using K.convert(a) 

To compute the sign of  a , one looks for an interval  (u, v)  containing  theta  and such that the polynomial 
p(t)  has no roots in the interval. Then the sign of  a  is the sign of  p(t)  at any  t  in the interval.

One can start with an interval isolating the root  theta  of the minimal polynomial and then refine it until
it contains no roots of  p. Such an interval exists if  p != 0 , since then also  p(theta) != 0. One can do the
refinement by simply bisecting the interval, or else following the method of  dup_isolate_real_roots_sqf
based on representing the root by a continued fraction.

The classical method of determining the number of roots in an interval is by means of Sturm's theorem.
The method used in  dup_isolate_real_roots_sqf  is simpler in principle. It is based on using a
Moebius mapping to transform the interval to the positive real axis. If the coefficients of the transformed
polynomial have the same sign (as reported by dup_sign_variations), it has no positive roots.
(This is special case of Descartes' rule of signs which is trivial to verify.) Then the original polynomial
has no roots in the refined isolating interval.
An extra bonus is that the coefficients of the Moebius mapping are readily obtained from the convergents
of the continued fractions.
Thanks Kalevi for giving such a detailed description. 
In summary: Implementation of inequalities in algebraic number fields is a substantial task. You may want
to reconsider your timeline on that part.
The current status of root isolation makes it pretty easy to implement the above idea and I have given sufficient amount of time to this part of implementation in my timeline. What reconsideration are you talking about?
 
Best wishes,
Kalevi

Kalevi Suominen

unread,
Mar 26, 2015, 3:32:09 AM3/26/15
to sy...@googlegroups.com


On Thursday, March 26, 2015 at 12:27:13 AM UTC+2, Luv Agarwal wrote:
 
 What reconsideration are you talking about?

I was just pondering if you really have time to complete also  isolate_roots  in week 1 already.
Even  proj  might take some time... 

Kalevi

Luv Agarwal

unread,
Mar 26, 2015, 3:45:18 AM3/26/15
to sy...@googlegroups.com
Yeah it's big to complete in one week but it's not just one week but also community bonding period which is about 4 weeks long. My exams are ending at 30 April and I will be starting right after that.

Thanks
Kalevi

Luv Agarwal

unread,
Mar 26, 2015, 7:22:26 PM3/26/15
to sy...@googlegroups.com
I have updated my proposal at https://github.com/sympy/sympy/wiki/GSoC-2015-Application-Luv-Agarwal:-Cylindrical-Algebraic-Decomposition.
Please have a look. Am I allowed to edit my proposal at wiki after deadline?

Thanks

Aaron Meurer

unread,
Mar 26, 2015, 7:34:43 PM3/26/15
to sy...@googlegroups.com
You should finish your proposal by the deadline. If you have any
comments to make after that, you can comment on the proposal in
Melange, or here. We will be reading the proposals in Melange, and
likely won't see any updates on the wiki.

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/5132c7ec-eb2c-4326-8598-9777576d19e8%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages