--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 EFUPXS ElementaryFunctionsUnivariatePuiseuxSeries
++ This package provides elementary functions on Puiseux series.
++ Author: Clifton J. Williamson
++ Date Created: 20 February 1990
++ Date Last Updated: 20 February 1990
++ Keywords: elementary function, Laurent series
++ Examples:
++ References:
ElementaryFunctionsUnivariatePuiseuxSeries(Coef,ULS,UPXS,EFULS):_
 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
  ULS    : UnivariateLaurentSeriesCategory Coef
  UPXS   : UnivariatePuiseuxSeriesConstructorCategory(Coef,ULS)
  EFULS  : PartialTranscendentalFunctions(ULS)
  I    ==> Integer
  NNI  ==> NonNegativeInteger
  RN   ==> Fraction Integer
 
  Exports ==> PartialTranscendentalFunctions(UPXS) with
 
    if Coef has Field then
      **: (UPXS,RN) -> UPXS
        ++ z ** r raises a Puiseaux series z to a rational power r
 
--% Exponentials and Logarithms
 
    exp: UPXS -> UPXS
      ++ exp(z) returns the exponential of a Puiseux series z.
    log: UPXS -> UPXS
      ++ log(z) returns the logarithm of a Puiseux series z.
 
--% TrigonometricFunctionCategory
 
    sin: UPXS -> UPXS
      ++ sin(z) returns the sine of a Puiseux series z.
    cos: UPXS -> UPXS
      ++ cos(z) returns the cosine of a Puiseux series z.
    tan: UPXS -> UPXS
      ++ tan(z) returns the tangent of a Puiseux series z.
    cot: UPXS -> UPXS
      ++ cot(z) returns the cotangent of a Puiseux series z.
    sec: UPXS -> UPXS
      ++ sec(z) returns the secant of a Puiseux series z.
    csc: UPXS -> UPXS
      ++ csc(z) returns the cosecant of a Puiseux series z.
 
--% ArcTrigonometricFunctionCategory
 
    asin: UPXS -> UPXS
      ++ asin(z) returns the arc-sine of a Puiseux series z.
    acos: UPXS -> UPXS
      ++ acos(z) returns the arc-cosine of a Puiseux series z.
    atan: UPXS -> UPXS
      ++ atan(z) returns the arc-tangent of a Puiseux series z.
    acot: UPXS -> UPXS
      ++ acot(z) returns the arc-cotangent of a Puiseux series z.
    asec: UPXS -> UPXS
      ++ asec(z) returns the arc-secant of a Puiseux series z.
    acsc: UPXS -> UPXS
      ++ acsc(z) returns the arc-cosecant of a Puiseux series z.
 
--% HyperbolicFunctionCategory
 
    sinh: UPXS -> UPXS
      ++ sinh(z) returns the hyperbolic sine of a Puiseux series z.
    cosh: UPXS -> UPXS
      ++ cosh(z) returns the hyperbolic cosine of a Puiseux series z.
    tanh: UPXS -> UPXS
      ++ tanh(z) returns the hyperbolic tangent of a Puiseux series z.
    coth: UPXS -> UPXS
      ++ coth(z) returns the hyperbolic cotangent of a Puiseux series z.
    sech: UPXS -> UPXS
      ++ sech(z) returns the hyperbolic secant of a Puiseux series z.
    csch: UPXS -> UPXS
      ++ csch(z) returns the hyperbolic cosecant of a Puiseux series z.
 
--% ArcHyperbolicFunctionCategory
 
    asinh: UPXS -> UPXS
      ++ asinh(z) returns the inverse hyperbolic sine of a Puiseux series z.
    acosh: UPXS -> UPXS
      ++ acosh(z) returns the inverse hyperbolic cosine of a Puiseux series z.
    atanh: UPXS -> UPXS
      ++ atanh(z) returns the inverse hyperbolic tangent of a Puiseux series z.
    acoth: UPXS -> UPXS
      ++ acoth(z) returns the inverse hyperbolic cotangent of a Puiseux series z.
    asech: UPXS -> UPXS
      ++ asech(z) returns the inverse hyperbolic secant of a Puiseux series z.
    acsch: UPXS -> UPXS
      ++ acsch(z) returns the inverse hyperbolic cosecant of a Puiseux series z.
 
  Implementation ==> add

    TRANSFCN : Boolean := Coef has TranscendentalFunctionCategory
 
