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

)abbrev package NUMERIC Numeric
++ Author: Manuel Bronstein
++ Date Created: 21 Feb 1990
++ Date Last Updated: 17 August 1995, Mike Dewar
++                    24 January 1997, Miked Dewar (added partial operators)
++ Basic Operations: numeric, complexNumeric, numericIfCan, complexNumericIfCan
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: 
++ References:
++ Description: Numeric provides real and complex numerical evaluation
++ functions for various symbolic types.
 
Numeric(S:ConvertibleTo Float): with
  numeric: S -> Float
    ++ numeric(x) returns a real approximation of x.
  numeric: (S, PositiveInteger) -> Float
    ++ numeric(x, n) returns a real approximation of x up to n decimal
    ++ places.
  complexNumeric: S -> Complex Float
    ++ complexNumeric(x) returns a complex approximation of x.
  complexNumeric: (S, PositiveInteger) -> Complex Float
    ++ complexNumeric(x, n) returns a complex approximation of x up
    ++ to n decimal places.
  if S has CommutativeRing then
    complexNumeric: Complex S -> Complex Float
      ++ complexNumeric(x) returns a complex approximation of x.
    complexNumeric: (Complex S, PositiveInteger) -> Complex Float
      ++ complexNumeric(x, n) returns a complex approximation of x up
      ++ to n decimal places.
    complexNumeric: Polynomial Complex S -> Complex Float
      ++ complexNumeric(x) returns a complex approximation of x.
    complexNumeric: (Polynomial Complex S, PositiveInteger) -> Complex Float
      ++ complexNumeric(x, n) returns a complex approximation of x up
      ++ to n decimal places.
  if S has Ring then
    numeric: Polynomial S -> Float
      ++ numeric(x) returns a real approximation of x.
    numeric: (Polynomial S, PositiveInteger) -> Float
      ++ numeric(x,n) returns a real approximation of x up to n decimal
      ++ places.
    complexNumeric: Polynomial S -> Complex Float
      ++ complexNumeric(x) returns a complex approximation of x.
    complexNumeric: (Polynomial S, PositiveInteger) -> Complex Float
      ++ complexNumeric(x, n) returns a complex approximation of x
      ++ up to n decimal places.
  if S has IntegralDomain then
    numeric: Fraction Polynomial S -> Float
      ++ numeric(x) returns a real approximation of x.
    numeric: (Fraction Polynomial S, PositiveInteger) -> Float
      ++ numeric(x,n) returns a real approximation of x up to n decimal
      ++ places.
    complexNumeric: Fraction Polynomial S -> Complex Float
      ++ complexNumeric(x) returns a complex approximation of x.
    complexNumeric: (Fraction Polynomial S, PositiveInteger) -> Complex Float
      ++ complexNumeric(x, n) returns a complex approximation of x
    complexNumeric: Fraction Polynomial Complex S -> Complex Float
      ++ complexNumeric(x) returns a complex approximation of x.
    complexNumeric: (Fraction Polynomial Complex S, PositiveInteger) -> 
                                                                Complex Float
      ++ complexNumeric(x, n) returns a complex approximation of x
      ++ up to n decimal places.
    if S has OrderedSet then
      numeric: Expression S -> Float
        ++ numeric(x) returns a real approximation of x.
      numeric: (Expression S, PositiveInteger) -> Float
        ++ numeric(x, n) returns a real approximation of x up to n
        ++ decimal places.
      complexNumeric: Expression S -> Complex Float
        ++ complexNumeric(x) returns a complex approximation of x.
      complexNumeric: (Expression S, PositiveInteger) -> Complex Float
        ++ complexNumeric(x, n) returns a complex approximation of x
        ++ up to n decimal places.
      complexNumeric: Expression Complex S -> Complex Float
        ++ complexNumeric(x) returns a complex approximation of x.
      complexNumeric: (Expression Complex S, PositiveInteger) -> Complex Float
        ++ complexNumeric(x, n) returns a complex approximation of x
        ++ up to n decimal places.
  if S has CommutativeRing then
    complexNumericIfCan: Polynomial Complex S -> Union(Complex Float,"failed")
      ++ complexNumericIfCan(x) returns a complex approximation of x,
      ++ or "failed" if \axiom{x} is not constant.
    complexNumericIfCan: (Polynomial Complex S, PositiveInteger) -> Union(Complex Float,"failed")
      ++ complexNumericIfCan(x, n) returns a complex approximation of x up
      ++ to n decimal places, or "failed" if \axiom{x} is not a constant.
  if S has Ring then
    numericIfCan: Polynomial S -> Union(Float,"failed")
      ++ numericIfCan(x) returns a real approximation of x,
      ++ or "failed" if \axiom{x} is not a constant.
    numericIfCan: (Polynomial S, PositiveInteger) -> Union(Float,"failed")
      ++ numericIfCan(x,n) returns a real approximation of x up to n decimal
      ++ places, or "failed" if \axiom{x} is not a constant.
    complexNumericIfCan: Polynomial S -> Union(Complex Float,"failed")
      ++ complexNumericIfCan(x) returns a complex approximation of x,
      ++ or "failed" if \axiom{x} is not a constant.
    complexNumericIfCan: (Polynomial S, PositiveInteger) -> Union(Complex Float,"failed")
      ++ complexNumericIfCan(x, n) returns a complex approximation of x
      ++ up to n decimal places, or "failed" if \axiom{x} is not a constant.
  if S has IntegralDomain then
    numericIfCan: Fraction Polynomial S -> Union(Float,"failed")
      ++ numericIfCan(x) returns a real approximation of x,
      ++ or "failed" if \axiom{x} is not a constant.
    numericIfCan: (Fraction Polynomial S, PositiveInteger) -> Union(Float,"failed")
      ++ numericIfCan(x,n) returns a real approximation of x up to n decimal
      ++ places, or "failed" if \axiom{x} is not a constant.
    complexNumericIfCan: Fraction Polynomial S -> Union(Complex Float,"failed")
      ++ complexNumericIfCan(x) returns a complex approximation of x,
      ++ or "failed" if \axiom{x} is not a constant.
    complexNumericIfCan: (Fraction Polynomial S, PositiveInteger) -> Union(Complex Float,"failed")
      ++ complexNumericIfCan(x, n) returns a complex approximation of x,
      ++ or "failed" if \axiom{x} is not a constant.
    complexNumericIfCan: Fraction Polynomial Complex S -> Union(Complex Float,"failed")
      ++ complexNumericIfCan(x) returns a complex approximation of x,
      ++ or "failed" if \axiom{x} is not a constant.
    complexNumericIfCan: (Fraction Polynomial Complex S, PositiveInteger) -> 
                                                  Union(Complex Float,"failed")
      ++ complexNumericIfCan(x, n) returns a complex approximation of x
      ++ up to n decimal places, or "failed" if \axiom{x} is not a constant.
    if S has OrderedSet then
      numericIfCan: Expression S -> Union(Float,"failed")
        ++ numericIfCan(x) returns a real approximation of x,
        ++ or "failed" if \axiom{x} is not a constant.
      numericIfCan: (Expression S, PositiveInteger) -> Union(Float,"failed")
        ++ numericIfCan(x, n) returns a real approximation of x up to n
        ++ decimal places, or "failed" if \axiom{x} is not a constant.
      complexNumericIfCan: Expression S -> Union(Complex Float,"failed")
        ++ complexNumericIfCan(x) returns a complex approximation of x,
        ++ or "failed" if \axiom{x} is not a constant.
      complexNumericIfCan: (Expression S, PositiveInteger) ->
                                                  Union(Complex Float,"failed")
        ++ complexNumericIfCan(x, n) returns a complex approximation of x
        ++ up to n decimal places, or "failed" if \axiom{x} is not a constant.
      complexNumericIfCan: Expression Complex S -> Union(Complex Float,"failed")
        ++ complexNumericIfCan(x) returns a complex approximation of x,
        ++ or "failed" if \axiom{x} is not a constant.
      complexNumericIfCan: (Expression Complex S, PositiveInteger) ->
                                                   Union(Complex Float,"failed")
        ++ complexNumericIfCan(x, n) returns a complex approximation of x
        ++ up to n decimal places, or "failed" if \axiom{x} is not a constant.
 == add
 
  if S has CommutativeRing then
    complexNumericIfCan(p:Polynomial Complex S) ==
      p' : Union(Complex(S),"failed") := retractIfCan p
      p' case "failed" => "failed"
      complexNumeric(p')

    complexNumericIfCan(p:Polynomial Complex S,n:PositiveInteger) ==
      p' : Union(Complex(S),"failed") := retractIfCan p
      p' case "failed" => "failed"
      complexNumeric(p',n)
 
  if S has Ring then
    numericIfCan(p:Polynomial S) ==
      p' : Union(S,"failed") := retractIfCan p
      p' case "failed" => "failed"
      numeric(p')

    complexNumericIfCan(p:Polynomial S) ==
      p' : Union(S,"failed") := retractIfCan p
      p' case "failed" => "failed"
      complexNumeric(p')
 
    complexNumericIfCan(p:Polynomial S, n:PositiveInteger) ==
      p' : Union(S,"failed") := retractIfCan p
      p' case "failed" => "failed"
      complexNumeric(p', n)
 
    numericIfCan(p:Polynomial S, n:PositiveInteger) ==
      old := digits(n)$Float
      ans := numericIfCan p
      digits(old)$Float
      ans
 
  if S has IntegralDomain then
    numericIfCan(f:Fraction Polynomial S)==
      num := numericIfCan(numer(f))
      num case "failed" => "failed"
      den := numericIfCan(denom f)
      den case "failed" => "failed"
      num/den
 
    complexNumericIfCan(f:Fraction Polynomial S) ==
      num := complexNumericIfCan(numer f)
      num case "failed" => "failed"
      den := complexNumericIfCan(denom f)
      den case "failed" => "failed"
      num/den
 
    complexNumericIfCan(f:Fraction Polynomial S, n:PositiveInteger) ==
      num := complexNumericIfCan(numer f, n)
      num case "failed" => "failed"
      den := complexNumericIfCan(denom f, n)
      den case "failed" => "failed"
      num/den
 
    numericIfCan(f:Fraction Polynomial S, n:PositiveInteger) ==
      old := digits(n)$Float
      ans := numericIfCan f
      digits(old)$Float
      ans

    complexNumericIfCan(f:Fraction Polynomial Complex S) ==
      num := complexNumericIfCan(numer f)
      num case "failed" => "failed"
      den := complexNumericIfCan(denom f)
      den case "failed" => "failed"
      num/den
 
    complexNumericIfCan(f:Fraction Polynomial Complex S, n:PositiveInteger) ==
      num := complexNumericIfCan(numer f, n)
      num case "failed" => "failed"
      den := complexNumericIfCan(denom f, n)
      den case "failed" => "failed"
      num/den
 
    if S has OrderedSet then
      numericIfCan(x:Expression S) ==
        retractIfCan(map(convert, x)$ExpressionFunctions2(S, Float))
 
      --s2cs(u:S):Complex(S) == complex(u,0)

      complexNumericIfCan(x:Expression S) ==
         complexNumericIfCan map(coerce, x)$ExpressionFunctions2(S,Complex S)
 
      numericIfCan(x:Expression S, n:PositiveInteger) ==
        old := digits(n)$Float
        x' : Expression Float := map(convert, x)$ExpressionFunctions2(S, Float)
        ans : Union(Float,"failed") := retractIfCan x'
        digits(old)$Float
        ans
 
      complexNumericIfCan(x:Expression S, n:PositiveInteger) ==
        old := digits(n)$Float
        x' : Expression Complex S := map(coerce, x)$ExpressionFunctions2(S, Complex S)
        ans : Union(Complex Float,"failed") := complexNumericIfCan(x')
        digits(old)$Float
        ans

      if S has RealConstant then
        complexNumericIfCan(x:Expression Complex S) ==
          retractIfCan(map(convert, x)$ExpressionFunctions2(Complex S,Complex Float))
 
        complexNumericIfCan(x:Expression Complex S, n:PositiveInteger) ==
          old := digits(n)$Float
          x' : Expression Complex Float :=
           map(convert, x)$ExpressionFunctions2(Complex S,Complex Float)
          ans : Union(Complex Float,"failed") := retractIfCan x'
          digits(old)$Float
          ans
      else
        convert(x:Complex S):Complex(Float)==map(convert,x)$ComplexFunctions2(S,Float)

        complexNumericIfCan(x:Expression Complex S) ==
          retractIfCan(map(convert, x)$ExpressionFunctions2(Complex S,Complex Float))
 
        complexNumericIfCan(x:Expression Complex S, n:PositiveInteger) ==
          old := digits(n)$Float
          x' : Expression Complex Float :=
           map(convert, x)$ExpressionFunctions2(Complex S,Complex Float)
          ans : Union(Complex Float,"failed") := retractIfCan x'
          digits(old)$Float
          ans
  numeric(s:S) == convert(s)@Float
 
  if S has ConvertibleTo Complex Float then
    complexNumeric(s:S) == convert(s)@Complex(Float)
 
    complexNumeric(s:S, n:PositiveInteger) ==
      old := digits(n)$Float
      ans := complexNumeric s
      digits(old)$Float
      ans
 
  else
    complexNumeric(s:S) == convert(s)@Float :: Complex(Float)
 
    complexNumeric(s:S,n:PositiveInteger) ==
      numeric(s, n)::Complex(Float)

  if S has CommutativeRing then
    complexNumeric(p:Polynomial Complex S) ==
      p' : Union(Complex(S),"failed") := retractIfCan p
      p' case "failed" => 
        error "Cannot compute the numerical value of a non-constant polynomial"
      complexNumeric(p')

    complexNumeric(p:Polynomial Complex S,n:PositiveInteger) ==
      p' : Union(Complex(S),"failed") := retractIfCan p
      p' case "failed" => 
        error "Cannot compute the numerical value of a non-constant polynomial"
      complexNumeric(p',n)

    if S has RealConstant then
      complexNumeric(s:Complex S) == convert(s)$Complex(S)
  
      complexNumeric(s:Complex S, n:PositiveInteger) ==
        old := digits(n)$Float
        ans := complexNumeric s
        digits(old)$Float
        ans

    else if Complex(S) has ConvertibleTo(Complex Float) then
      complexNumeric(s:Complex S) == convert(s)@Complex(Float)
  
      complexNumeric(s:Complex S, n:PositiveInteger) ==
        old := digits(n)$Float
        ans := complexNumeric s
        digits(old)$Float
        ans

    else
      complexNumeric(s:Complex S) ==
        s' : Union(S,"failed") := retractIfCan s
        s' case "failed" =>
          error "Cannot compute the numerical value of a non-constant object"
        complexNumeric(s')
  
      complexNumeric(s:Complex S, n:PositiveInteger) ==
        s' : Union(S,"failed") := retractIfCan s
        s' case "failed" =>
          error "Cannot compute the numerical value of a non-constant object"
        old := digits(n)$Float
        ans := complexNumeric s'
        digits(old)$Float
        ans
 
  numeric(s:S, n:PositiveInteger) ==
    old := digits(n)$Float
    ans := numeric s
    digits(old)$Float
    ans
 
  if S has Ring then
    numeric(p:Polynomial S) ==
      p' : Union(S,"failed") := retractIfCan p
      p' case "failed" => error
       "Can only compute the numerical value of a constant, real-valued polynomial"
      numeric(p')

    complexNumeric(p:Polynomial S) ==
      p' : Union(S,"failed") := retractIfCan p
      p' case "failed" => 
        error "Cannot compute the numerical value of a non-constant polynomial"
      complexNumeric(p')
 
    complexNumeric(p:Polynomial S, n:PositiveInteger) ==
      p' : Union(S,"failed") := retractIfCan p
      p' case "failed" => 
        error "Cannot compute the numerical value of a non-constant polynomial"
      complexNumeric(p', n)
 
    numeric(p:Polynomial S, n:PositiveInteger) ==
      old := digits(n)$Float
      ans := numeric p
      digits(old)$Float
      ans
 
  if S has IntegralDomain then
    numeric(f:Fraction Polynomial S)==
        numeric(numer(f)) / numeric(denom f)
 
    complexNumeric(f:Fraction Polynomial S) ==
      complexNumeric(numer f)/complexNumeric(denom f)
 
    complexNumeric(f:Fraction Polynomial S, n:PositiveInteger) ==
      complexNumeric(numer f, n)/complexNumeric(denom f, n)
 
    numeric(f:Fraction Polynomial S, n:PositiveInteger) ==
      old := digits(n)$Float
      ans := numeric f
      digits(old)$Float
      ans

    complexNumeric(f:Fraction Polynomial Complex S) ==
      complexNumeric(numer f)/complexNumeric(denom f)
 
    complexNumeric(f:Fraction Polynomial Complex S, n:PositiveInteger) ==
      complexNumeric(numer f, n)/complexNumeric(denom f, n)
 
    if S has OrderedSet then
      numeric(x:Expression S) ==
        x' : Union(Float,"failed") := 
         retractIfCan(map(convert, x)$ExpressionFunctions2(S, Float))
        x' case "failed" => error
         "Can only compute the numerical value of a constant, real-valued Expression"
        x'
 
      complexNumeric(x:Expression S) ==
        x' : Union(Complex Float,"failed") := retractIfCan(
         map(complexNumeric, x)$ExpressionFunctions2(S,Complex Float))
        x' case "failed" =>
         error "Cannot compute the numerical value of a non-constant expression"
        x'
 
      numeric(x:Expression S, n:PositiveInteger) ==
        old := digits(n)$Float
        x' : Expression Float := map(convert, x)$ExpressionFunctions2(S, Float)
        ans : Union(Float,"failed") := retractIfCan x'
        digits(old)$Float
        ans case "failed" => error
         "Can only compute the numerical value of a constant, real-valued Expression"
        ans
 
      complexNumeric(x:Expression S, n:PositiveInteger) ==
        old := digits(n)$Float
        x' : Expression Complex Float :=
         map(complexNumeric, x)$ExpressionFunctions2(S,Complex Float)
        ans : Union(Complex Float,"failed") := retractIfCan x'
        digits(old)$Float
        ans case "failed" =>
         error "Cannot compute the numerical value of a non-constant expression"
        ans

      complexNumeric(x:Expression Complex S) ==
        x' : Union(Complex Float,"failed") := retractIfCan(
         map(complexNumeric, x)$ExpressionFunctions2(Complex S,Complex Float))
        x' case "failed" =>
         error "Cannot compute the numerical value of a non-constant expression"
        x'
 
      complexNumeric(x:Expression Complex S, n:PositiveInteger) ==
        old := digits(n)$Float
        x' : Expression Complex Float :=
         map(complexNumeric, x)$ExpressionFunctions2(Complex S,Complex Float)
        ans : Union(Complex Float,"failed") := retractIfCan x'
        digits(old)$Float
        ans case "failed" =>
         error "Cannot compute the numerical value of a non-constant expression"
        ans

)abbrev package DRAWHACK DrawNumericHack
++ Author: Manuel Bronstein
++ Date Created: 21 Feb 1990
++ Date Last Updated: 21 Feb 1990
++ Basic Operations: coerce
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: 
++ References:
++ Description: Hack for the draw interface. DrawNumericHack provides
++ a "coercion" from something of the form \spad{x = a..b} where \spad{a} 
++ and b are
++ formal expressions to a binding of the form \spad{x = c..d} where c and d
++ are the numerical values of \spad{a} and b. This "coercion" fails if
++ \spad{a} and b contains symbolic variables, but is meant for expressions
++ involving %pi.  
++ NOTE:  This is meant for internal use only.
 
DrawNumericHack(R:Join(OrderedSet,IntegralDomain,ConvertibleTo Float)):
 with coerce: SegmentBinding Expression R -> SegmentBinding Float
        ++ coerce(x = a..b) returns \spad{x = c..d} where c and d are the
        ++ numerical values of \spad{a} and b.
  == add
   coerce s ==
     map(numeric$Numeric(R),s)$SegmentBindingFunctions2(Expression R, Float)

