interval.spad

)abbrev category INTCAT IntervalCategory
+++ Author: Mike Dewar
+++ Date Created: November 1996
+++ Basic Functions:
+++ Related Constructors:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ References:
+++ Description:
+++ This category implements of interval arithmetic and transcendental
+++ functions over intervals.
IntervalCategory(R : Join(FloatingPointSystem, TranscendentalFunctionCategory)):
 Category == Join(GcdDomain, OrderedSet, TranscendentalFunctionCategory,
                  RadicalCategory, RetractableTo(Integer), Approximate) with
  interval : (R, R) -> %
    ++ interval(inf, sup) creates a new interval, either \spad{[inf, sup]} if
    ++ \spad{inf <= sup} or \spad{[sup, inf]} otherwise.
  qinterval : (R, R) -> %
    ++ qinterval(inf, sup) creates a new interval \spad{[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 \spad{u}.
  sup : % -> R
    ++ sup(u) returns the supremum of \spad{u}.
  width : % -> R
    ++ width(u) returns \spad{sup(u) - inf(u)}.
  positive? : % -> Boolean
    ++ positive?(u) returns \spad{true} if every element of u is positive,
    ++ \spad{false} otherwise.
  negative? : % -> Boolean
    ++ negative?(u) returns \spad{true} if every element of u is negative,
    ++ \spad{false} otherwise.
  contains? : (%, R) -> Boolean
    ++ contains?(i, f) returns true if \spad{f} is contained within the interval
    ++ \spad{i}, false otherwise.

)abbrev domain INTRVL Interval
+++ Author: Mike Dewar
+++ Date Created: November 1996
+++ 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 : Join(FloatingPointSystem, TranscendentalFunctionCategory)) : IntervalCategory(R) == add

  import from Integer
--  import from R

  Rep := Record(Inf : R, Sup : R)

  roundDown(u : R) : R ==
      zero?(u) => u
      prec := bits()
      m := mantissa u
      l := length(m)
      eu := exponent(u)
      l = prec => float(m - 1, eu)
      dig : Integer := base()
      l < prec =>
          float(m*dig^(qcoerce(prec - l)@PositiveInteger) - 1, eu - prec + l)
      float(m - dig^(qcoerce(l - prec)@ PositiveInteger), eu)

  roundUp(u : R) : R   ==
      zero?(u) => u
      prec := bits()
      m := mantissa u
      l := length(m)
      eu := exponent(u)
      l = prec => float(m + 1, eu)
      dig : Integer := base()
      l < prec =>
          float(m*dig^(qcoerce(prec - l)@PositiveInteger) + 1, eu - prec + l)
      float(m + dig^(qcoerce(l - prec)@PositiveInteger), eu)

  -- 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.
  normaliseFloat(u : R) : R ==
    zero? u => u
    m : Integer := mantissa u
    b : Integer := bits() :: Integer
    l : Integer := length(m)
    if l < b then
      BASE : Integer := base()$R :: Integer
      float(m*BASE^(qcoerce(b - l)@PositiveInteger), exponent(u) - b + l)
    else
      u

  interval(i : R, s : R) : % ==
    i > s =>  [roundDown s, roundUp i]
    [roundDown i, roundUp s]

  interval(f : R) : % ==
    zero?(f) => 0
    one?(f)  => 1
    -- This next part is necessary to allow e.g. mapping between Expressions:
    -- FriCAS assumes that Integers stay as Integers! FIXME
