--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 EFULS ElementaryFunctionsUnivariateLaurentSeries
++ This package provides elementary functions on Laurent series.
++ Author: Clifton J. Williamson
++ Date Created: 6 February 1990
++ Date Last Updated: 25 February 1990
++ Keywords: elementary function, Laurent series
++ Examples:
++ References:
ElementaryFunctionsUnivariateLaurentSeries(Coef,UTS,ULS):_
 Exports == Implementation where
  ++ This package provides elementary functions on any Laurent series
  ++ domain over a field which was constructed from a Taylor series
  ++ domain.  These functions are implemented by calling the
  ++ corresponding functions on the Taylor series domain.  We also
  ++ provide 'partial functions' which compute transcendental
  ++ functions of Laurent series when possible and return "failed"
  ++ when this is not possible.
  Coef   : Algebra Fraction Integer
  UTS    : UnivariateTaylorSeriesCategory Coef
  ULS    : UnivariateLaurentSeriesConstructorCategory(Coef,UTS)
  I    ==> Integer
  NNI  ==> NonNegativeInteger
  RN   ==> Fraction Integer
  S    ==> String
  STTF ==> StreamTranscendentalFunctions(Coef)
 
  Exports ==> PartialTranscendentalFunctions(ULS) with
 
    if Coef has Field then
      **: (ULS,RN) -> ULS
        ++ s ** r raises a Laurent series s to a rational power r
 
--% Exponentials and Logarithms
 
    exp: ULS -> ULS
      ++ exp(z) returns the exponential of Laurent series z.
    log: ULS -> ULS
      ++ log(z) returns the logarithm of Laurent series z.
 
--% TrigonometricFunctionCategory
 
    sin: ULS -> ULS
      ++ sin(z) returns the sine of Laurent series z.
    cos: ULS -> ULS
      ++ cos(z) returns the cosine of Laurent series z.
    tan: ULS -> ULS
      ++ tan(z) returns the tangent of Laurent series z.
    cot: ULS -> ULS
      ++ cot(z) returns the cotangent of Laurent series z.
    sec: ULS -> ULS
      ++ sec(z) returns the secant of Laurent series z.
    csc: ULS -> ULS
      ++ csc(z) returns the cosecant of Laurent series z.
 
--% ArcTrigonometricFunctionCategory
 
    asin: ULS -> ULS
      ++ asin(z) returns the arc-sine of Laurent series z.
    acos: ULS -> ULS
      ++ acos(z) returns the arc-cosine of Laurent series z.
    atan: ULS -> ULS
      ++ atan(z) returns the arc-tangent of Laurent series z.
    acot: ULS -> ULS
      ++ acot(z) returns the arc-cotangent of Laurent series z.
    asec: ULS -> ULS
      ++ asec(z) returns the arc-secant of Laurent series z.
    acsc: ULS -> ULS
      ++ acsc(z) returns the arc-cosecant of Laurent series z.
 
--% HyperbolicFunctionCategory
 
    sinh: ULS -> ULS
      ++ sinh(z) returns the hyperbolic sine of Laurent series z.
    cosh: ULS -> ULS
      ++ cosh(z) returns the hyperbolic cosine of Laurent series z.
    tanh: ULS -> ULS
      ++ tanh(z) returns the hyperbolic tangent of Laurent series z.
    coth: ULS -> ULS
      ++ coth(z) returns the hyperbolic cotangent of Laurent series z.
    sech: ULS -> ULS
      ++ sech(z) returns the hyperbolic secant of Laurent series z.
    csch: ULS -> ULS
      ++ csch(z) returns the hyperbolic cosecant of Laurent series z.
 
--% ArcHyperbolicFunctionCategory
 
    asinh: ULS -> ULS
      ++ asinh(z) returns the inverse hyperbolic sine of Laurent series z.
    acosh: ULS -> ULS
      ++ acosh(z) returns the inverse hyperbolic cosine of Laurent series z.
    atanh: ULS -> ULS
      ++ atanh(z) returns the inverse hyperbolic tangent of Laurent series z.
    acoth: ULS -> ULS
      ++ acoth(z) returns the inverse hyperbolic cotangent of Laurent series z.
    asech: ULS -> ULS
      ++ asech(z) returns the inverse hyperbolic secant of Laurent series z.
    acsch: ULS -> ULS
      ++ acsch(z) returns the inverse hyperbolic cosecant of Laurent series z.
 
  Implementation ==> add
 
--% roots
 
    RATPOWERS : Boolean := Coef has "**":(Coef,RN) -> Coef
    TRANSFCN  : Boolean := Coef has TranscendentalFunctionCategory
    RATS      : Boolean := Coef has retractIfCan: Coef -> Union(RN,"failed")
 
    nthRootUTS:(UTS,I) -> Union(UTS,"failed")
    nthRootUTS(uts,n) ==
      -- assumed: n > 1, uts has non-zero constant term
      one? coefficient(uts,0) => uts ** inv(n::RN)
      RATPOWERS => uts ** inv(n::RN)
      "failed"
 
    nthRootIfCan(uls,nn) ==
      (n := nn :: I) < 1 => error "nthRootIfCan: n must be positive"
      n = 1 => uls
      deg := degree uls
      if zero? (coef := coefficient(uls,deg)) then
        uls := removeZeroes(1000,uls); deg := degree uls
        zero? (coef := coefficient(uls,deg)) =>
          error "root of series with many leading zero coefficients"
      (k := deg exquo n) case "failed" => "failed"
      uts := taylor(uls * monomial(1,-deg))
      (root := nthRootUTS(uts,n)) case "failed" => "failed"
      monomial(1,k :: I) * (root :: UTS :: ULS)
 
    if Coef has Field then
       (uls:ULS) ** (r:RN) ==
         num := numer r; den := denom r
         one? den => uls ** num
         deg := degree uls
         if zero? (coef := coefficient(uls,deg)) then
           uls := removeZeroes(1000,uls); deg := degree uls
           zero? (coef := coefficient(uls,deg)) =>
             error "power of series with many leading zero coefficients"
         (k := deg exquo den) case "failed" =>
           error "**: rational power does not exist"
         uts := taylor(uls * monomial(1,-deg)) ** r
         monomial(1,(k :: I) * num) * (uts :: ULS)
 
