On Sun, Aug 26, 2012 at 11:07 AM, Dima Pasechnik <dimp...@gmail.com> wrote: > On 2012-08-26, Geoffrey Irving <irv...@naml.us> wrote: >> On Sun, Aug 26, 2012 at 3:38 AM, Dima Pasechnik <dimp...@gmail.com> wrote: >>> On 2012-08-25, Geoffrey Irving <geoffrey.irv...@gmail.com> wrote: >>>> On Aug 24, 2012, at 6:19 PM, Geoffrey Irving <irv...@naml.us> wrote: >>>> >>>>> Hello, >>>>> >>>>> I'm implementing a code generator for simulation of simplicity (a way >>>>> to handle degeneracies in exact geometric computation). This >>>>> involves constructing and analyzing polynomials over a number of >>>>> coordinate variables plus one special infinitesimal variable 'e'. >>>>> For example, here's the polynomial for a 2x2 determinant (e.g., for >>>>> checking triangle areas). We add a different power of e to each >>>>> coordinate variable in an attempt to make sure our expressions are >>>>> never exactly zero: >>>>> >>>>> sage: R.<a,b,c,d,e> = >>>>> PolynomialRing(PolynomialRing(QQ,'a,b,c,d'),'e',sparse=True) sage: p >>>>> = (a+e**2**(0+1))*(d+e**2**(10+2))-(b+e**2**(10+1))*(c+e**2**(0+2)); >>>>> p e^4098 + a*e^4096 - e^2052 - c*e^2048 - b*e^4 + d*e^2 - b*c + a*d >>>>> >>>>> In order to avoid all degeneracies, this function must be nonzero in >>>>> the limit of small but nonzero e regardless of the values of the >>>>> other variables. To check this, I need to know whether the system of >>>>> polynomial equations defined by the coefficients of the distinct >>>>> powers of e is solvable over the other variables. What is the >>>>> easiest way to do that? Specifically? >>>>> >>>>> 1. How do I extract the coefficients of e as an ordered list? I >>>>> thought p.coefficients() would do it since I constructed the >>>>> polynomial ring nested, but p.coefficients() treats the ring as >>>>> flattened. Is that an artifact of the special assignment notation I >>>>> used to generate the ring? >>>> >>>> Yep, it works fine if I avoid the .< stuff. >>>> >>>>> 2. Once I get this list of polynomials, how I check whether it's >>>>> solvable over the reals? Ideally it will be unsolvable, but for >>>>> debugging purposes I want to know some of the solutions if any do >>>>> exist. >>>> >>>> Actually, it may turn out that for all practical cases one of the >>>> coefficients of e is a nonzero constant, in which case a solve is >>>> entirely unnecessary. >>>> >>>> I would still like to know how to handle the all nonconstant case, >>>> though. >>> >>> This is what semi-algebraic geometry is for. I don't think Sage >>> implements much of it. >>> How much do you know about the mathematics in question? >> >> Not much: I just looked up the relevant theorems now (Hilbert's >> Nullstellensatz and such). >> >>> There are basically two things one can do: one is to use a refinement of >>> the classical quantifier elimination over the reals approach, as >>> described in e.g. >>> http://perso.univ-rennes1.fr/marie-francoise.roy/bpr-ed2-posted1.html >>> (there is some Maxima and Singular code available that implements some >>> pieces of it...) >>> the other to use semidefinite programming based appoach, which boils down to >>> approximating nonnegative polynomials by sums of squares of polynomials: >>> http://books.google.com/books/about/Moments_Positive_Polynomials_and_Their_A.html?id=VY6imTsdIrEC&redir_esc=y >>> http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Tutorials.MomentRelaxations >> >> In the cases I've tried so far one of my equations is always a unit, >> so unsolvability is trivial. If I run across a more complicated case >> it is likely that unsolvability over C will be sufficient, so it looks >> like I'm unlikely to need the semi-algebraic case. > > unsolvability over C can be dealt with by Groebner bases. > This is available in Sage.
Yep, Volker's link above gave a good example of that. Thanks for the help! Geoffrey -- You received this message because you are subscribed to the Google Groups "sage-support" group. To post to this group, send email to sage-support@googlegroups.com. To unsubscribe from this group, send email to sage-support+unsubscr...@googlegroups.com. Visit this group at http://groups.google.com/group/sage-support?hl=en.
