aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/interval.as.pamphlet
diff options
context:
space:
mode:
Diffstat (limited to 'src/algebra/interval.as.pamphlet')
-rw-r--r--src/algebra/interval.as.pamphlet564
1 files changed, 564 insertions, 0 deletions
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}
+<<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}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-29 11:36:56 +0000