diff options
Diffstat (limited to 'src/algebra/realzero.spad.pamphlet')
| -rw-r--r-- | src/algebra/realzero.spad.pamphlet | 347 |
1 files changed, 347 insertions, 0 deletions
diff --git a/src/algebra/realzero.spad.pamphlet b/src/algebra/realzero.spad.pamphlet new file mode 100644 index 00000000..27374130 --- /dev/null +++ b/src/algebra/realzero.spad.pamphlet @@ -0,0 +1,347 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra realzero.spad} +\author{Andy Neff} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{package REAL0 RealZeroPackage} +<<package REAL0 RealZeroPackage>>= +)abbrev package REAL0 RealZeroPackage +++ Author: Andy Neff +++ Date Created: +++ Date Last Updated: +++ Basic Functions: +++ Related Constructors: UnivariatePolynomial, RealZeroPackageQ +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This package provides functions for finding the real zeros +++ of univariate polynomials over the integers to arbitrary user-specified +++ precision. The results are returned as a list of +++ isolating intervals which are expressed as records with "left" and "right" rational number +++ components. + +RealZeroPackage(Pol): T == C where + Pol: UnivariatePolynomialCategory Integer + RN ==> Fraction Integer + Interval ==> Record(left : RN, right : RN) + isoList ==> List(Interval) + T == with + -- next two functions find isolating intervals + realZeros: (Pol) -> isoList + ++ realZeros(pol) returns a list of isolating intervals for + ++ all the real zeros of the univariate polynomial pol. + realZeros: (Pol, Interval) -> isoList + ++ realZeros(pol, range) returns a list of isolating intervals + ++ for all the real zeros of the univariate polynomial pol which + ++ lie in the interval expressed by the record range. + -- next two functions return intervals smaller then tolerence + realZeros: (Pol, RN) -> isoList + ++ realZeros(pol, eps) returns a list of intervals of length less + ++ than the rational number eps for all the real roots of the + ++ polynomial pol. + realZeros: (Pol, Interval, RN) -> isoList + ++ realZeros(pol, int, eps) returns a list of intervals of length + ++ less than the rational number eps for all the real roots of the + ++ polynomial pol which lie in the interval expressed by the + ++ record int. + refine: (Pol, Interval, RN) -> Interval + ++ refine(pol, int, eps) refines the interval int containing + ++ exactly one root of the univariate polynomial pol to size less + ++ than the rational number eps. + refine: (Pol, Interval, Interval) -> Union(Interval,"failed") + ++ refine(pol, int, range) takes a univariate polynomial pol and + ++ and isolating interval int containing exactly one real + ++ root of pol; the operation returns an isolating interval which + ++ is contained within range, or "failed" if no such isolating interval exists. + midpoint: Interval -> RN + ++ midpoint(int) returns the midpoint of the interval int. + midpoints: isoList -> List RN + ++ midpoints(isolist) returns the list of midpoints for the list + ++ of intervals isolist. + C == add + --Local Functions + makeSqfr: Pol -> Pol + ReZeroSqfr: (Pol) -> isoList + PosZero: (Pol) -> isoList + Zero1: (Pol) -> isoList + transMult: (Integer, Pol) -> Pol + transMultInv: (Integer, Pol) -> Pol + transAdd1: (Pol) -> Pol + invert: (Pol) -> Pol + minus: (Pol) -> Pol + negate: Interval -> Interval + rootBound: (Pol) -> Integer + var: (Pol) -> Integer + + negate(int : Interval):Interval == [-int.right,-int.left] + + midpoint(i : Interval):RN == (1/2)*(i.left + i.right) + + midpoints(li : isoList) : List RN == + [midpoint x for x in li] + + makeSqfr(F : Pol):Pol == + sqfr := squareFree F + F := */[s.factor for s in factors(sqfr)] + + realZeros(F : Pol) == + ReZeroSqfr makeSqfr F + + realZeros(F : Pol, rn : RN) == + F := makeSqfr F + [refine(F,int,rn) for int in ReZeroSqfr(F)] + + realZeros(F : Pol, bounds : Interval) == + F := makeSqfr F + [rint::Interval for int in ReZeroSqfr(F) | + (rint:=refine(F,int,bounds)) case Interval] + + realZeros(F : Pol, bounds : Interval, rn : RN) == + F := makeSqfr F + [refine(F,int,rn) for int in realZeros(F,bounds)] + + ReZeroSqfr(F : Pol) == + F = 0 => error "ReZeroSqfr: zero polynomial" + L : isoList := [] + degree(F) = 0 => L + if (r := minimumDegree(F)) > 0 then + L := [[0,0]$Interval] + tempF := F exquo monomial(1, r) + if not (tempF case "failed") then + F := tempF + J:isoList := [negate int for int in reverse(PosZero(minus(F)))] + K : isoList := PosZero(F) + append(append(J, L), K) + + PosZero(F : Pol) == --F is square free, primitive + --and F(0) ^= 0; returns isoList for positive + --roots of F + + b : Integer := rootBound(F) + F := transMult(b,F) + L : isoList := Zero1(F) + int : Interval + L := [[b*int.left, b*int.right]$Interval for int in L] + + Zero1(F : Pol) == --returns isoList for roots of F in (0,1) + J : isoList + K : isoList + L : isoList + L := [] + (v := var(transAdd1(invert(F)))) = 0 => [] + v = 1 => L := [[0,1]$Interval] + G : Pol := transMultInv(2, F) + H : Pol := transAdd1(G) + if minimumDegree H > 0 then + -- H has a root at 0 => F has one at 1/2, and G at 1 + L := [[1/2,1/2]$Interval] + Q : Pol := monomial(1, 1) + tempH : Union(Pol, "failed") := H exquo Q + if not (tempH case "failed") then H := tempH + Q := Q + monomial(-1, 0) + tempG : Union(Pol, "failed") := G exquo Q + if not (tempG case "failed") then G := tempG + int : Interval + J := [[(int.left+1)* (1/2),(int.right+1) * (1/2)]$Interval + for int in Zero1(H)] + K := [[int.left * (1/2), int.right * (1/2)]$Interval + for int in Zero1(G)] + append(append(J, L), K) + + rootBound(F : Pol) == --returns power of 2 that is a bound + --for the positive roots of F + if leadingCoefficient(F) < 0 then F := -F + lcoef := leadingCoefficient(F) + F := reductum(F) + i : Integer := 0 + while not (F = 0) repeat + if (an := leadingCoefficient(F)) < 0 then i := i - an + F := reductum(F) + b : Integer := 1 + while (b * lcoef) <= i repeat + b := 2 * b + b + + transMult(c : Integer, F : Pol) == + --computes Pol G such that G(x) = F(c*x) + G : Pol := 0 + while not (F = 0) repeat + n := degree(F) + G := G + monomial((c**n) * leadingCoefficient(F), n) + F := reductum(F) + G + + transMultInv(c : Integer, F : Pol) == + --computes Pol G such that G(x) = (c**n) * F(x/c) + d := degree(F) + cc : Integer := 1 + G : Pol := monomial(leadingCoefficient F,d) + while (F:=reductum(F)) ^= 0 repeat + n := degree(F) + cc := cc*(c**(d-n):NonNegativeInteger) + G := G + monomial(cc * leadingCoefficient(F), n) + d := n + G + +-- otransAdd1(F : Pol) == +-- --computes Pol G such that G(x) = F(x+1) +-- G : Pol := F +-- n : Integer := 1 +-- while (F := differentiate(F)) ^= 0 repeat +-- if not ((tempF := F exquo n) case "failed") then F := tempF +-- G := G + F +-- n := n + 1 +-- G + + transAdd1(F : Pol) == + --computes Pol G such that G(x) = F(x+1) + n := degree F + v := vectorise(F, n+1) + for i in 0..(n-1) repeat + for j in (n-i)..n repeat + qsetelt_!(v,j, qelt(v,j) + qelt(v,(j+1))) + ans : Pol := 0 + for i in 0..n repeat + ans := ans + monomial(qelt(v,(i+1)),i) + ans + + + minus(F : Pol) == + --computes Pol G such that G(x) = F(-x) + G : Pol := 0 + while not (F = 0) repeat + n := degree(F) + coef := leadingCoefficient(F) + odd? n => + G := G + monomial(-coef, n) + F := reductum(F) + G := G + monomial(coef, n) + F := reductum(F) + G + + invert(F : Pol) == + --computes Pol G such that G(x) = (x**n) * F(1/x) + G : Pol := 0 + n := degree(F) + while not (F = 0) repeat + G := G + monomial(leadingCoefficient(F), + (n-degree(F))::NonNegativeInteger) + F := reductum(F) + G + + var(F : Pol) == --number of sign variations in coefs of F + i : Integer := 0 + LastCoef : Boolean + next : Boolean + LastCoef := leadingCoefficient(F) < 0 + while not ((F := reductum(F)) = 0) repeat + next := leadingCoefficient(F) < 0 + if ((not LastCoef) and next) or + ((not next) and LastCoef) then i := i+1 + LastCoef := next + i + + refine(F : Pol, int : Interval, bounds : Interval) == + lseg := min(int.right,bounds.right) - max(int.left,bounds.left) + lseg < 0 => "failed" + lseg = 0 => + pt := + int.left = bounds.right => int.left + int.right + elt(transMultInv(denom(pt),F),numer pt) = 0 => [pt,pt] + "failed" + lseg = int.right - int.left => int + refine(F, refine(F, int, lseg), bounds) + + refine(F : Pol, int : Interval, eps : RN) == + a := int.left + b := int.right + a=b => [a,b]$Interval + an : Integer := numer(a) + ad : Integer := denom(a) + bn : Integer := numer(b) + bd : Integer := denom(b) + xfl : Boolean := false + if (u:=elt(transMultInv(ad, F), an)) = 0 then + F := (F exquo (monomial(ad,1)-monomial(an,0)))::Pol + u:=elt(transMultInv(ad, F), an) + if (v:=elt(transMultInv(bd, F), bn)) = 0 then + F := (F exquo (monomial(bd,1)-monomial(bn,0)))::Pol + v:=elt(transMultInv(bd, F), bn) + u:=elt(transMultInv(ad, F), an) + if u > 0 then (F:=-F;v:=-v) + if v < 0 then + error [int, "is not a valid isolation interval for", F] + if eps <= 0 then error "precision must be positive" + while (b - a) >= eps repeat + mid : RN := (b + a) * (1/2) + midn : Integer := numer(mid) + midd : Integer := denom(mid) + (v := elt(transMultInv(midd, F), midn)) < 0 => + a := mid + an := midn + ad := midd + v > 0 => + b := mid + bn := midn + bd := midd + v = 0 => + a := mid + b := mid + an := midn + ad := midd + xfl := true + [a, b]$Interval + +@ +\section{License} +<<license>>= +--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. +--All rights reserved. +-- +--Redistribution and use in source and binary forms, with or without +--modification, are permitted provided that the following conditions are +--met: +-- +-- - Redistributions of source code must retain the above copyright +-- notice, this list of conditions and the following disclaimer. +-- +-- - Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in +-- the documentation and/or other materials provided with the +-- distribution. +-- +-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the +-- names of its contributors may be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +@ +<<*>>= +<<license>> + +<<package REAL0 RealZeroPackage>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |