--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 domain SULS SparseUnivariateLaurentSeries
++ Author: Clifton J. Williamson
++ Date Created: 11 November 1994
++ Date Last Updated: June 18, 2010
++ Basic Operations:
++ Related Domains: InnerSparseUnivariatePowerSeries,
++   SparseUnivariateTaylorSeries, SparseUnivariatePuiseuxSeries
++ Also See:
++ AMS Classifications:
++ Keywords: sparse, series
++ Examples:
++ References:
++ Description: Sparse Laurent series in one variable
++   \spadtype{SparseUnivariateLaurentSeries} is a domain representing Laurent
++   series in one variable with coefficients in an arbitrary ring.  The
++   parameters of the type specify the coefficient ring, the power series
++   variable, and the center of the power series expansion.  For example,
++   \spad{SparseUnivariateLaurentSeries(Integer,x,3)} represents Laurent
++   series in \spad{(x - 3)} with integer coefficients.
SparseUnivariateLaurentSeries(Coef,var,cen): Exports == Implementation where
  Coef : Ring
  var  : Symbol
  cen  : Coef
  I     ==> Integer
  NNI   ==> NonNegativeInteger
  OUT   ==> OutputForm
  P     ==> Polynomial Coef
  RF    ==> Fraction Polynomial Coef
  RN    ==> Fraction Integer
  S     ==> String
  SUTS  ==> SparseUnivariateTaylorSeries(Coef,var,cen)
  EFULS ==> ElementaryFunctionsUnivariateLaurentSeries(Coef,SUTS,%)

  Exports ==> Join(UnivariateLaurentSeriesConstructorCategory(Coef,SUTS),_
                   PartialDifferentialDomain(%,Variable var)) with
    coerce: Variable(var) -> %
      ++ \spad{coerce(var)} converts the series variable \spad{var} into a
      ++ Laurent series.
    if Coef has Algebra Fraction Integer then
      integrate: (%,Variable(var)) -> %
        ++ \spad{integrate(f(x))} returns an anti-derivative of the power
        ++ series \spad{f(x)} with constant coefficient 0.
        ++ We may integrate a series when we can divide coefficients
        ++ by integers.

  Implementation ==> InnerSparseUnivariatePowerSeries(Coef) add

    Rep := InnerSparseUnivariatePowerSeries(Coef)

    variable x == var
    center   x == cen

    coerce(v: Variable(var)) ==
      zero? cen => monomial(1,1)
      monomial(1,1) + monomial(cen,0)

    pole? x == negative? order(x,0)

--% operations with Taylor series

    coerce(uts:SUTS) == uts pretend %

    taylorIfCan uls ==
      pole? uls => "failed"
      uls pretend SUTS

    taylor uls ==
      (uts := taylorIfCan uls) case "failed" =>
        error "taylor: Laurent series has a pole"
      uts :: SUTS

    retractIfCan(x:%):Union(SUTS,"failed") == taylorIfCan x

    laurent(n,uts) == monomial(1,n) * (uts :: %)

    removeZeroes uls    == uls
    removeZeroes(n,uls) == uls

    taylorRep uls == taylor(monomial(1,-order(uls,0)) * uls)
    degree uls    == order(uls,0)

    numer uls == taylorRep uls
    denom uls == monomial(1,(-order(uls,0)) :: NNI)$SUTS

    (uts:SUTS) * (uls:%) == (uts :: %) * uls
    (uls:%) * (uts:SUTS) == uls * (uts :: %)

    if Coef has Field then
      (uts1:SUTS) / (uts2:SUTS) == (uts1 :: %) / (uts2 :: %)

    recip(uls) == iExquo(1,uls,false)

    if Coef has IntegralDomain then
      uls1 exquo uls2 == iExquo(uls1,uls2,false)

    if Coef has Field then
      uls1:% / uls2:% ==
        (q := uls1 exquo uls2) case "failed" =>
          error "quotient cannot be computed"
        q :: %

    differentiate(uls:%,v:Variable(var)) == differentiate uls

    elt(uls1:%,uls2:%) ==
      order(uls2,1) < 1 =>
        error "elt: second argument must have positive order"
      negative?(ord := order(uls1,0)) =>
        (recipr := recip uls2) case "failed" =>
          error "elt: second argument not invertible"
        uls3 := uls1 * monomial(1,-ord)
        iCompose(uls3,uls2) * (recipr :: %) ** ((-ord) :: NNI)
      iCompose(uls1,uls2)

    if Coef has IntegralDomain then
      rationalFunction(uls,n) ==
        zero?(e := order(uls,0)) =>
          negative? n => 0
          polynomial(taylor uls,n :: NNI) :: RF
        negative?(m := n - e) => 0
        poly := polynomial(taylor(monomial(1,-e) * uls),m :: NNI) :: RF
        v := variable(uls) :: RF; c := center(uls) :: P :: RF
        poly / (v - c) ** ((-e) :: NNI)

      rationalFunction(uls,n1,n2) == rationalFunction(truncate(uls,n1,n2),n2)

    if Coef has Algebra Fraction Integer then

      integrate uls ==
        zero? coefficient(uls,-1) =>
          error "integrate: series has term of order -1"
        integrate(uls)$Rep

      integrate(uls:%,v:Variable(var)) == integrate uls

      (uls1:%) ** (uls2:%) == exp(log(uls1) * uls2)

      exp uls   == exp(uls)$EFULS
      log uls   == log(uls)$EFULS
      sin uls   == sin(uls)$EFULS
      cos uls   == cos(uls)$EFULS
      tan uls   == tan(uls)$EFULS
      cot uls   == cot(uls)$EFULS
      sec uls   == sec(uls)$EFULS
      csc uls   == csc(uls)$EFULS
      asin uls  == asin(uls)$EFULS
      acos uls  == acos(uls)$EFULS
      atan uls  == atan(uls)$EFULS
      acot uls  == acot(uls)$EFULS
      asec uls  == asec(uls)$EFULS
      acsc uls  == acsc(uls)$EFULS
      sinh uls  == sinh(uls)$EFULS
      cosh uls  == cosh(uls)$EFULS
      tanh uls  == tanh(uls)$EFULS
      coth uls  == coth(uls)$EFULS
      sech uls  == sech(uls)$EFULS
      csch uls  == csch(uls)$EFULS
      asinh uls == asinh(uls)$EFULS
      acosh uls == acosh(uls)$EFULS
      atanh uls == atanh(uls)$EFULS
      acoth uls == acoth(uls)$EFULS
      asech uls == asech(uls)$EFULS
      acsch uls == acsch(uls)$EFULS

      if Coef has CommutativeRing then

        (uls:%) ** (r:RN) == cRationalPower(uls,r)

      else

        (uls:%) ** (r:RN) ==
          negative?(ord0 := order(uls,0)) =>
            order := ord0 :: I
            (n := order exquo denom(r)) case "failed" =>
              error "**: rational power does not exist"
            uts := retract(uls * monomial(1,-order))@SUTS
            utsPow := (uts ** r) :: %
            monomial(1,(n :: I) * numer(r)) * utsPow
          uts := retract(uls)@SUTS
          (uts ** r) :: %

--% OutputForms

    coerce(uls:%): OUT ==
      st := getStream uls
      if not(explicitlyEmpty? st or explicitEntries? st) _
        and (nx := retractIfCan(deref getRef uls))@Union(I,"failed") case I then
        count : NNI := _$streamCount$Lisp
        degr := min(count,(nx :: I) + count + 1)
        extend(uls,degr)
      seriesToOutputForm(st,getRef uls,variable uls,center uls,1)

