diff options
Diffstat (limited to 'src/algebra/reclos.spad.pamphlet')
| -rw-r--r-- | src/algebra/reclos.spad.pamphlet | 1242 |
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} |