--    import from Union(value1:Integer,failed:"failed")
    retractIfCan(f)@Union(Integer,"failed") case "failed" =>
        [roundDown f, roundUp f]
    [f, f]

  qinterval(i : R, s : R) : % ==
    [roundDown i, roundUp s]

  exactInterval(i : R, s : R) : % == [i, s]
  exactSupInterval(i : R, s : R) : % == [roundDown i, s]
  exactInfInterval(i : R, s : R) : % == [i, roundUp s]

  inf(u : %) : R == u.Inf
  sup(u : %) : R == u.Sup
  width(u : %) : R == u.Sup - 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 blatant 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)
    qinterval(inf(a) - sup(b), sup(a) - inf(b))


  _* (a : %, b : %) : % ==
    -- A couple of blatant 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@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
        [one, one]

  0 == [0, 0]

  recip(u:%):Union(%,"failed") ==
   contains?(u,0) => "failed"
   vals : List R := sort [1/inf(u), 1/sup(u)]$List(R)
   qinterval(first vals, last vals)


  unit?(u : %) : Boolean == contains?(u, 0)

  _exquo(u:%,v:%):Union(%,"failed") ==
   contains?(v,0) => "failed"
   one?(v) => u
   zero?(u) => 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)]$List(R)
   qinterval(first vals, last vals)


  gcd(u : %, v : %) : % == 1

  coerce(u : Integer) : % ==
    ur := normaliseFloat(u::R)
    exactInterval(ur, ur)


  interval(u : Fraction Integer) : % ==
--    import from   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
    bin : Union(Integer,"failed") := retractIfCan(log2(denom(u)::Float))
    bin case Integer and length(numer u)$Integer < bits() =>
      flt := normaliseFloat flt
      exactInterval(flt, flt)

    qinterval(flt, flt)


  retractIfCan(u:%):Union(Integer,"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'.
  hasTwoPiMultiple(offset : R, ipi : R, i : %) : Boolean ==
    next : Integer := retract ceiling( (inf(i) - offset)/(2*ipi) )
    contains?(i, offset+2*next*ipi)


  -- This function checks whether an interval contains a value of the form
  -- `offset + n pi'.
  hasPiMultiple(offset : R, ipi : R, i : %) : Boolean ==
    next : Integer := retract ceiling( (inf(i) - offset)/ipi )
    contains?(i, offset+next*ipi)


  sin(u : %) : % ==
    ipi : R := pi()$R
    hasOne? : Boolean := hasTwoPiMultiple(ipi/(2::R), ipi, u)
    hasMinusOne? : Boolean := hasTwoPiMultiple(3*ipi/(2::R), ipi, 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 : %) : % ==
    ipi : R := pi()
    hasOne? : Boolean := hasTwoPiMultiple(0, ipi, u)
    hasMinusOne? : Boolean := hasTwoPiMultiple(ipi, ipi, 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 : %) : % ==
    ipi : R := pi()
    if width(u) > ipi 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 : %) : % ==
    ipi : R := pi()
    if width(u) > ipi then
      error "Interval contains a singularity"
    else
--      import from Integer
      -- singularities are at multiples of Pi
      if hasPiMultiple(0,ipi,u) then error "Interval contains a singularity"
      vals : List R := sort [csc inf u, csc sup u]
      if hasTwoPiMultiple(ipi/(2::R), ipi, u) then
        exactInfInterval(1, last vals)
      else if hasTwoPiMultiple(3*ipi/(2::R), ipi, u) then
        exactSupInterval(first vals, -1)
      else
        qinterval(first vals, last vals)



  sec(u : %) : % ==
    ipi : R := pi()
    if width(u) > ipi then
      error "Interval contains a singularity"
    else
--      import from Integer
      -- singularities are at Pi/2 + n Pi
      if hasPiMultiple(ipi/(2::R), ipi, u) then
        error "Interval contains a singularity"
      vals : List R := sort [sec inf u, sec sup u]
      if hasTwoPiMultiple(0, ipi, u) then
        exactInfInterval(1, last vals)
      else if hasTwoPiMultiple(ipi, ipi, u) then
        exactSupInterval(first vals, -1)
      else
        qinterval(first vals, last vals)




  cot(u : %) : % ==
    ipi : R := pi()
    if width(u) > ipi 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) : % ==
      lo : R := inf(u)
      hi : R := sup(u)
      lo <= 0 =>
           even?(denom(n)) =>
               error "fractional power only defined for x > 0"
           even?(numer(n)) =>
               hi < 0 =>
                   interval(0, lo^n)
               interval(0, max(lo^n, hi^n))
           interval(lo^n, hi^n)
      interval(lo^n, hi^n)

--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.