axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / INTRVL.spad

)abbrev domain INTRVL Interval
++ Author: Mike Dewar
++ Date Created: November 1996
++ Description:
++ This domain is an implementation of interval arithmetic and transcendental
++ functions over intervals.

Interval(R) : SIG == CODE where
 R : Join(FloatingPointSystem,TranscendentalFunctionCategory)

 SIG ==> IntervalCategory(R)

 CODE ==> add 

  import Integer

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

  roundDown(u:R):R == 
    if zero?(u) then float(-1,-(bits()@Integer))
                else float(mantissa(u) - 1,exponent(u))

  roundUp(u:R):R   == 
    if zero?(u) then float(1, -(bits())@Integer)
                else float(mantissa(u) + 1,exponent(u))

  -- Sometimes the float representation does not use all the bits (for when
  -- example, 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**((b-l) pretend PositiveInteger),exponent(u)-b+l)
    else
      u

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

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

  qinterval(i:R,s:R):% ==
    [roundDown normaliseFloat i,roundUp normaliseFloat 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 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@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@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
   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):% == 
    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()@Integer) => 
      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):% == interval(inf(u)**n,sup(u)**n)