--% transcendental functions
 
    applyIfCan: (UTS -> UTS,ULS) -> Union(ULS,"failed")
    applyIfCan(fcn,uls) ==
      uts := taylorIfCan uls
      uts case "failed" => "failed"
      fcn(uts :: UTS) :: ULS
 
    expIfCan   uls == applyIfCan(exp,uls)
    sinIfCan   uls == applyIfCan(sin,uls)
    cosIfCan   uls == applyIfCan(cos,uls)
    asinIfCan  uls == applyIfCan(asin,uls)
    acosIfCan  uls == applyIfCan(acos,uls)
    asecIfCan  uls == applyIfCan(asec,uls)
    acscIfCan  uls == applyIfCan(acsc,uls)
    sinhIfCan  uls == applyIfCan(sinh,uls)
    coshIfCan  uls == applyIfCan(cosh,uls)
    asinhIfCan uls == applyIfCan(asinh,uls)
    acoshIfCan uls == applyIfCan(acosh,uls)
    atanhIfCan uls == applyIfCan(atanh,uls)
    acothIfCan uls == applyIfCan(acoth,uls)
    asechIfCan uls == applyIfCan(asech,uls)
    acschIfCan uls == applyIfCan(acsch,uls)
 
    logIfCan uls ==
      uts := taylorIfCan uls
      uts case "failed" => "failed"
      zero? coefficient(ts := uts :: UTS,0) => "failed"
      log(ts) :: ULS
 
    tanIfCan uls ==
      -- don't call 'tan' on a UTS (tan(uls) may have a singularity)
      uts := taylorIfCan uls
      uts case "failed" => "failed"
      sc := sincos(coefficients(uts :: UTS))$STTF
      (cosInv := recip(series(sc.cos) :: ULS)) case "failed" => "failed"
      (series(sc.sin) :: ULS) * (cosInv :: ULS)
 
    cotIfCan uls ==
      -- don't call 'cot' on a UTS (cot(uls) may have a singularity)
      uts := taylorIfCan uls
      uts case "failed" => "failed"
      sc := sincos(coefficients(uts :: UTS))$STTF
      (sinInv := recip(series(sc.sin) :: ULS)) case "failed" => "failed"
      (series(sc.cos) :: ULS) * (sinInv :: ULS)
 
    secIfCan uls ==
      cos := cosIfCan uls
      cos case "failed" => "failed"
      (cosInv := recip(cos :: ULS)) case "failed" => "failed"
      cosInv :: ULS
 
    cscIfCan uls ==
      sin := sinIfCan uls
      sin case "failed" => "failed"
      (sinInv := recip(sin :: ULS)) case "failed" => "failed"
      sinInv :: ULS

    atanIfCan uls ==
      coef := coefficient(uls,0)
      (ord := order(uls,0)) = 0 and coef * coef = -1 => "failed"
      cc : Coef := 
        negative? ord =>
          TRANSFCN =>
            RATS =>
              lc := coefficient(uls,ord)
              (rat := retractIfCan(lc)@Union(RN,"failed")) case "failed" =>
                (1/2) * pi()
              positive?(rat :: RN) => (1/2) * pi()
              (-1/2) * pi()
            (1/2) * pi()
          return "failed"
        coef = 0 => 0
        TRANSFCN => atan coef
        return "failed"
      (z := recip(1 + uls*uls)) case "failed" => "failed"
      (cc :: ULS) + integrate(differentiate(uls) * (z :: ULS))

    acotIfCan uls ==
      coef := coefficient(uls,0)
      (ord := order(uls,0)) = 0 and coef * coef = -1 => "failed"
      cc : Coef := 
        negative? ord =>
          RATS =>
            lc := coefficient(uls,ord)
            (rat := retractIfCan(lc)@Union(RN,"failed")) case "failed" => 0
            positive?(rat :: RN) => 0
            TRANSFCN => pi()
            return "failed"
          0
        TRANSFCN => acot coef
        return "failed"
      (z := recip(1 + uls*uls)) case "failed" => "failed"
      (cc :: ULS) - integrate(differentiate(uls) * (z :: ULS))
 
    tanhIfCan uls ==
      -- don't call 'tanh' on a UTS (tanh(uls) may have a singularity)
      uts := taylorIfCan uls
      uts case "failed" => "failed"
      sc := sinhcosh(coefficients(uts :: UTS))$STTF
      (coshInv := recip(series(sc.cosh) :: ULS)) case "failed" =>
        "failed"
      (series(sc.sinh) :: ULS) * (coshInv :: ULS)
 
    cothIfCan uls ==
      -- don't call 'coth' on a UTS (coth(uls) may have a singularity)
      uts := taylorIfCan uls
      uts case "failed" => "failed"
      sc := sinhcosh(coefficients(uts :: UTS))$STTF
      (sinhInv := recip(series(sc.sinh) :: ULS)) case "failed" =>
        "failed"
      (series(sc.cosh) :: ULS) * (sinhInv :: ULS)
 
    sechIfCan uls ==
      cosh := coshIfCan uls
      cosh case "failed" => "failed"
      (coshInv := recip(cosh :: ULS)) case "failed" => "failed"
      coshInv :: ULS
 
    cschIfCan uls ==
      sinh := sinhIfCan uls
      sinh case "failed" => "failed"
      (sinhInv := recip(sinh :: ULS)) case "failed" => "failed"
      sinhInv :: ULS
 
    applyOrError:(ULS -> Union(ULS,"failed"),S,ULS) -> ULS
    applyOrError(fcn,name,uls) ==
      ans := fcn uls
      ans case "failed" =>
        error concat(name," of function with singularity")
      ans :: ULS
 
    exp uls   == applyOrError(expIfCan,"exp",uls)
    log uls   == applyOrError(logIfCan,"log",uls)
    sin uls   == applyOrError(sinIfCan,"sin",uls)
    cos uls   == applyOrError(cosIfCan,"cos",uls)
    tan uls   == applyOrError(tanIfCan,"tan",uls)
    cot uls   == applyOrError(cotIfCan,"cot",uls)
    sec uls   == applyOrError(secIfCan,"sec",uls)
    csc uls   == applyOrError(cscIfCan,"csc",uls)
    asin uls  == applyOrError(asinIfCan,"asin",uls)
    acos uls  == applyOrError(acosIfCan,"acos",uls)
    asec uls  == applyOrError(asecIfCan,"asec",uls)
    acsc uls  == applyOrError(acscIfCan,"acsc",uls)
    sinh uls  == applyOrError(sinhIfCan,"sinh",uls)
    cosh uls  == applyOrError(coshIfCan,"cosh",uls)
    tanh uls  == applyOrError(tanhIfCan,"tanh",uls)
    coth uls  == applyOrError(cothIfCan,"coth",uls)
    sech uls  == applyOrError(sechIfCan,"sech",uls)
    csch uls  == applyOrError(cschIfCan,"csch",uls)
    asinh uls == applyOrError(asinhIfCan,"asinh",uls)
    acosh uls == applyOrError(acoshIfCan,"acosh",uls)
    atanh uls == applyOrError(atanhIfCan,"atanh",uls)
    acoth uls == applyOrError(acothIfCan,"acoth",uls)
    asech uls == applyOrError(asechIfCan,"asech",uls)
    acsch uls == applyOrError(acschIfCan,"acsch",uls)

    atan uls ==
    -- code is duplicated so that correct error messages will be returned
      coef := coefficient(uls,0)
      (ord := order(uls,0)) = 0 and coef * coef = -1 =>
        error "atan: series expansion has logarithmic term"
      cc : Coef := 
        negative? ord =>
          TRANSFCN =>
            RATS =>
              lc := coefficient(uls,ord)
              (rat := retractIfCan(lc)@Union(RN,"failed")) case "failed" =>
                (1/2) * pi()
              positive?(rat :: RN) => (1/2) * pi()
              (-1/2) * pi()
            (1/2) * pi()
          error "atan: series expansion involves transcendental constants"
        coef = 0 => 0
        TRANSFCN => atan coef
        error "atan: series expansion involves transcendental constants"
      (z := recip(1 + uls*uls)) case "failed" =>
        error "atan: leading coefficient not invertible"
      (cc :: ULS) + integrate(differentiate(uls) * (z :: ULS))

    acot uls ==
    -- code is duplicated so that correct error messages will be returned
      coef := coefficient(uls,0)
      (ord := order(uls,0)) = 0 and coef * coef = -1 =>
        error "acot: series expansion has logarithmic term"
      cc : Coef := 
        negative? ord =>
          RATS =>
            lc := coefficient(uls,ord)
            (rat := retractIfCan(lc)@Union(RN,"failed")) case "failed" => 0
            positive?(rat :: RN) => 0
            TRANSFCN => pi()
            error "acot: series expansion involves transcendental constants"
          0
        TRANSFCN => acot coef
        error "acot: series expansion involves transcendental constants"
      (z := recip(1 + uls*uls)) case "failed" =>
        error "acot: leading coefficient not invertible"
      (cc :: ULS) - integrate(differentiate(uls) * (z :: ULS))

