diff options
Diffstat (limited to 'src/algebra/interval.as.pamphlet')
| -rw-r--r-- | src/algebra/interval.as.pamphlet | 564 |
1 files changed, 0 insertions, 564 deletions
diff --git a/src/algebra/interval.as.pamphlet b/src/algebra/interval.as.pamphlet deleted file mode 100644 index 123f1d97..00000000 --- a/src/algebra/interval.as.pamphlet +++ /dev/null @@ -1,564 +0,0 @@ -\documentclass{article} -\usepackage{open-axiom} -\begin{document} -\title{\$SPAD/src/algebra interval.as} -\author{Mike Dewar} -\maketitle -\begin{abstract} -\end{abstract} -\eject -\tableofcontents -\eject -\section{IntervalCategory} -<<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} -<<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} -<<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>> - -<<IntervalCategory>> -<<Interval>> -@ -\eject -\begin{thebibliography}{99} -\bibitem{1} nothing -\end{thebibliography} -\end{document} |