Edit detail for Interval Arithmetic revision 3 of 3

changed:
-Author -- Mike Dewar
-
-Date Created -- November 1996
-
-Description: --
-This category is an implementation of interval arithmetic and
-transcendental functions over intervals.
-
-**Note:** -- It appears to be nearly identical to the SPAD category
-'interval.spad'.
-
-\begin{aldor}[interval2]
-
-#include "axiom.as"
-
-FUNCAT ==> Join(FloatingPointSystem,TranscendentalFunctionCategory);
-
-define IntervalCategory2(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.
-}
-
-Interval2(R:FUNCAT): IntervalCategory2(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);
-  
-}
-\end{aldor}
Interval arithmetic for Axiom was implemented by Mike Dewar in form of Aldor files.
Juergen Weiss converted the files to Spad and they are included in current FriCAS
distribution.



changed:
-x:=interval(1.1,2.2)$Interval2(Float)
x:=interval(1.1,2.2)$Interval(Float)

removed:
-From unknown Fri Sep 2 02:26:28 -0500 2005
-From: unknown
-Date: Fri, 02 Sep 2005 02:26:28 -0500
-Subject: 
-Message-ID: <20050902022628-0500@www.axiom-developer.org>
-
-Reason that it looks like interval.spad is that I rewrote interval.as for the SPAD compiler, so that we are able to run Axiom without Aldor. Juergen Weiss
-
-From BillPage Fri Sep 2 22:14:49 -0500 2005
-From: Bill Page
-Date: Fri, 02 Sep 2005 22:14:49 -0500
-Subject: Thanks
-Message-ID: <20050902221449-0500@www.axiom-developer.org>
-
-Juergen Weiss wrote:
-
->  I rewrote interval.as for the SPAD compiler
-
-Oh yes, I am sorry that I forgot. I think your important contribution
-should be so noted in the interval.spad.pamphlet file! Here's a
-patch::
-
-  $ diff -Nar interval.spad.pamphlet.orig interval.spad.pamphlet
-  5c5
-  < \author{Mike Dewar}
-  ---
-  > \author{Juergen Weiss}
-  7a8,9
-  > This is a re-write of `interval.as' in SPAD to enable all of
-  > Axiom's library to be compiled without the Aldor compiler.
-  15c17
-  < +++ Author: Mike Dewar
-  ---
-  > +++ Original Author: Mike Dewar
-
---------

Interval arithmetic for Axiom was implemented by Mike Dewar in form of Aldor files. Juergen Weiss converted the files to Spad and they are included in current FriCAS distribution.

fricas

(1) -> x:=interval(1.1,2.2)$Interval(Float)

(1)

Type: Interval(Float)

fricas

sin(x)

(2)

Type: Interval(Float)