--% roots
 
    nthRootIfCan(upxs,n) ==
      one? n => upxs
      r := rationalPower upxs; uls := laurentRep upxs
      deg := degree uls
      if zero?(coef := coefficient(uls,deg)) then
        deg := order(uls,deg + 1000)
        zero?(coef := coefficient(uls,deg)) =>
          error "root of series with many leading zero coefficients"
      uls := uls * monomial(1,-deg)$ULS
      (ulsRoot := nthRootIfCan(uls,n)) case "failed" => "failed"
      puiseux(r,ulsRoot :: ULS) * monomial(1,deg * r * inv(n :: RN))
 
    if Coef has Field then
       (upxs:UPXS) ** (q:RN) ==
         num := numer q; den := denom q
         one? den => upxs ** num
         r := rationalPower upxs; uls := laurentRep upxs
         deg := degree uls
         if zero?(coef := coefficient(uls,deg)) then
           deg := order(uls,deg + 1000)
           zero?(coef := coefficient(uls,deg)) =>
             error "power of series with many leading zero coefficients"
         ulsPow := (uls * monomial(1,-deg)$ULS) ** q
         puiseux(r,ulsPow) * monomial(1,deg*q*r)
 
--% transcendental functions
 
    applyIfCan: (ULS -> Union(ULS,"failed"),UPXS) -> Union(UPXS,"failed")
    applyIfCan(fcn,upxs) ==
      uls := fcn laurentRep upxs
      uls case "failed" => "failed"
      puiseux(rationalPower upxs,uls :: ULS)
 
    expIfCan   upxs == applyIfCan(expIfCan,upxs)
    logIfCan   upxs == applyIfCan(logIfCan,upxs)
    sinIfCan   upxs == applyIfCan(sinIfCan,upxs)
    cosIfCan   upxs == applyIfCan(cosIfCan,upxs)
    tanIfCan   upxs == applyIfCan(tanIfCan,upxs)
    cotIfCan   upxs == applyIfCan(cotIfCan,upxs)
    secIfCan   upxs == applyIfCan(secIfCan,upxs)
    cscIfCan   upxs == applyIfCan(cscIfCan,upxs)
    atanIfCan  upxs == applyIfCan(atanIfCan,upxs)
    acotIfCan  upxs == applyIfCan(acotIfCan,upxs)
    sinhIfCan  upxs == applyIfCan(sinhIfCan,upxs)
    coshIfCan  upxs == applyIfCan(coshIfCan,upxs)
    tanhIfCan  upxs == applyIfCan(tanhIfCan,upxs)
    cothIfCan  upxs == applyIfCan(cothIfCan,upxs)
    sechIfCan  upxs == applyIfCan(sechIfCan,upxs)
    cschIfCan  upxs == applyIfCan(cschIfCan,upxs)
    atanhIfCan upxs == applyIfCan(atanhIfCan,upxs)
    acothIfCan upxs == applyIfCan(acothIfCan,upxs)

    asinIfCan upxs ==
      negative? order(upxs,0) => "failed"
      (coef := coefficient(upxs,0)) = 0 =>
        integrate((1 - upxs*upxs)**(-1/2) * (differentiate upxs))
      TRANSFCN =>
        cc := asin(coef) :: UPXS
        cc + integrate((1 - upxs*upxs)**(-1/2) * (differentiate upxs))
      "failed"

    acosIfCan upxs ==
      negative? order(upxs,0) => "failed"
      TRANSFCN =>
        cc := acos(coefficient(upxs,0)) :: UPXS
        cc + integrate(-(1 - upxs*upxs)**(-1/2) * (differentiate upxs))
      "failed"

    asecIfCan upxs ==
      negative? order(upxs,0) => "failed"
      TRANSFCN =>
        cc := asec(coefficient(upxs,0)) :: UPXS
        f := (upxs*upxs - 1)**(-1/2) * (differentiate upxs)
        (rec := recip upxs) case "failed" => "failed"
        cc + integrate(f * (rec :: UPXS))
      "failed"

    acscIfCan upxs ==
      negative? order(upxs,0) => "failed"
      TRANSFCN =>
        cc := acsc(coefficient(upxs,0)) :: UPXS
        f := -(upxs*upxs - 1)**(-1/2) * (differentiate upxs)
        (rec := recip upxs) case "failed" => "failed"
        cc + integrate(f * (rec :: UPXS))
      "failed"

    asinhIfCan upxs ==
      negative? order(upxs,0) => "failed"
      TRANSFCN or (coefficient(upxs,0) = 0) =>
        log(upxs + (1 + upxs*upxs)**(1/2))
      "failed"

    acoshIfCan upxs ==
      TRANSFCN =>
        negative? order(upxs,0) => "failed"
        log(upxs + (upxs*upxs - 1)**(1/2))
      "failed"

    asechIfCan upxs ==
      TRANSFCN =>
        negative? order(upxs,0) => "failed"
        (rec := recip upxs) case "failed" => "failed"
        log((1 + (1 - upxs*upxs)*(1/2)) * (rec :: UPXS))
      "failed"

    acschIfCan upxs ==
      TRANSFCN =>
        negative? order(upxs,0) => "failed"
        (rec := recip upxs) case "failed" => "failed"
        log((1 + (1 + upxs*upxs)*(1/2)) * (rec :: UPXS))
      "failed"
 
    applyOrError:(UPXS -> Union(UPXS,"failed"),String,UPXS) -> UPXS
    applyOrError(fcn,name,upxs) ==
      ans := fcn upxs
      ans case "failed" =>
        error concat(name," of function with singularity")
      ans :: UPXS
 
    exp upxs   == applyOrError(expIfCan,"exp",upxs)
    log upxs   == applyOrError(logIfCan,"log",upxs)
    sin upxs   == applyOrError(sinIfCan,"sin",upxs)
    cos upxs   == applyOrError(cosIfCan,"cos",upxs)
    tan upxs   == applyOrError(tanIfCan,"tan",upxs)
    cot upxs   == applyOrError(cotIfCan,"cot",upxs)
    sec upxs   == applyOrError(secIfCan,"sec",upxs)
    csc upxs   == applyOrError(cscIfCan,"csc",upxs)
    asin upxs  == applyOrError(asinIfCan,"asin",upxs)
    acos upxs  == applyOrError(acosIfCan,"acos",upxs)
    atan upxs  == applyOrError(atanIfCan,"atan",upxs)
    acot upxs  == applyOrError(acotIfCan,"acot",upxs)
    asec upxs  == applyOrError(asecIfCan,"asec",upxs)
    acsc upxs  == applyOrError(acscIfCan,"acsc",upxs)
    sinh upxs  == applyOrError(sinhIfCan,"sinh",upxs)
    cosh upxs  == applyOrError(coshIfCan,"cosh",upxs)
    tanh upxs  == applyOrError(tanhIfCan,"tanh",upxs)
    coth upxs  == applyOrError(cothIfCan,"coth",upxs)
    sech upxs  == applyOrError(sechIfCan,"sech",upxs)
    csch upxs  == applyOrError(cschIfCan,"csch",upxs)
    asinh upxs == applyOrError(asinhIfCan,"asinh",upxs)
    acosh upxs == applyOrError(acoshIfCan,"acosh",upxs)
    atanh upxs == applyOrError(atanhIfCan,"atanh",upxs)
    acoth upxs == applyOrError(acothIfCan,"acoth",upxs)
    asech upxs == applyOrError(asechIfCan,"asech",upxs)
    acsch upxs == applyOrError(acschIfCan,"acsch",upxs)

