aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/reclos.spad.pamphlet
diff options
context:
space:
mode:
authordos-reis <gdr@axiomatics.org>2007-08-14 05:14:52 +0000
committerdos-reis <gdr@axiomatics.org>2007-08-14 05:14:52 +0000
commitab8cc85adde879fb963c94d15675783f2cf4b183 (patch)
treec202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/reclos.spad.pamphlet
downloadopen-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz
Initial population.
Diffstat (limited to 'src/algebra/reclos.spad.pamphlet')
-rw-r--r--src/algebra/reclos.spad.pamphlet1242
1 files changed, 1242 insertions, 0 deletions
diff --git a/src/algebra/reclos.spad.pamphlet b/src/algebra/reclos.spad.pamphlet
new file mode 100644
index 00000000..b0400ad1
--- /dev/null
+++ b/src/algebra/reclos.spad.pamphlet
@@ -0,0 +1,1242 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra reclos.spad}
+\author{Renaud Rioboo}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+This file describes the Real Closure 1.0 package which consists of different
+packages, categoris and domains :
+
+- the package RealPolynomialUtilitiesPackage whichs receives a field and a
+univariate polynomial domain with coefficients in the field. It computes some
+simple functions such as Strum and Sylvester sequences.
+
+- The category RealRootCharacterizationCategory provides abstarct
+functionalities to work with "real roots" of univariate polynomials. These
+resemble variables with some functionalities needed to compute important
+operations.
+
+- RealClosedField is a category with provides comon operations available over
+real closed fiels. These include finding all the roots of univariate
+polynomial, taking square roots, ...
+
+- The domain RightOpenIntervalRootCharacterization is the main code that
+provides the functionalities of RealRootCharacterizationCategory for the case
+of archimedean fileds. Abstract roots are encoded with a left closed right
+open interval containing the root together with a defining polynomial for the
+root.
+
+- The RealClosure domain is the end-user code, it provides usual arithmetics
+with real algebraic numbers, along with the functionalities of a real closed
+field. It also provides functions to approximate a real algebraic number by an
+element of the base field. This approximation may either be absolute
+(approximate) or relative (realtivApprox).
+
+
+CAVEEATS
+
+Since real algebraic expressions are stored as depending on "real roots" which
+are managed like variables, there is an ordering on these. This ordering is
+dynamical in the sense that any new algebraic takes precedence over older
+ones. In particular every cretaion function raises a new "real root". This has
+the effect that when you type something like sqrt(2) + sqrt(2) you have two
+new variables which happen to be equal. To avoid this name the expression such
+as in s2 := sqrt(2) ; s2 + s2
+
+Also note that computing times depend strongly on the ordering you implicitly
+provide. Please provide algebraics in the order which most natural to you.
+
+LIMITATIONS
+
+The file reclos.input show some basic use of the package. This packages uses
+algorithms which are published in [1] and [2] which are based on field
+arithmetics, inparticular for polynomial gcd related algorithms. This can be
+quite slow for high degree polynomials and subresultants methods usually work
+best. Betas versions of the package try to use these techniques in a better
+way and work significantly faster. These are mostly based on unpublished
+algorithms and cannot be distributed. Please contact the author if you have a
+particular problem to solve or want to use these versions.
+
+Be aware that approximations behave as post-processing and that all
+computations are done excatly. They can thus be quite time consuming when
+depending on several "real roots".
+\section{package POLUTIL RealPolynomialUtilitiesPackage}
+<<package POLUTIL RealPolynomialUtilitiesPackage>>=
+)abbrev package POLUTIL RealPolynomialUtilitiesPackage
+++ Author: Renaud Rioboo
+++ Date Created: summer 1992
+++ Basic Functions: provides polynomial utilities
+++ Related Constructors: RealClosure,
+++ Date Last Updated: July 2004
+++ Also See:
+++ AMS Classifications:
+++ Keywords: Sturm sequences
+++ References:
+++ Description:
+++ \axiomType{RealPolynomialUtilitiesPackage} provides common functions used
+++ by interval coding.
+RealPolynomialUtilitiesPackage(TheField,ThePols) : PUB == PRIV where
+
+ TheField : Field
+ ThePols : UnivariatePolynomialCategory(TheField)
+
+ Z ==> Integer
+ N ==> NonNegativeInteger
+ P ==> ThePols
+
+ PUB == with
+
+ sylvesterSequence : (ThePols,ThePols) -> List ThePols
+ ++ \axiom{sylvesterSequence(p,q)} is the negated remainder sequence
+ ++ of p and q divided by the last computed term
+ sturmSequence : ThePols -> List ThePols
+ ++ \axiom{sturmSequence(p) = sylvesterSequence(p,p')}
+ if TheField has OrderedRing then
+ boundOfCauchy : ThePols -> TheField
+ ++ \axiom{boundOfCauchy(p)} bounds the roots of p
+ sturmVariationsOf : List TheField -> N
+ ++ \axiom{sturmVariationsOf(l)} is the number of sign variations
+ ++ in the list of numbers l,
+ ++ note that the first term counts as a sign
+ lazyVariations : (List(TheField), Z, Z) -> N
+ ++ \axiom{lazyVariations(l,s1,sn)} is the number of sign variations
+ ++ in the list of non null numbers [s1::l]@sn,
+
+
+ PRIV == add
+
+ sturmSequence(p) ==
+ sylvesterSequence(p,differentiate(p))
+
+ sylvesterSequence(p1,p2) ==
+ res : List(ThePols) := [p1]
+ while (p2 ^= 0) repeat
+ res := cons(p2 , res)
+ (p1 , p2) := (p2 , -(p1 rem p2))
+ if degree(p1) > 0
+ then
+ p1 := unitCanonical(p1)
+ res := [ term quo p1 for term in res ]
+ reverse! res
+
+ if TheField has OrderedRing
+ then
+
+ boundOfCauchy(p) ==
+ c :TheField := inv(leadingCoefficient(p))
+ l := [ c*term for term in rest(coefficients(p))]
+ null(l) => 1
+ 1 + ("max" / [ abs(t) for t in l ])
+
+-- sturmVariationsOf(l) ==
+-- res : N := 0
+-- lsg := sign(first(l))
+-- for term in l repeat
+-- if ^( (sg := sign(term) ) = 0 ) then
+-- if (sg ^= lsg) then res := res + 1
+-- lsg := sg
+-- res
+
+ sturmVariationsOf(l) ==
+ null(l) => error "POLUTIL: sturmVariationsOf: empty list !"
+ l1 := first(l)
+ -- first 0 counts as a sign
+ ll : List(TheField) := []
+ for term in rest(l) repeat
+ -- zeros don't count
+ if not(zero?(term)) then ll := cons(term,ll)
+ -- if l1 is not zero then ll = reverse(l)
+ null(ll) => error "POLUTIL: sturmVariationsOf: Bad sequence"
+ ln := first(ll)
+ ll := reverse(rest(ll))
+ -- if l1 is not zero then first(l) = first(ll)
+ -- if l1 is zero then first zero should count as a sign
+ zero?(l1) => 1 + lazyVariations(rest(ll),sign(first(ll)),sign(ln))
+ lazyVariations(ll, sign(l1), sign(ln))
+
+ lazyVariations(l,sl,sh) ==
+ zero?(sl) or zero?(sh) => error "POLUTIL: lazyVariations: zero sign!"
+ null(l) =>
+ if sl = sh then 0 else 1
+ null(rest(l)) =>
+ if zero?(first(l))
+ then error "POLUTIL: lazyVariations: zero sign!"
+ else
+ if sl = sh
+ then
+ if (sl = sign(first(l)))
+ then 0
+ else 2
+ -- in this case we save one test
+ else 1
+ s := sign(l.2)
+ lazyVariations([first(l)],sl,s) +
+ lazyVariations(rest(rest(l)),s,sh)
+
+@
+\section{category RRCC RealRootCharacterizationCategory}
+<<category RRCC RealRootCharacterizationCategory>>=
+)abbrev category RRCC RealRootCharacterizationCategory
+++ Author: Renaud Rioboo
+++ Date Created: summer 1992
+++ Date Last Updated: January 2004
+++ Basic Functions: provides operations with generic real roots of
+++ polynomials
+++ Related Constructors: RealClosure, RightOpenIntervalRootCharacterization
+++ Also See:
+++ AMS Classifications:
+++ Keywords: Real Algebraic Numbers
+++ References:
+++ Description:
+++ \axiomType{RealRootCharacterizationCategory} provides common acces
+++ functions for all real root codings.
+RealRootCharacterizationCategory(TheField, ThePols ) : Category == PUB where
+
+ TheField : Join(OrderedRing, Field)
+ ThePols : UnivariatePolynomialCategory(TheField)
+
+ Z ==> Integer
+ N ==> PositiveInteger
+
+ PUB ==>
+ SetCategory with
+
+ sign: ( ThePols, $ ) -> Z
+ ++ \axiom{sign(pol,aRoot)} gives the sign of \axiom{pol}
+ ++ interpreted as \axiom{aRoot}
+ zero? : ( ThePols, $ ) -> Boolean
+ ++ \axiom{zero?(pol,aRoot)} answers if \axiom{pol}
+ ++ interpreted as \axiom{aRoot} is \axiom{0}
+ negative?: ( ThePols, $ ) -> Boolean
+ ++ \axiom{negative?(pol,aRoot)} answers if \axiom{pol}
+ ++ interpreted as \axiom{aRoot} is negative
+ positive?: ( ThePols, $ ) -> Boolean
+ ++ \axiom{positive?(pol,aRoot)} answers if \axiom{pol}
+ ++ interpreted as \axiom{aRoot} is positive
+ recip: ( ThePols, $ ) -> Union(ThePols,"failed")
+ ++ \axiom{recip(pol,aRoot)} tries to inverse \axiom{pol}
+ ++ interpreted as \axiom{aRoot}
+ definingPolynomial: $ -> ThePols
+ ++ \axiom{definingPolynomial(aRoot)} gives a polynomial
+ ++ such that \axiom{definingPolynomial(aRoot).aRoot = 0}
+ allRootsOf: ThePols -> List $
+ ++ \axiom{allRootsOf(pol)} creates all the roots of \axiom{pol}
+ ++ in the Real Closure, assumed in order.
+ rootOf: ( ThePols, N ) -> Union($,"failed")
+ ++ \axiom{rootOf(pol,n)} gives the nth root for the order of the
+ ++ Real Closure
+ approximate : (ThePols,$,TheField) -> TheField
+ ++ \axiom{approximate(term,root,prec)} gives an approximation
+ ++ of \axiom{term} over \axiom{root} with precision \axiom{prec}
+
+ relativeApprox : (ThePols,$,TheField) -> TheField
+ ++ \axiom{approximate(term,root,prec)} gives an approximation
+ ++ of \axiom{term} over \axiom{root} with precision \axiom{prec}
+
+ add
+
+ zero?(toTest, rootChar) ==
+ sign(toTest, rootChar) = 0
+
+ negative?(toTest, rootChar) ==
+ sign(toTest, rootChar) < 0
+
+ positive?(toTest, rootChar) ==
+ sign(toTest, rootChar) > 0
+
+ rootOf(pol,n) ==
+ liste:List($):= allRootsOf(pol)
+ # liste > n => "failed"
+ liste.n
+
+ recip(toInv,rootChar) ==
+ degree(toInv) = 0 =>
+ res := recip(leadingCoefficient(toInv))
+ if (res case "failed") then "failed" else (res::TheField::ThePols)
+ defPol := definingPolynomial(rootChar)
+ d := principalIdeal([defPol,toInv])
+ zero?(d.generator,rootChar) => "failed"
+ if (degree(d.generator) ^= 0 )
+ then
+ defPol := (defPol exquo (d.generator))::ThePols
+ d := principalIdeal([defPol,toInv])
+ d.coef.2
+
+@
+\section{category RCFIELD RealClosedField}
+<<category RCFIELD RealClosedField>>=
+)abbrev category RCFIELD RealClosedField
+++ Author: Renaud Rioboo
+++ Date Created: may 1993
+++ Date Last Updated: January 2004
+++ Basic Functions: provides computations with generic real roots of
+++ polynomials
+++ Related Constructors: SimpleOrderedAlgebraicExtension, RealClosure
+++ Also See:
+++ AMS Classifications:
+++ Keywords: Real Algebraic Numbers
+++ References:
+++ Description:
+++ \axiomType{RealClosedField} provides common acces
+++ functions for all real closed fields.
+RealClosedField : Category == PUB where
+
+ E ==> OutputForm
+ SUP ==> SparseUnivariatePolynomial
+ OFIELD ==> Join(OrderedRing,Field)
+ PME ==> SUP($)
+ N ==> NonNegativeInteger
+ PI ==> PositiveInteger
+ RN ==> Fraction(Integer)
+ Z ==> Integer
+ POLY ==> Polynomial
+ PACK ==> SparseUnivariatePolynomialFunctions2
+
+ PUB == Join(CharacteristicZero,
+ OrderedRing,
+ CommutativeRing,
+ Field,
+ FullyRetractableTo(Fraction(Integer)),
+ Algebra Integer,
+ Algebra(Fraction(Integer)),
+ RadicalCategory) with
+
+ mainForm : $ -> Union(E,"failed")
+ ++ \axiom{mainForm(x)} is the main algebraic quantity name of
+ ++ \axiom{x}
+
+ mainDefiningPolynomial : $ -> Union(PME,"failed")
+ ++ \axiom{mainDefiningPolynomial(x)} is the defining
+ ++ polynomial for the main algebraic quantity of \axiom{x}
+
+ mainValue : $ -> Union(PME,"failed")
+ ++ \axiom{mainValue(x)} is the expression of \axiom{x} in terms
+ ++ of \axiom{SparseUnivariatePolynomial($)}
+
+ rootOf: (PME,PI,E) -> Union($,"failed")
+ ++ \axiom{rootOf(pol,n,name)} creates the nth root for the order
+ ++ of \axiom{pol} and names it \axiom{name}
+
+ rootOf: (PME,PI) -> Union($,"failed")
+ ++ \axiom{rootOf(pol,n)} creates the nth root for the order
+ ++ of \axiom{pol} and gives it unique name
+
+ allRootsOf: PME -> List $
+ ++ \axiom{allRootsOf(pol)} creates all the roots
+ ++ of \axiom{pol} naming each uniquely
+
+ allRootsOf: (SUP(RN)) -> List $
+ ++ \axiom{allRootsOf(pol)} creates all the roots
+ ++ of \axiom{pol} naming each uniquely
+
+ allRootsOf: (SUP(Z)) -> List $
+ ++ \axiom{allRootsOf(pol)} creates all the roots
+ ++ of \axiom{pol} naming each uniquely
+
+ allRootsOf: (POLY($)) -> List $
+ ++ \axiom{allRootsOf(pol)} creates all the roots
+ ++ of \axiom{pol} naming each uniquely
+
+ allRootsOf: (POLY(RN)) -> List $
+ ++ \axiom{allRootsOf(pol)} creates all the roots
+ ++ of \axiom{pol} naming each uniquely
+
+ allRootsOf: (POLY(Z)) -> List $
+ ++ \axiom{allRootsOf(pol)} creates all the roots
+ ++ of \axiom{pol} naming each uniquely
+
+ sqrt: ($,N) -> $
+ ++ \axiom{sqrt(x,n)} is \axiom{x ** (1/n)}
+
+ sqrt: $ -> $
+ ++ \axiom{sqrt(x)} is \axiom{x ** (1/2)}
+
+ sqrt: RN -> $
+ ++ \axiom{sqrt(x)} is \axiom{x ** (1/2)}
+
+ sqrt: Z -> $
+ ++ \axiom{sqrt(x)} is \axiom{x ** (1/2)}
+
+ rename! : ($,E) -> $
+ ++ \axiom{rename!(x,name)} changes the way \axiom{x} is printed
+
+ rename : ($,E) -> $
+ ++ \axiom{rename(x,name)} gives a new number that prints as name
+
+ approximate: ($,$) -> RN
+ ++ \axiom{approximate(n,p)} gives an approximation of \axiom{n}
+ ++ that has precision \axiom{p}
+
+ add
+
+ sqrt(a:$):$ == sqrt(a,2)
+
+ sqrt(a:RN):$ == sqrt(a::$,2)
+
+ sqrt(a:Z):$ == sqrt(a::$,2)
+
+ characteristic() == 0
+
+ rootOf(pol,n,o) ==
+ r := rootOf(pol,n)
+ r case "failed" => "failed"
+ rename!(r,o)
+
+ rootOf(pol,n) ==
+ liste:List($):= allRootsOf(pol)
+ # liste > n => "failed"
+ liste.n
+
+
+ sqrt(x,n) ==
+ n = 0 => 1
+ n = 1 => x
+ zero?(x) => 0
+ one?(x) => 1
+ if odd?(n)
+ then
+ r := rootOf(monomial(1,n) - (x :: PME), 1)
+ else
+ r := rootOf(monomial(1,n) - (x :: PME), 2)
+ r case "failed" => error "no roots"
+ n = 2 => rename(r,root(x::E)$E)
+ rename(r,root(x :: E, n :: E)$E)
+
+ (x : $) ** (rn : RN) == sqrt(x**numer(rn),denom(rn)::N)
+
+ nthRoot(x, n) ==
+ zero?(n) => x
+ negative?(n) => inv(sqrt(x,(-n) :: N))
+ sqrt(x,n :: N)
+
+ allRootsOf(p:SUP(RN)) == allRootsOf(map(#1 :: $ ,p)$PACK(RN,$))
+
+ allRootsOf(p:SUP(Z)) == allRootsOf(map(#1 :: $ ,p)$PACK(Z,$))
+
+ allRootsOf(p:POLY($)) == allRootsOf(univariate(p))
+
+ allRootsOf(p:POLY(RN)) == allRootsOf(univariate(p))
+
+ allRootsOf(p:POLY(Z)) == allRootsOf(univariate(p))
+
+@
+\section{domain ROIRC RightOpenIntervalRootCharacterization}
+\subsection{makeChar performance problem}
+The following lines of code, which check for a possible error,
+cause major performance problems and were removed by Renaud Rioboo,
+the original author. They were originally inserted for debugging.
+\begin{verbatim}
+ right <= left => error "ROIRC: makeChar: Bad interval"
+ (pol.left * pol.right) > 0 => error "ROIRC: makeChar: Bad pol"
+\end{verbatim}
+<<performance problem>>=
+@
+<<domain ROIRC RightOpenIntervalRootCharacterization>>=
+)abbrev domain ROIRC RightOpenIntervalRootCharacterization
+++ Author: Renaud Rioboo
+++ Date Created: summer 1992
+++ Date Last Updated: January 2004
+++ Basic Functions: provides computations with real roots of olynomials
+++ Related Constructors: RealRootCharacterizationCategory, RealClosure
+++ Also See:
+++ AMS Classifications:
+++ Keywords: Real Algebraic Numbers
+++ References:
+++ Description:
+++ \axiomType{RightOpenIntervalRootCharacterization} provides work with
+++ interval root coding.
+RightOpenIntervalRootCharacterization(TheField,ThePolDom) : PUB == PRIV where
+
+ TheField : Join(OrderedRing,Field)
+ ThePolDom : UnivariatePolynomialCategory(TheField)
+
+
+ Z ==> Integer
+ P ==> ThePolDom
+ N ==> NonNegativeInteger
+ B ==> Boolean
+ UTIL ==> RealPolynomialUtilitiesPackage(TheField,ThePolDom)
+ RRCC ==> RealRootCharacterizationCategory
+ O ==> OutputForm
+ TwoPoints ==> Record(low:TheField , high:TheField)
+
+ PUB == RealRootCharacterizationCategory(TheField, ThePolDom) with
+
+ left : $ -> TheField
+ ++ \axiom{left(rootChar)} is the left bound of the isolating
+ ++ interval
+ right : $ -> TheField
+ ++ \axiom{right(rootChar)} is the right bound of the isolating
+ ++ interval
+ size : $ -> TheField
+ ++ The size of the isolating interval
+ middle : $ -> TheField
+ ++ \axiom{middle(rootChar)} is the middle of the isolating
+ ++ interval
+ refine : $ -> $
+ ++ \axiom{refine(rootChar)} shrinks isolating interval around
+ ++ \axiom{rootChar}
+ mightHaveRoots : (P,$) -> B
+ ++ \axiom{mightHaveRoots(p,r)} is false if \axiom{p.r} is not 0
+ relativeApprox : (P,$,TheField) -> TheField
+ ++ \axiom{relativeApprox(exp,c,p) = a} is relatively close to exp
+ ++ as a polynomial in c ip to precision p
+
+ PRIV == add
+
+
+
+-- local functions
+
+
+ makeChar: (TheField,TheField,ThePolDom) -> $
+ refine! : $ -> $
+ sturmIsolate : (List(P), TheField, TheField,N,N) -> List TwoPoints
+ isolate : List(P) -> List TwoPoints
+ rootBound : P -> TheField
+-- varStar : P -> N
+ linearRecip : ( P , $) -> Union(P, "failed")
+ linearZero? : (TheField,$) -> B
+ linearSign : (P,$) -> Z
+ sturmNthRoot : (List(P), TheField, TheField,N,N,N) -> Union(TwoPoints,"failed")
+ addOne : P -> P
+ minus : P -> P
+ translate : (P,TheField) -> P
+ dilate : (P,TheField) -> P
+ invert : P -> P
+ evalOne : P -> TheField
+ hasVarsl: List(TheField) -> B
+ hasVars: P -> B
+
+-- Representation
+
+ Rep:= Record(low:TheField,high:TheField,defPol:ThePolDom)
+
+-- and now the code !
+
+
+ size(rootCode) ==
+ rootCode.high - rootCode.low
+
+ relativeApprox(pval,rootCode,prec) ==
+ -- beurk !
+ dPol := rootCode.defPol
+ degree(dPol) = 1 =>
+ c := -coefficient(dPol,0)/leadingCoefficient(dPol)
+ pval.c
+ pval := pval rem dPol
+ degree(pval) = 0 => leadingCoefficient(pval)
+ zero?(pval,rootCode) => 0
+ while mightHaveRoots(pval,rootCode) repeat
+ rootCode := refine(rootCode)
+ dpval := differentiate(pval)
+ degree(dpval) = 0 =>
+ l := left(rootCode)
+ r := right(rootCode)
+ a := pval.l
+ b := pval.r
+ while ( abs(2*(a-b)/(a+b)) > prec ) repeat
+ rootCode := refine(rootCode)
+ l := left(rootCode)
+ r := right(rootCode)
+ a := pval.l
+ b := pval.r
+ (a+b)/(2::TheField)
+ zero?(dpval,rootCode) =>
+ relativeApprox(pval,
+ [left(rootCode),
+ right(rootCode),
+ gcd(dpval,rootCode.defPol)]$Rep,
+ prec)
+ while mightHaveRoots(dpval,rootCode) repeat
+ rootCode := refine(rootCode)
+ l := left(rootCode)
+ r := right(rootCode)
+ a := pval.l
+ b := pval.r
+ while ( abs(2*(a-b)/(a+b)) > prec ) repeat
+ rootCode := refine(rootCode)
+ l := left(rootCode)
+ r := right(rootCode)
+ a := pval.l
+ b := pval.r
+ (a+b)/(2::TheField)
+
+ approximate(pval,rootCode,prec) ==
+ -- glurp
+ dPol := rootCode.defPol
+ degree(dPol) = 1 =>
+ c := -coefficient(dPol,0)/leadingCoefficient(dPol)
+ pval.c
+ pval := pval rem dPol
+ degree(pval) = 0 => leadingCoefficient(pval)
+ dpval := differentiate(pval)
+ degree(dpval) = 0 =>
+ l := left(rootCode)
+ r := right(rootCode)
+ while ( abs((a := pval.l) - (b := pval.r)) > prec ) repeat
+ rootCode := refine(rootCode)
+ l := left(rootCode)
+ r := right(rootCode)
+ (a+b)/(2::TheField)
+ zero?(dpval,rootCode) =>
+ approximate(pval,
+ [left(rootCode),
+ right(rootCode),
+ gcd(dpval,rootCode.defPol)]$Rep,
+ prec)
+ while mightHaveRoots(dpval,rootCode) repeat
+ rootCode := refine(rootCode)
+ l := left(rootCode)
+ r := right(rootCode)
+ while ( abs((a := pval.l) - (b := pval.r)) > prec ) repeat
+ rootCode := refine(rootCode)
+ l := left(rootCode)
+ r := right(rootCode)
+ (a+b)/(2::TheField)
+
+
+ addOne(p) == p.(monomial(1,1)+(1::P))
+
+ minus(p) == p.(monomial(-1,1))
+
+ translate(p,a) == p.(monomial(1,1)+(a::P))
+
+ dilate(p,a) == p.(monomial(a,1))
+
+ evalOne(p) == "+" / coefficients(p)
+
+ invert(p) ==
+ d := degree(p)
+ mapExponents((d-#1)::N, p)
+
+ rootBound(p) ==
+ res : TheField := 1
+ raw :TheField := 1+boundOfCauchy(p)$UTIL
+ while (res < raw) repeat
+ res := 2*(res)
+ res
+
+ sturmNthRoot(lp,l,r,vl,vr,n) ==
+ nv := (vl - vr)::N
+ nv < n => "failed"
+ ((nv = 1) and (n = 1)) => [l,r]
+ int := (l+r)/(2::TheField)
+ lt:List(TheField):=[]
+ for t in lp repeat
+ lt := cons(t.int , lt)
+ vi := sturmVariationsOf(reverse! lt)$UTIL
+ o :Z := n - vl + vi
+ if o > 0
+ then
+ sturmNthRoot(lp,int,r,vi,vr,o::N)
+ else
+ sturmNthRoot(lp,l,int,vl,vi,n)
+
+ sturmIsolate(lp,l,r,vl,vr) ==
+ r <= l => error "ROIRC: sturmIsolate: bad bounds"
+ n := (vl - vr)::N
+ zero?(n) => []
+ one?(n) => [[l,r]]
+ int := (l+r)/(2::TheField)
+ vi := sturmVariationsOf( [t.int for t in lp ] )$UTIL
+ append(sturmIsolate(lp,l,int,vl,vi),sturmIsolate(lp,int,r,vi,vr))
+
+ isolate(lp) ==
+ b := rootBound(first(lp))
+ l1,l2 : List(TheField)
+ (l1,l2) := ([] , [])
+ for t in reverse(lp) repeat
+ if odd?(degree(t))
+ then
+ (l1,l2):= (cons(-leadingCoefficient(t),l1),
+ cons(leadingCoefficient(t),l2))
+ else
+ (l1,l2):= (cons(leadingCoefficient(t),l1),
+ cons(leadingCoefficient(t),l2))
+ sturmIsolate(lp,
+ -b,
+ b,
+ sturmVariationsOf(l1)$UTIL,
+ sturmVariationsOf(l2)$UTIL)
+
+ rootOf(pol,n) ==
+ ls := sturmSequence(pol)$UTIL
+ pol := unitCanonical(first(ls)) -- this one is SqFR
+ degree(pol) = 0 => "failed"
+ numberOfMonomials(pol) = 1 => ([0,1,monomial(1,1)]$Rep)::$
+ b := rootBound(pol)
+ l1,l2 : List(TheField)
+ (l1,l2) := ([] , [])
+ for t in reverse(ls) repeat
+ if odd?(degree(t))
+ then
+ (l1,l2):= (cons(leadingCoefficient(t),l1),
+ cons(-leadingCoefficient(t),l2))
+ else
+ (l1,l2):= (cons(leadingCoefficient(t),l1),
+ cons(leadingCoefficient(t),l2))
+ res := sturmNthRoot(ls,
+ -b,
+ b,
+ sturmVariationsOf(l2)$UTIL,
+ sturmVariationsOf(l1)$UTIL,
+ n)
+ res case "failed" => "failed"
+ makeChar(res.low,res.high,pol)
+
+ allRootsOf(pol) ==
+ ls := sturmSequence(unitCanonical pol)$UTIL
+ pol := unitCanonical(first(ls)) -- this one is SqFR
+ degree(pol) = 0 => []
+ numberOfMonomials(pol) = 1 => [[0,1,monomial(1,1)]$Rep]
+ [ makeChar(term.low,term.high,pol) for term in isolate(ls) ]
+
+
+ hasVarsl(l:List(TheField)) ==
+ null(l) => false
+ f := sign(first(l))
+ for term in rest(l) repeat
+ if f*term < 0 then return(true)
+ false
+
+ hasVars(p:P) ==
+ zero?(p) => error "ROIRC: hasVars: null polynonial"
+ zero?(coefficient(p,0)) => true
+ hasVarsl(coefficients(p))
+
+
+ mightHaveRoots(p,rootChar) ==
+ a := rootChar.low
+ q := translate(p,a)
+ not(hasVars(q)) => false
+-- varStar(q) = 0 => false
+ a := (rootChar.high) - a
+ q := dilate(q,a)
+ sign(coefficient(q,0))*sign(evalOne(q)) <= 0 => true
+ q := minus(addOne(q))
+ not(hasVars(q)) => false
+-- varStar(q) = 0 => false
+ q := invert(q)
+ hasVars(addOne(q))
+-- ^(varStar(addOne(q)) = 0)
+
+ coerce(rootChar:$):O ==
+ commaSeparate([ hconcat("[" :: O , (rootChar.low)::O),
+ hconcat((rootChar.high)::O,"[" ::O ) ])
+
+ c1 = c2 ==
+ mM := max(c1.low,c2.low)
+ Mm := min(c1.high,c2.high)
+ mM >= Mm => false
+ rr : ThePolDom := gcd(c1.defPol,c2.defPol)
+ degree(rr) = 0 => false
+ sign(rr.mM) * sign(rr.Mm) <= 0
+
+ makeChar(left,right,pol) ==
+<<performance problem>>
+ res :$ := [left,right,leadingMonomial(pol)+reductum(pol)]$Rep -- safe copy
+ while zero?(pol.(res.high)) repeat refine!(res)
+ while (res.high * res.low < 0 ) repeat refine!(res)
+ zero?(pol.(res.low)) => [res.low,res.high,monomial(1,1)-(res.low)::P]
+ res
+
+ definingPolynomial(rootChar) == rootChar.defPol
+
+ linearRecip(toTest,rootChar) ==
+ c := - inv(leadingCoefficient(toTest)) * coefficient(toTest,0)
+ r := recip(rootChar.defPol.c)
+ if (r case "failed")
+ then
+ if (c - rootChar.low) * (c - rootChar.high) <= 0
+ then
+ "failed"
+ else
+ newPol := (rootChar.defPol exquo toTest)::P
+ ((1$ThePolDom - inv(newPol.c)*newPol) exquo toTest)::P
+ else
+ ((1$ThePolDom - (r::TheField)*rootChar.defPol) exquo toTest)::P
+
+ recip(toTest,rootChar) ==
+ degree(toTest) = 0 or degree(rootChar.defPol) <= degree(toTest) =>
+ error "IRC: recip: Not reduced"
+ degree(rootChar.defPol) = 1 =>
+ error "IRC: recip: Linear Defining Polynomial"
+ degree(toTest) = 1 =>
+ linearRecip(toTest, rootChar)
+ d := extendedEuclidean((rootChar.defPol),toTest)
+ (degree(d.generator) = 0 ) =>
+ d.coef2
+ d.generator := unitCanonical(d.generator)
+ (d.generator.(rootChar.low) *
+ d.generator.(rootChar.high)<= 0) => "failed"
+ newPol := (rootChar.defPol exquo (d.generator))::P
+ degree(newPol) = 1 =>
+ c := - inv(leadingCoefficient(newPol)) * coefficient(newPol,0)
+ inv(toTest.c)::P
+ degree(toTest) = 1 =>
+ c := - coefficient(toTest,0)/ leadingCoefficient(toTest)
+ ((1$ThePolDom - inv(newPol.(c))*newPol) exquo toTest)::P
+ d := extendedEuclidean(newPol,toTest)
+ d.coef2
+
+ linearSign(toTest,rootChar) ==
+ c := - inv(leadingCoefficient(toTest)) * coefficient(toTest,0)
+ ev := sign(rootChar.defPol.c)
+ if zero?(ev)
+ then
+ if (c - rootChar.low) * (c - rootChar.high) <= 0
+ then
+ 0
+ else
+ sign(toTest.(rootChar.high))
+ else
+ if (ev*sign(rootChar.defPol.(rootChar.high)) <= 0 )
+ then
+ sign(toTest.(rootChar.high))
+ else
+ sign(toTest.(rootChar.low))
+
+ sign(toTest,rootChar) ==
+ degree(toTest) = 0 or degree(rootChar.defPol) <= degree(toTest) =>
+ error "IRC: sign: Not reduced"
+ degree(rootChar.defPol) = 1 =>
+ error "IRC: sign: Linear Defining Polynomial"
+ degree(toTest) = 1 =>
+ linearSign(toTest, rootChar)
+ s := sign(leadingCoefficient(toTest))
+ toTest := monomial(1,degree(toTest))+
+ inv(leadingCoefficient(toTest))*reductum(toTest)
+ delta := gcd(toTest,rootChar.defPol)
+ newChar := [rootChar.low,rootChar.high,rootChar.defPol]$Rep
+ if degree(delta) > 0
+ then
+ if sign(delta.(rootChar.low) * delta.(rootChar.high)) <= 0
+ then
+ return(0)
+ else
+ newChar.defPol := (newChar.defPol exquo delta) :: P
+ toTest := toTest rem (newChar.defPol)
+ degree(toTest) = 0 => s * sign(leadingCoefficient(toTest))
+ degree(toTest) = 1 => s * linearSign(toTest, newChar)
+ while mightHaveRoots(toTest,newChar) repeat
+ newChar := refine(newChar)
+ s*sign(toTest.(newChar.low))
+
+ linearZero?(c,rootChar) ==
+ zero?((rootChar.defPol).c) and
+ (c - rootChar.low) * (c - rootChar.high) <= 0
+
+ zero?(toTest,rootChar) ==
+ degree(toTest) = 0 or degree(rootChar.defPol) <= degree(toTest) =>
+ error "IRC: zero?: Not reduced"
+ degree(rootChar.defPol) = 1 =>
+ error "IRC: zero?: Linear Defining Polynomial"
+ degree(toTest) = 1 =>
+ linearZero?(- inv(leadingCoefficient(toTest)) * coefficient(toTest,0),
+ rootChar)
+ toTest := monomial(1,degree(toTest))+
+ inv(leadingCoefficient(toTest))*reductum(toTest)
+ delta := gcd(toTest,rootChar.defPol)
+ degree(delta) = 0 => false
+ sign(delta.(rootChar.low) * delta.(rootChar.high)) <= 0
+
+
+ refine!(rootChar) ==
+ -- this is not a safe function, it can work with badly created object
+ -- we do not assume (rootChar.defPol).(rootChar.high) <> 0
+ int := middle(rootChar)
+ s1 := sign((rootChar.defPol).(rootChar.low))
+ zero?(s1) =>
+ rootChar.high := int
+ rootChar.defPol := monomial(1,1) - (rootChar.low)::P
+ rootChar
+ s2 := sign((rootChar.defPol).int)
+ zero?(s2) =>
+ rootChar.low := int
+ rootChar.defPol := monomial(1,1) - int::P
+ rootChar
+ if (s1*s2 < 0)
+ then
+ rootChar.high := int
+ else
+ rootChar.low := int
+ rootChar
+
+ refine(rootChar) ==
+ -- we assume (rootChar.defPol).(rootChar.high) <> 0
+ int := middle(rootChar)
+ s:= (rootChar.defPol).int * (rootChar.defPol).(rootChar.high)
+ zero?(s) => [int,rootChar.high,monomial(1,1)-int::P]
+ if s < 0
+ then
+ [int,rootChar.high,rootChar.defPol]
+ else
+ [rootChar.low,int,rootChar.defPol]
+
+ left(rootChar) == rootChar.low
+
+ right(rootChar) == rootChar.high
+
+ middle(rootChar) == (rootChar.low + rootChar.high)/(2::TheField)
+
+-- varStar(p) == -- if 0 no roots in [0,:infty[
+-- res : N := 0
+-- lsg := sign(coefficient(p,0))
+-- l := [ sign(i) for i in reverse!(coefficients(p))]
+-- for sg in l repeat
+-- if (sg ^= lsg) then res := res + 1
+-- lsg := sg
+-- res
+@
+\section{domain RECLOS RealClosure}
+The domain constructore {\bf RealClosure} by Renaud Rioboo (University
+of Paris 6, France) provides the real closure of an ordered field.
+The implementation is based on interval arithmetic. Moreover, the
+design of this constructor and its related packages allows an easy
+use of other codings for real algebraic numbers.
+ordered field
+<<domain RECLOS RealClosure>>=
+)abbrev domain RECLOS RealClosure
+++ Author: Renaud Rioboo
+++ Date Created: summer 1988
+++ Date Last Updated: January 2004
+++ Basic Functions: provides computations in an ordered real closure
+++ Related Constructors: RightOpenIntervalRootCharacterization
+++ Also See:
+++ AMS Classifications:
+++ Keywords: Real Algebraic Numbers
+++ References:
+++ Description:
+++ This domain implements the real closure of an ordered field.
+++ Note:
+++ The code here is generic i.e. it does not depend of the way the operations
+++ are done. The two macros PME and SEG should be passed as functorial
+++ arguments to the domain. It does not help much to write a category
+++ since non trivial methods cannot be placed there either.
+++
+RealClosure(TheField): PUB == PRIV where
+
+ TheField : Join(OrderedRing, Field, RealConstant)
+
+-- ThePols : UnivariatePolynomialCategory($)
+-- PME ==> ThePols
+-- TheCharDom : RealRootCharacterizationCategory($, ThePols )
+-- SEG ==> TheCharDom
+-- this does not work yet
+
+ E ==> OutputForm
+ Z ==> Integer
+ SE ==> Symbol
+ B ==> Boolean
+ SUP ==> SparseUnivariatePolynomial($)
+ N ==> PositiveInteger
+ RN ==> Fraction Z
+ LF ==> ListFunctions2($,N)
+
+-- *****************************************************************
+-- *****************************************************************
+-- PUT YOUR OWN PREFERENCE HERE
+-- *****************************************************************
+-- *****************************************************************
+ PME ==> SparseUnivariatePolynomial($)
+ SEG ==> RightOpenIntervalRootCharacterization($,PME)
+-- *****************************************************************
+-- *****************************************************************
+
+
+ PUB == Join(RealClosedField,
+ FullyRetractableTo TheField,
+ Algebra TheField) with
+
+ algebraicOf : (SEG,E) -> $
+ ++ \axiom{algebraicOf(char)} is the external number
+
+ mainCharacterization : $ -> Union(SEG,"failed")
+ ++ \axiom{mainCharacterization(x)} is the main algebraic
+ ++ quantity of \axiom{x} (\axiom{SEG})
+
+ relativeApprox : ($,$) -> RN
+ ++ \axiom{relativeApprox(n,p)} gives a relative
+ ++ approximation of \axiom{n}
+ ++ that has precision \axiom{p}
+
+ PRIV == add
+
+-- local functions
+
+ lessAlgebraic : $ -> $
+ newElementIfneeded : (SEG,E) -> $
+
+-- Representation
+
+ Rec := Record(seg: SEG, val:PME, outForm:E, order:N)
+ Rep := Union(TheField,Rec)
+
+-- global (mutable) variables
+
+ orderOfCreation : N := 1$N
+ -- it is internally used to sort the algebraic levels
+
+ instanceName : Symbol := new()$Symbol
+ -- this used to print the results, thus different instanciations
+ -- use different names
+
+-- now the code
+
+ relativeApprox(nbe,prec) ==
+ nbe case TheField => retract(nbe)
+ appr := relativeApprox(nbe.val, nbe.seg, prec)
+ -- now appr has the good exact precision but is $
+ relativeApprox(appr,prec)
+
+
+ approximate(nbe,prec) ==
+ abs(nbe) < prec => 0
+ nbe case TheField => retract(nbe)
+ appr := approximate(nbe.val, nbe.seg, prec)
+ -- now appr has the good exact precision but is $
+ approximate(appr,prec)
+
+ newElementIfneeded(s,o) ==
+ p := definingPolynomial(s)
+ degree(p) = 1 =>
+ - coefficient(p,0) / leadingCoefficient(p)
+ res := [s, monomial(1,1), o, orderOfCreation ]$Rec
+ orderOfCreation := orderOfCreation + 1
+ res :: $
+
+ algebraicOf(s,o) ==
+ pol := definingPolynomial(s)
+ degree(pol) = 1 =>
+ -coefficient(pol,0) / leadingCoefficient(pol)
+ res := [s, monomial(1,1), o, orderOfCreation ]$Rec
+ orderOfCreation := orderOfCreation + 1
+ res :: $
+
+ rename!(x,o) ==
+ x.outForm := o
+ x
+
+ rename(x,o) ==
+ [x.seg, x.val, o, x.order]$Rec
+
+ rootOf(pol,n) ==
+ degree(pol) = 0 => "failed"
+ degree(pol) = 1 =>
+ if n=1
+ then
+ -coefficient(pol,0) / leadingCoefficient(pol)
+ else
+ "failed"
+ r := rootOf(pol,n)$SEG
+ r case "failed" => "failed"
+ o := hconcat(instanceName :: E , orderOfCreation :: E)$E
+ algebraicOf(r,o)
+
+ allRootsOf(pol:SUP):List($) ==
+ degree(pol)=0 => []
+ degree(pol)=1 => [-coefficient(pol,0) / leadingCoefficient(pol)]
+ liste := allRootsOf(pol)$SEG
+ res : List $ := []
+ for term in liste repeat
+ o := hconcat(instanceName :: E , orderOfCreation :: E)$E
+ res := cons(algebraicOf(term,o), res)
+ reverse! res
+
+ coerce(x:$):$ ==
+ x case TheField => x
+ [x.seg,x.val rem$PME definingPolynomial(x.seg),x.outForm,x.order]$Rec
+
+ positive?(x) ==
+ x case TheField => positive?(x)$TheField
+ positive?(x.val,x.seg)$SEG
+
+ negative?(x) ==
+ x case TheField => negative?(x)$TheField
+ negative?(x.val,x.seg)$SEG
+
+ abs(x) == sign(x)*x
+
+ sign(x) ==
+ x case TheField => sign(x)$TheField
+ sign(x.val,x.seg)$SEG
+
+ x < y == positive?(y-x)
+
+ x = y == zero?(x-y)
+
+ mainCharacterization(x) ==
+ x case TheField => "failed"
+ x.seg
+
+ mainDefiningPolynomial(x) ==
+ x case TheField => "failed"
+ definingPolynomial x.seg
+
+ mainForm(x) ==
+ x case TheField => "failed"
+ x.outForm
+
+ mainValue(x) ==
+ x case TheField => "failed"
+ x.val
+
+ coerce(x:$):E ==
+ x case TheField => x::TheField :: E
+ xx:$ := coerce(x)
+ outputForm(univariate(xx.val),x.outForm)$SUP
+
+
+ inv(x) ==
+ (res:= recip x) case "failed" => error "Division by 0"
+ res :: $
+
+ recip(x) ==
+ x case TheField =>
+ if ((r := recip(x)$TheField) case TheField)
+ then r::$
+ else "failed"
+ if ((r := recip(x.val,x.seg)$SEG) case "failed")
+ then "failed"
+ else lessAlgebraic([x.seg,r::PME,x.outForm,x.order]$Rec)
+
+ (n:Z * x:$):$ ==
+ x case TheField => n *$TheField x
+ zero?(n) => 0
+ one?(n) => x
+ [x.seg,map(n * #1, x.val),x.outForm,x.order]$Rec
+
+ (rn:TheField * x:$):$ ==
+ x case TheField => rn *$TheField x
+ zero?(rn) => 0
+ one?(rn) => x
+ [x.seg,map(rn * #1, x.val),x.outForm,x.order]$Rec
+
+ (x:$ * y:$):$ ==
+ (x case TheField) and (y case TheField) => x *$TheField y
+ (x case TheField) => x::TheField * y
+ -- x is no longer TheField
+ (y case TheField) => y::TheField * x
+ -- now both are algebraic
+ y.order > x.order =>
+ [y.seg,map(x * #1 , y.val),y.outForm,y.order]$Rec
+ x.order > y.order =>
+ [x.seg,map( #1 * y , x.val),x.outForm,x.order]$Rec
+ -- now x.exp = y.exp
+ -- we will multiply the polynomials and then reduce
+ -- however wee need to call lessAlgebraic
+ lessAlgebraic([x.seg,
+ (x.val * y.val) rem definingPolynomial(x.seg),
+ x.outForm,
+ x.order]$Rec)
+
+ nonNull(rep:Rec):$ ==
+ degree(rep.val)=0 => leadingCoefficient(rep.val)
+ numberOfMonomials(rep.val) = 1 => rep
+ zero?(rep.val,rep.seg)$SEG => 0
+ rep
+
+-- zero?(x) ==
+-- x case TheField => zero?(x)$TheField
+-- zero?(x.val,x.seg)$SEG
+
+ zero?(x) ==
+ x case TheField => zero?(x)$TheField
+ false
+
+ x + y ==
+ (x case TheField) and (y case TheField) => x +$TheField y
+ (x case TheField) =>
+ if zero?(x)
+ then
+ y
+ else
+ nonNull([y.seg,x::PME+(y.val),y.outForm,y.order]$Rec)
+ -- x is no longer TheField
+ (y case TheField) =>
+ if zero?(y)
+ then
+ x
+ else
+ nonNull([x.seg,(x.val)+y::PME,x.outForm,x.order]$Rec)
+ -- now both are algebraic
+ y.order > x.order =>
+ nonNull([y.seg,x::PME+y.val,y.outForm,y.order]$Rec)
+ x.order > y.order =>
+ nonNull([x.seg,(x.val)+y::PME,x.outForm,x.order]$Rec)
+ -- now x.exp = y.exp
+ -- we simply add polynomials (since degree cannot increase)
+ -- however wee need to call lessAlgebraic
+ nonNull([x.seg,x.val + y.val,x.outForm,x.order])
+
+
+ -x ==
+ x case TheField => -$TheField (x::TheField)
+ [x.seg,-$PME x.val,x.outForm,x.order]$Rec
+
+
+ retractIfCan(x:$):Union(TheField,"failed") ==
+ x case TheField => x
+ o := x.order
+ res := lessAlgebraic x
+ res case TheField => res
+ o = res.order => "failed"
+ retractIfCan res
+
+ retract(x:$):TheField ==
+ x case TheField => x
+ o := x.order
+ res := lessAlgebraic x
+ res case TheField => res
+ o = res.order => error "Can't retract"
+ retract res
+
+
+ lessAlgebraic(x) ==
+ x case TheField => x
+ degree(x.val) = 0 => leadingCoefficient(x.val)
+ def := definingPolynomial(x.seg)
+ degree(def) = 1 =>
+ x.val.(- coefficient(def,0) / leadingCoefficient(def))
+ x
+
+ 0 == (0$TheField) :: $
+
+ 1 == (1$TheField) :: $
+
+ coerce(rn:TheField):$ == rn :: $
+
+@
+\section{License}
+<<license>>=
+-----------------------------------------------------------------------------
+-- This software was written by Renaud Rioboo (Computer Algebra group of
+-- Laboratoire d'Informatique de Paris 6) and is the property of university
+-- Paris 6.
+-----------------------------------------------------------------------------
+@
+<<*>>=
+<<license>>
+<<package POLUTIL RealPolynomialUtilitiesPackage>>
+<<category RRCC RealRootCharacterizationCategory>>
+<<category RCFIELD RealClosedField>>
+<<domain ROIRC RightOpenIntervalRootCharacterization>>
+<<domain RECLOS RealClosure>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} R. Rioboo,
+{\sl Real Algebraic Closure of an ordered Field : Implementation in Axiom.},
+In proceedings of the ISSAC'92 Conference, Berkeley 1992 pp. 206-215.
+\bibitem{2} Z. Ligatsikas, R. Rioboo, M. F. Roy
+{\sl Generic computation of the real closure of an ordered field.},
+In Mathematics and Computers in Simulation Volume 42, Issue 4-6,
+November 1996.
+\end{thebibliography}
+\end{document}