[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-developer] [Cylindrical Algebraic Decomposition] (new)
From: |
Bill Page |
Subject: |
[Axiom-developer] [Cylindrical Algebraic Decomposition] (new) |
Date: |
Wed, 12 Oct 2005 16:49:34 -0500 |
Changes http://wiki.axiom-developer.org/CylindricalAlgebraicDecomposition/diff
--
**On Wednesday, October 12, 2005 5:29 PM Renaud Rioboo wrote:**
Dear Axiom fans,
In the last few years there have been a regain of interest in some Computer
Algebra users circles about Cylindrical Algebraic Decomposition (aka CAD). As
some of you may know my thesis was about CAD and I am regularly receiving
queries about CAD and the Axiom implementation I wrote.
I already privately gave sources or ran particular problems for some people
but I have never officially released my sources because I do not think they
are at production level and many work has still to be done on them.
Recently Martin Rubey thought it could be a good idea to release them and I
have of course no objection on that point nor on sharing these sources with
the community. If there is some interest on it you may include it into the
Axiom distribution with the same permissions than the rest of the software.
While I wait for my lab to give me the necessary permissions to be able to use
Axiom and export this software you may find it's sources at the unlisted url
http://rioboo.free.fr/CadPub/
I would like to add just a few words about that work which is a
straightforward implementation of the standard papers by Arnon, Collins and
McCallum which are more than 20 years old now.
My work on Cad is 15 years old and was developped for the defense of my thesis
where CAD was an sample application for real algebraic numbers
manipulations. By the time I wrote it I also needed some subresultant
calculations which were not in Axiom (Scratchpad by that time) and the
required machinery to work with real algebraic numbers.
Along the times, real algebraic numbers were included into Axiom and
subresultant calculations which are used in several places of the Axiom
Library were modified to all use Lionel Ducos's package which enables a
performance gain. While making these modifications to Axiom I tried to keep my
CAD package up-to-date with the Axiom Library and have used it as a test for
other packages. It thus should compile under recent Axiom versions and should
compute something that resembles a CAD.
Don't expect this package to be the absolute CAD program, I never found time
to further work on it in order to include many enhancements which are present
in other CAD packages. While Axiom versions evolved I had to remove what I
thought were enhancements and remove some support for generalization. As of
today it only includes one projection method which is the McCallum projection
and no reconstruction nor adjacency's are taken into account.
However it did reasonabily compare with more recent techniques such as
Rational Univariate Representation (aka RUR) for simple cases though it could
not produce results for some more complicated cases. By the time this
comparison was made Axiom had some severe memory restrictions which seem to
have disappear now. I thus think that many objections to using CAD could now
be revised.
For hard problems this package should thus not be worse than any other CAD
package except that you cannot fall into an optimzation drawback :-) There is
no optimization !
I will of course try to maintain the package to the best of my ability. One
thing I should do is include comments describing the different functions, but
even that requires going into the code and the algorithms which are quite old
for me now.
Let me know if there is some interest on it.
Best regards,
Renaud
<hr />
\begin{spad}
)ab pack CADU CylindricalAlgebraicDecompositionUtilities
CylindricalAlgebraicDecompositionUtilities(R,P) : PUB == PRIV where
-- Tese are some standard tools which are needed to compute with univariate
-- polynomials.
--
-- A gcd basis for a set of polynomials is a set of pairwise relatively prime
-- polynomials which all divide the original set and whose product is the
-- same than the product of the original set.
--
-- A square free basis for a list of polynomials is just a list
-- of square free polynomials which are pairwise relatively primes and have
-- the same roots than the original polynomials.
R : GcdDomain
P : UnivariatePolynomialCategory(R)
PUB == with
squareFreeBasis : List(P) -> List(P)
++
gcdBasis : List(P) -> List(P)
++ decompose a list of polynomials into pairwise relatively
++ prime polynomials
gcdBasisAdd : (P,List(P)) -> List(P)
++ add one polynomial to list of pairwise relatively prime polynomials
PRIV == add
squareFreeBasis(lpols) ==
lpols = [] => []
pol := first(lpols)
sqpol := unitCanonical(squareFreePart(pol))
gcdBasis(cons(sqpol,squareFreeBasis(rest(lpols))))
gcdBasisAdd(p,lpols) ==
(degree(p) = 0) => lpols
null lpols => [unitCanonical p]
p1 := first(lpols)
g := gcd(p,p1)
(degree(g) = 0) => cons(p1,gcdBasisAdd(p,rest lpols))
p := (p exquo g)::P
p1 := (p1 exquo g)::P
basis := gcdBasisAdd(p,rest(lpols))
if degree(p1) > 0 then basis := cons(p1,basis)
gcdBasisAdd(g,basis)
gcdBasis(lpols) ==
(#lpols <= 1) => lpols
basis := gcdBasis(rest lpols)
gcdBasisAdd(first(lpols),basis)
\end{spad}
\begin{spad}
)ab domain SCELL SimpleCell
-- This domain is made to work with one-dimensional semi-algebraic cells
-- ie either an algebraic point, either an interval. The point is specified as
-- specification of an algebraic value.
SimpleCell(TheField,ThePols) : PUB == PRIV where
TheField : RealClosedField
ThePols : UnivariatePolynomialCategory(TheField)
O ==> OutputForm
B ==> Boolean
Z ==> Integer
N ==> NonNegativeInteger
-- VARS ==> VariationsPackage(TheField,ThePols)
VARS ==> RealPolynomialUtilitiesPackage(TheField,ThePols)
LF ==> List(TheField)
PUB == CoercibleTo(O) with
allSimpleCells : (ThePols,Symbol) -> List %
allSimpleCells : (List(ThePols),Symbol) -> List %
hasDimension? : % -> B
samplePoint : % -> TheField
stablePol : % -> ThePols
variableOf : % -> Symbol
separe : (LF,TheField,TheField) -> LF
pointToCell : (TheField,B,Symbol) -> %
PRIV == add
Rep := Record(samplePoint:TheField,
hasDim:B,
varOf:Symbol)
samplePoint(c) == c.samplePoint
stablePol(c) == error "Prout"
hasDimension?(c) == c.hasDim
variableOf(c) == c.varOf
coerce(c:%):O ==
o : O := ((c.varOf)::O) = ((c.samplePoint)::O)
brace [o,(c.hasDim)::O]
separe(liste,gauche,droite) ==
milieu : TheField := (gauche + droite) / (2::TheField)
liste = [] => [milieu]
#liste = 1 => [gauche,first(liste),droite]
nbe := first(liste)
lg :List(TheField) := []
ld :List(TheField) := rest(liste)
sg := sign(milieu-nbe)
while sg > 0 repeat
lg := cons(nbe,lg)
ld = [] => return(separe(reverse(lg),gauche,milieu))
nbe := first(ld)
sg := sign(milieu-nbe)
ld := rest(ld)
sg < 0 =>
append(separe(reverse(lg),gauche,milieu),
rest(separe(cons(nbe,ld),milieu,droite)))
newDroite := (gauche+milieu)/(2::TheField)
null lg =>
newGauche := (milieu+droite)/(2::TheField)
while newGauche >= first(ld) repeat
newGauche := (milieu+newGauche)/(2::TheField)
append([gauche,milieu],separe(ld,newGauche,droite))
while newDroite <= first(lg) repeat
newDroite := (newDroite+milieu)/(2::TheField)
newGauche := (milieu+droite)/(2::TheField)
null ld => append(separe(reverse(lg),gauche,newDroite),[milieu,droite])
while newGauche >= first(ld) repeat
newGauche := (milieu+newGauche)/(2::TheField)
append(separe(reverse(lg),gauche,newDroite),
cons(milieu,separe(ld,newGauche,droite)))
pointToCell(sp,hasDim?,varName) ==
[sp,hasDim?,varName]$Rep
allSimpleCells(p:ThePols,var:Symbol) ==
allSimpleCells([p],var)
PACK ==> CylindricalAlgebraicDecompositionUtilities(TheField,ThePols)
allSimpleCells(lp:List(ThePols),var:Symbol) ==
lp1 := gcdBasis(lp)$PACK
null(lp1) => [pointToCell(0,true,var)]
b := ("max" / [ boundOfCauchy(p)$VARS for p in lp1 ])::TheField
l := "append" / [allRootsOf(makeSUP(unitCanonical(p))) for p in lp1]
l := sort(l)
l1 := separe(l,-b,b)
res : List(%) := [pointToCell(first(l1),true,var)]
l1 := rest(l1)
while not(null(l1)) repeat
res := cons(pointToCell(first(l1),false,var),res)
l1 := rest(l1)
l1 = [] => return(error "Liste vide")
res := cons(pointToCell(first(l1),true,var),res)
l1 := rest(l1)
reverse! res
\end{spad}
\begin{spad}
)ab domain CELL Cell
Cell(TheField) : PUB == PRIV where
TheField : RealClosedField
ThePols ==> Polynomial(TheField)
O ==> OutputForm
B ==> Boolean
Z ==> Integer
N ==> NonNegativeInteger
BUP ==> SparseUnivariatePolynomial(TheField)
SCELL ==> SimpleCell(TheField,BUP)
PUB == CoercibleTo(O) with
samplePoint : % -> List(TheField)
dimension : % -> N
hasDimension? : (%,Symbol) -> B
makeCell : List(SCELL) -> %
makeCell : (SCELL,%) -> %
mainVariableOf : % -> Symbol
variablesOf : % -> List Symbol
projection : % -> Union(%,"failed")
PRIV == add
Rep := List(SCELL)
coerce(c:%):O ==
paren [sc::O for sc in c]
projection(cell) ==
null cell => error "projection: should not appear"
r := rest(cell)
null r => "failed"
r
makeCell(l:List(SCELL)) == l
makeCell(scell,toAdd) == cons(scell,toAdd)
mainVariableOf(cell) ==
null(cell) =>
error "Should not appear"
variableOf(first(cell))
variablesOf(cell) ==
null(cell) => []
cons(mainVariableOf(cell),variablesOf(rest(cell)::%))
dimension(cell) ==
null(cell) => 0
hasDimension?(first(cell)) => 1+dimension(rest(cell))
dimension(rest(cell))
hasDimension?(cell,var) ==
null(cell) =>
error "Should not appear"
sc : SCELL := first(cell)
v := variableOf(sc)
v = var => hasDimension?(sc)
v < var => false
v > var => true
error "Caca Prout"
samplePoint(cell) ==
null(cell) => []
cons(samplePoint(first(cell)),samplePoint(rest(cell)))
\end{spad}
\begin{spad}
)ab pack CAD CylindricalAlgebraicDecompositionPackage
CylindricalAlgebraicDecompositionPackage(TheField) : PUB == PRIV where
TheField : RealClosedField
ThePols ==> Polynomial(TheField)
P ==> ThePols
BUP ==> SparseUnivariatePolynomial(TheField)
RUP ==> SparseUnivariatePolynomial(ThePols)
Z ==> Integer
N ==> NonNegativeInteger
CELL ==> Cell(TheField)
SCELL ==> SimpleCell(TheField,BUP)
PUB == with
cylindricalDecomposition: List P -> List CELL
cylindricalDecomposition: (List(P),List(Symbol)) -> List CELL
projectionSet: (List RUP) -> List P
coefficientSet: RUP -> List P
discriminantSet : List RUP -> List(P)
resultantSet : List RUP -> List P
principalSubResultantSet : (RUP,RUP) -> List P
specialise : (List(ThePols),CELL) -> List(BUP)
PRIV == add
cylindricalDecomposition(lpols) ==
lv : List(Symbol) := []
for pol in lpols repeat
ground?(pol) => "next pol"
lv := removeDuplicates(append(variables(pol),lv))
lv := reverse(sort(lv))
cylindricalDecomposition(lpols,lv)
cylindricalDecomposition(lpols,lvars) ==
lvars = [] => error("CAD: cylindricalDecomposition: empty list of vars")
mv := first(lvars)
lv := rest(lvars)
lv = [] =>
lp1 := [ univariate(pol) for pol in lpols ]
scells := allSimpleCells(lp1,mv)$SCELL
[ makeCell([scell]) for scell in scells ]
lpols1 := projectionSet [univariate(pol,mv) for pol in lpols]
previousCad := cylindricalDecomposition(lpols1,lv)
res : List(CELL) := []
for cell in previousCad repeat
lspec := specialise(lpols,cell)
scells := allSimpleCells(lspec,mv)
res := append(res,[makeCell(scell,cell) for scell in scells])
res
PACK1 ==> CylindricalAlgebraicDecompositionUtilities(ThePols,RUP)
PACK2 ==> CylindricalAlgebraicDecompositionUtilities(TheField,BUP)
specialise(lpols,cell) ==
lpols = [] => error("CAD: specialise: empty list of pols")
sp := samplePoint(cell)
vl := variablesOf(cell)
res : List(BUP) := []
for pol in lpols repeat
p1 := univariate(eval(pol,vl,sp))
-- zero?(p1) => return(error "Bad sepcialise")
degree(p1) = 0 => "next pol"
res := cons(p1,res)
res
coefficientSet(pol) ==
res : List(ThePols) := []
for c in coefficients(pol) repeat
ground?(c) => return(res)
res := cons(c,res)
res
SUBRES ==> SubResultantPackage(ThePols,RUP)
discriminantSet(lpols) ==
res : List(ThePols) := []
for p in lpols repeat
v := subresultantVector(p,differentiate(p))$SUBRES
not(zero?(degree(v.0))) => return(error "Bad discriminant")
d : ThePols := leadingCoefficient(v.0)
-- d := discriminant p
zero?(d) => return(error "Non Square Free polynomial")
if not(ground? d) then res := cons(d,res)
res
principalSubResultantSet(p,q) ==
if degree(p) < degree(q)
then (p,q) := (q,p)
if degree(p) = degree(q)
then (p,q) := (q,pseudoRemainder(p, q))
v := subresultantVector(p,q)$SUBRES
[coefficient(v.i,i) for i in 0..(((#v)-2)::N)]
resultantSet(lpols) ==
res : List(ThePols) := []
laux := lpols
for p in lpols repeat
laux := rest(laux)
for q in laux repeat
r : ThePols := first(principalSubResultantSet(p,q))
-- r := resultant(p,q)
zero?(r) => return(error "Non relatively prime polynomials")
if not(ground? r) then res := cons(r,res)
res
projectionSet(lpols) ==
res : List(ThePols) := []
for p in lpols repeat
c := content(p)
ground?(c) => "next p"
res := cons(c,res)
lp1 := [primitivePart p for p in lpols]
f : ((RUP,RUP) -> Boolean) := (degree(#1) <= degree(#2))
lp1 := sort(f,lp1)
lsqfrb := squareFreeBasis(lp1)$PACK1
lsqfrb := sort(f,lsqfrb)
for p in lp1 repeat
res := append(res,coefficientSet(p))
res := append(res,discriminantSet(lsqfrb))
append(res,resultantSet(lsqfrb))
\end{spad}
\begin{axiom}
)lib )dir .
)lib SCELL CELL CAD
Ran := RECLOS(FRAC INT)
Cad := CAD(Ran)
p1 : POLY(Ran) := x^2+y^2-1
p2 : POLY(Ran) := y-x^2
lpols : List(POLY(Ran)) := [p1,p2]
cad := cylindricalDecomposition(lpols)$Cad
precision(30)
ls := [ samplePoint cell for cell in cad]
lf := [ [relativeApprox(coord,2^(-precision()))::Float for coord in point] for
point in ls]
lp := [ point(term::List(DoubleFloat))$Point(DoubleFloat) for term in lf ]
[ sign(eval(p1,variables(p1),samplePoint(cell))::Ran) for cell in cad ]
--[ sign(eval(p2,variables(p2),samplePoint(cell))::Ran) for cell in cad ]
\end{axiom}
--
forwarded from http://wiki.axiom-developer.org/address@hidden
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Axiom-developer] [Cylindrical Algebraic Decomposition] (new),
Bill Page <=