From ab8cc85adde879fb963c94d15675783f2cf4b183 Mon Sep 17 00:00:00 2001 From: dos-reis Date: Tue, 14 Aug 2007 05:14:52 +0000 Subject: Initial population. --- src/algebra/interval.as.pamphlet | 564 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 564 insertions(+) create mode 100644 src/algebra/interval.as.pamphlet (limited to 'src/algebra/interval.as.pamphlet') diff --git a/src/algebra/interval.as.pamphlet b/src/algebra/interval.as.pamphlet new file mode 100644 index 00000000..765d7327 --- /dev/null +++ b/src/algebra/interval.as.pamphlet @@ -0,0 +1,564 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra interval.as} +\author{Mike Dewar} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{IntervalCategory} +<>= +#include "axiom.as" + ++++ Author: Mike Dewar ++++ Date Created: November 1996 ++++ Date Last Updated: ++++ Basic Functions: ++++ Related Constructors: ++++ Also See: ++++ AMS Classifications: ++++ Keywords: ++++ References: ++++ Description: ++++ This category is an implementation of interval arithmetic and transcendental ++++ functions over intervals. + +FUNCAT ==> Join(FloatingPointSystem,TranscendentalFunctionCategory); + +define IntervalCategory(R:FUNCAT): Category == + Join(GcdDomain, OrderedSet, TranscendentalFunctionCategory, RadicalCategory, + RetractableTo(Integer)) + with { + approximate; + interval : (R,R) -> %; + ++ interval(inf,sup) creates a new interval, either \axiom{[inf,sup]} if + ++ \axiom{inf <= sup} or \axiom{[sup,in]} otherwise. + qinterval : (R,R) -> %; + ++ qinterval(inf,sup) creates a new interval \axiom{[inf,sup]}, without + ++ checking the ordering on the elements. + interval : R -> %; + ++ interval(f) creates a new interval around f. + interval : Fraction Integer -> %; + ++ interval(f) creates a new interval around f. + inf : % -> R; + ++ inf(u) returns the infinum of \axiom{u}. + sup : % -> R; + ++ sup(u) returns the supremum of \axiom{u}. + width : % -> R; + ++ width(u) returns \axiom{sup(u) - inf(u)}. + positive? : % -> Boolean; + ++ positive?(u) returns \axiom{true} if every element of u is positive, + ++ \axiom{false} otherwise. + negative? : % -> Boolean; + ++ negative?(u) returns \axiom{true} if every element of u is negative, + ++ \axiom{false} otherwise. + contains? : (%,R) -> Boolean; + ++ contains?(i,f) returns true if \axiom{f} is contained within the interval + ++ \axiom{i}, false otherwise. +} + +@ +\section{Interval} +<>= ++++ Author: Mike Dewar ++++ Date Created: November 1996 ++++ Date Last Updated: ++++ Basic Functions: ++++ Related Constructors: ++++ Also See: ++++ AMS Classifications: ++++ Keywords: ++++ References: ++++ Description: ++++ This domain is an implementation of interval arithmetic and transcendental ++++ functions over intervals. + +Interval(R:FUNCAT): IntervalCategory(R) == add { + + import from Integer; + import from R; + + Rep ==> Record(Inf:R, Sup:R); + + import from Rep; + + local roundDown(u:R):R == + if zero?(u) then float(-1,-(bits() pretend Integer)); + else float(mantissa(u) - 1,exponent(u)); + + local roundUp(u:R):R == + if zero?(u) then float(1, -(bits()) pretend Integer); + else float(mantissa(u) + 1,exponent(u)); + + -- Sometimes the float representation does not use all the bits (e.g. when + -- representing an integer in software using arbitrary-length Integers as + -- your mantissa it is convenient to keep them exact). This function + -- normalises things so that rounding etc. works as expected. It is only + -- called when creating new intervals. + local normaliseFloat(u:R):R == + if zero? u then u else { + m : Integer := mantissa u; + b : Integer := bits() pretend Integer; + l : Integer := length(m); + if (l < b) then { + BASE : Integer := base()$R pretend Integer; + float(m*BASE**((b-l) pretend PositiveInteger),exponent(u)-b+l); + } + else + u; + } + + interval(i:R,s:R):% == { + i > s => per [roundDown normaliseFloat s,roundUp normaliseFloat i]; + per [roundDown normaliseFloat i,roundUp normaliseFloat s]; + } + + interval(f:R):% == { + zero?(f) => 0; + one?(f) => 1; + -- This next part is necessary to allow e.g. mapping between Expressions: + -- AXIOM assumes that Integers stay as Integers! + import from Union(value1:Integer,failed:'failed'); + fnew : R := normaliseFloat f; + retractIfCan(f)@Union(value1:Integer,failed:'failed') case value1 => + per [fnew,fnew]; + per [roundDown fnew, roundUp fnew]; + } + + qinterval(i:R,s:R):% == + per [roundDown normaliseFloat i,roundUp normaliseFloat s]; + + local exactInterval(i:R,s:R):% == per [i,s]; + local exactSupInterval(i:R,s:R):% == per [roundDown i,s]; + local exactInfInterval(i:R,s:R):% == per [i,roundUp s]; + + inf(u:%):R == (rep u).Inf; + sup(u:%):R == (rep u).Sup; + width(u:%):R == (rep u).Sup - (rep u).Inf; + + contains?(u:%,f:R):Boolean == (f > inf(u)) and (f < sup(u)); + + positive?(u:%):Boolean == inf(u) > 0; + negative?(u:%):Boolean == sup(u) < 0; + + (<)(a:%,b:%):Boolean == + if inf(a) < inf(b) then + true + else if inf(a) > inf(b) then + false + else + sup(a) < sup(b); + + (+)(a:%,b:%):% == { + -- A couple of blatent hacks to preserve the Ring Axioms! + if zero?(a) then return(b) else if zero?(b) then return(a); + if a=b then return qinterval(2*inf(a),2*sup(a)); + qinterval(inf(a) + inf(b), sup(a) + sup(b)); + } + + (-)(a:%,b:%):% == { + if zero?(a) then return(-b) else if zero?(b) then return(a); + if a=b then 0 else qinterval(inf(a) - sup(b), sup(a) - inf(b)); + } + + (*)(a:%,b:%):% == { + -- A couple of blatent hacks to preserve the Ring Axioms! + if one?(a) then return(b) else if one?(b) then return(a); + if zero?(a) then return(0) else if zero?(b) then return(0); + prods : List R := sort [inf(a)*inf(b),sup(a)*sup(b), + inf(a)*sup(b),sup(a)*inf(b)]; + qinterval(first prods, last prods); + } + + (*)(a:Integer,b:%):% == { + if (a > 0) then + qinterval(a*inf(b),a*sup(b)); + else if (a < 0) then + qinterval(a*sup(b),a*inf(b)); + else + 0; + } + + (*)(a:PositiveInteger,b:%):% == qinterval(a*inf(b),a*sup(b)); + + (**)(a:%,n:PositiveInteger):% == { + contains?(a,0) and zero?((n pretend Integer) rem 2) => + interval(0,max(inf(a)**n,sup(a)**n)); + interval(inf(a)**n,sup(a)**n); + } + + (^) (a:%,n:PositiveInteger):% == { + contains?(a,0) and zero?((n pretend Integer) rem 2) => + interval(0,max(inf(a)**n,sup(a)**n)); + interval(inf(a)**n,sup(a)**n); + } + + (-)(a:%):% == exactInterval(-sup(a),-inf(a)); + + (=)(a:%,b:%):Boolean == (inf(a)=inf(b)) and (sup(a)=sup(b)); + (~=)(a:%,b:%):Boolean == (inf(a)~=inf(b)) or (sup(a)~=sup(b)); + + 1:% == {one : R := normaliseFloat 1; per([one,one])}; + 0:% == per([0,0]); + + recip(u:%):Union(value1:%,failed:'failed') == { + contains?(u,0) => [failed]; + vals:List R := sort[1/inf(u),1/sup(u)]; + [qinterval(first vals, last vals)]; + } + + unit?(u:%):Boolean == contains?(u,0); + + exquo(u:%,v:%):Union(value1:%,failed:'failed') == { + contains?(v,0) => [failed]; + one?(v) => [u]; + u=v => [1]; + u=-v => [-1]; + vals:List R := sort[inf(u)/inf(v),inf(u)/sup(v),sup(u)/inf(v),sup(u)/sup(v)]; + [qinterval(first vals, last vals)]; + } + + gcd(u:%,v:%):% == 1; + + coerce(u:Integer):% == { + ur := normaliseFloat(u::R); + exactInterval(ur,ur); + } + + interval(u:Fraction Integer):% == { + import { log2 : % -> %; + coerce : Integer -> %; + retractIfCan : % -> Union(value1:Integer,failed:'failed');} + from Float; + flt := u::R; + + -- Test if the representation in R is exact + --den := denom(u)::Float; + local bin : Union(value1:Integer,failed:'failed'); + bin := retractIfCan(log2(denom(u)::Float)); + bin case value1 and length(numer u)$Integer < (bits() pretend Integer) => { + flt := normaliseFloat flt; + exactInterval(flt,flt); + } + + qinterval(flt,flt); + } + + retractIfCan(u:%):Union(value1:Integer,failed:'failed') == { + not zero? width(u) => [failed]; + retractIfCan inf u; + } + + retract(u:%):Integer == { + not zero? width(u) => + error "attempt to retract a non-Integer interval to an Integer"; + retract inf u; + } + + coerce(u:%):OutputForm == + bracket([coerce inf(u), coerce sup(u)]$List(OutputForm)); + + characteristic():NonNegativeInteger == 0; + + + -- Explicit export from TranscendentalFunctionCategory + pi():% == qinterval(pi(),pi()); + + -- From ElementaryFunctionCategory + log(u:%):% == { + positive?(u) => qinterval(log inf u, log sup u); + error "negative logs in interval"; + } + + exp(u:%):% == qinterval(exp inf u, exp sup u); + + (**)(u:%,v:%):% == { + zero?(v) => if zero?(u) then error "0**0 is undefined" else 1; + one?(u) => 1; + expts : List R := sort [inf(u)**inf(v),sup(u)**sup(v), + inf(u)**sup(v),sup(u)**inf(v)]; + qinterval(first expts, last expts); + } + + -- From TrigonometricFunctionCategory + + -- This function checks whether an interval contains a value of the form + -- `offset + 2 n pi'. + local hasTwoPiMultiple(offset:R,Pi:R,i:%):Boolean == { + import from Integer; + next : Integer := retract ceiling( (inf(i) - offset)/(2*Pi) ); + contains?(i,offset+2*next*Pi); + } + + -- This function checks whether an interval contains a value of the form + -- `offset + n pi'. + local hasPiMultiple(offset:R,Pi:R,i:%):Boolean == { + import from Integer; + next : Integer := retract ceiling( (inf(i) - offset)/Pi ); + contains?(i,offset+next*Pi); + } + + sin(u:%):% == { + import from Integer; + Pi : R := pi(); + hasOne? : Boolean := hasTwoPiMultiple(Pi/(2::R),Pi,u); + hasMinusOne? : Boolean := hasTwoPiMultiple(3*Pi/(2::R),Pi,u); + + if hasOne? and hasMinusOne? then + exactInterval(-1,1); + else { + vals : List R := sort [sin inf u, sin sup u]; + if hasOne? then + exactSupInterval(first vals, 1); + else if hasMinusOne? then + exactInfInterval(-1,last vals); + else + qinterval(first vals, last vals); + } + } + + cos(u:%):% == { + Pi : R := pi(); + hasOne? : Boolean := hasTwoPiMultiple(0,Pi,u); + hasMinusOne? : Boolean := hasTwoPiMultiple(Pi,Pi,u); + + if hasOne? and hasMinusOne? then + exactInterval(-1,1); + else { + vals : List R := sort [cos inf u, cos sup u]; + if hasOne? then + exactSupInterval(first vals, 1); + else if hasMinusOne? then + exactInfInterval(-1,last vals); + else + qinterval(first vals, last vals); + } + } + + tan(u:%):% == { + Pi : R := pi(); + if width(u) > Pi then + error "Interval contains a singularity" + else { + -- Since we know the interval is less than pi wide, monotonicity implies + -- that there is no singularity. If there is a singularity on a endpoint + -- of the interval the user will see the error generated by R. + lo : R := tan inf u; + hi : R := tan sup u; + + lo > hi => error "Interval contains a singularity"; + qinterval(lo,hi); + } + } + + csc(u:%):% == { + Pi : R := pi(); + if width(u) > Pi then + error "Interval contains a singularity" + else { + import from Integer; + -- singularities are at multiples of Pi + if hasPiMultiple(0,Pi,u) then error "Interval contains a singularity"; + vals : List R := sort [csc inf u, csc sup u]; + if hasTwoPiMultiple(Pi/(2::R),Pi,u) then + exactInfInterval(1,last vals); + else if hasTwoPiMultiple(3*Pi/(2::R),Pi,u) then + exactSupInterval(first vals,-1); + else + qinterval(first vals, last vals); + } + } + + sec(u:%):% == { + Pi : R := pi(); + if width(u) > Pi then + error "Interval contains a singularity" + else { + import from Integer; + -- singularities are at Pi/2 + n Pi + if hasPiMultiple(Pi/(2::R),Pi,u) then + error "Interval contains a singularity"; + vals : List R := sort [sec inf u, sec sup u]; + if hasTwoPiMultiple(0,Pi,u) then + exactInfInterval(1,last vals); + else if hasTwoPiMultiple(Pi,Pi,u) then + exactSupInterval(first vals,-1); + else + qinterval(first vals, last vals); + } + } + + + cot(u:%):% == { + Pi : R := pi(); + if width(u) > Pi then + error "Interval contains a singularity" + else { + -- Since we know the interval is less than pi wide, monotonicity implies + -- that there is no singularity. If there is a singularity on a endpoint + -- of the interval the user will see the error generated by R. + hi : R := cot inf u; + lo : R := cot sup u; + + lo > hi => error "Interval contains a singularity"; + qinterval(lo,hi); + } + } + + -- From ArcTrigonometricFunctionCategory + + asin(u:%):% == { + lo : R := inf(u); + hi : R := sup(u); + if (lo < -1) or (hi > 1) then error "asin only defined on the region -1..1"; + qinterval(asin lo,asin hi); + } + + acos(u:%):% == { + lo : R := inf(u); + hi : R := sup(u); + if (lo < -1) or (hi > 1) then error "acos only defined on the region -1..1"; + qinterval(acos hi,acos lo); + } + + atan(u:%):% == qinterval(atan inf u, atan sup u); + + acot(u:%):% == qinterval(acot sup u, acot inf u); + + acsc(u:%):% == { + lo : R := inf(u); + hi : R := sup(u); + if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then + error "acsc not defined on the region -1..1"; + qinterval(acsc hi, acsc lo); + } + + asec(u:%):% == { + lo : R := inf(u); + hi : R := sup(u); + if ((lo < -1) and (hi > -1)) or ((lo < 1) and (hi > 1)) then + error "asec not defined on the region -1..1"; + qinterval(asec lo, asec hi); + } + + -- From HyperbolicFunctionCategory + + tanh(u:%):% == qinterval(tanh inf u, tanh sup u); + + sinh(u:%):% == qinterval(sinh inf u, sinh sup u); + + sech(u:%):% == { + negative? u => qinterval(sech inf u, sech sup u); + positive? u => qinterval(sech sup u, sech inf u); + vals : List R := sort [sech inf u, sech sup u]; + exactSupInterval(first vals,1); + } + + cosh(u:%):% == { + negative? u => qinterval(cosh sup u, cosh inf u); + positive? u => qinterval(cosh inf u, cosh sup u); + vals : List R := sort [cosh inf u, cosh sup u]; + exactInfInterval(1,last vals); + } + + csch(u:%):% == { + contains?(u,0) => error "csch: singularity at zero"; + qinterval(csch sup u, csch inf u); + } + + coth(u:%):% == { + contains?(u,0) => error "coth: singularity at zero"; + qinterval(coth sup u, coth inf u); + } + + -- From ArcHyperbolicFunctionCategory + + acosh(u:%):% == { + inf(u)<1 => error "invalid argument: acosh only defined on the region 1.."; + qinterval(acosh inf u, acosh sup u); + } + + acoth(u:%):% == { + lo : R := inf(u); + hi : R := sup(u); + if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then + error "acoth not defined on the region -1..1"; + qinterval(acoth hi, acoth lo); + } + + acsch(u:%):% == { + contains?(u,0) => error "acsch: singularity at zero"; + qinterval(acsch sup u, acsch inf u); + } + + asech(u:%):% == { + lo : R := inf(u); + hi : R := sup(u); + if (lo <= 0) or (hi > 1) then + error "asech only defined on the region 0 < x <= 1"; + qinterval(asech hi, asech lo); + } + + asinh(u:%):% == qinterval(asinh inf u, asinh sup u); + + atanh(u:%):% == { + lo : R := inf(u); + hi : R := sup(u); + if (lo <= -1) or (hi >= 1) then + error "atanh only defined on the region -1 < x < 1"; + qinterval(atanh lo, atanh hi); + } + + -- From RadicalCategory + (**)(u:%,n:Fraction Integer):% == interval(inf(u)**n,sup(u)**n); + +} + +@ +\section{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. +@ +<<*>>= +<> + +<> +<> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} -- cgit v1.2.3