--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 EXPR2UPS ExpressionToUnivariatePowerSeries
++ Author: Clifton J. Williamson
++ Date Created: 9 May 1989
++ Date Last Updated: 20 September 1993
++ Basic Operations: taylor, laurent, puiseux, series
++ Related Domains: UnivariateTaylorSeries, UnivariateLaurentSeries,
++   UnivariatePuiseuxSeries, Expression
++ Also See: FunctionSpaceToUnivariatePowerSeries
++ AMS Classifications:
++ Keywords: Taylor series, Laurent series, Puiseux series
++ Examples:
++ References:
++ Description:
++   This package provides functions to convert functional expressions
++   to power series.
ExpressionToUnivariatePowerSeries(R,FE): Exports == Implementation where
  R  : Join(GcdDomain,RetractableTo Integer,_
            LinearlyExplicitRingOver Integer)
  FE : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_
            FunctionSpace R)

  EQ     ==> Equation
  I      ==> Integer
  NNI    ==> NonNegativeInteger
  RN     ==> Fraction Integer
  SY     ==> Symbol
  UTS    ==> UnivariateTaylorSeries
  ULS    ==> UnivariateLaurentSeries
  UPXS   ==> UnivariatePuiseuxSeries
  GSER   ==> GeneralUnivariatePowerSeries
  EFULS  ==> ElementaryFunctionsUnivariateLaurentSeries
  EFUPXS ==> ElementaryFunctionsUnivariatePuiseuxSeries
  FS2UPS ==> FunctionSpaceToUnivariatePowerSeries
  Prob   ==> Record(func:String,prob:String)
  ANY1   ==> AnyFunctions1

  Exports ==> with
    taylor: SY -> Any
      ++ \spad{taylor(x)} returns x viewed as a Taylor series.
    taylor: FE -> Any
      ++ \spad{taylor(f)} returns a Taylor expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable.
    taylor: (FE,NNI) -> Any
      ++ \spad{taylor(f,n)} returns a Taylor expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable and terms will be computed
      ++ up to order at least n.
    taylor: (FE,EQ FE) -> Any
      ++ \spad{taylor(f,x = a)} expands the expression f as a Taylor series
      ++ in powers of \spad{(x - a)}.
    taylor: (FE,EQ FE,NNI) -> Any
      ++ \spad{taylor(f,x = a)} expands the expression f as a Taylor series
      ++ in powers of \spad{(x - a)}; terms will be computed up to order
      ++ at least n.

    laurent: SY -> Any
      ++ \spad{laurent(x)} returns x viewed as a Laurent series.
    laurent: FE -> Any
      ++ \spad{laurent(f)} returns a Laurent expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable.
    laurent: (FE,I) -> Any
      ++ \spad{laurent(f,n)} returns a Laurent expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable and terms will be computed
      ++ up to order at least n.
    laurent: (FE,EQ FE) -> Any
      ++ \spad{laurent(f,x = a)} expands the expression f as a Laurent series
      ++ in powers of \spad{(x - a)}.
    laurent: (FE,EQ FE,I) -> Any
      ++ \spad{laurent(f,x = a,n)} expands the expression f as a Laurent
      ++ series in powers of \spad{(x - a)}; terms will be computed up to order
      ++ at least n.
    puiseux: SY -> Any
      ++ \spad{puiseux(x)} returns x viewed as a Puiseux series.
    puiseux: FE -> Any
      ++ \spad{puiseux(f)} returns a Puiseux expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable.
    puiseux: (FE,RN) -> Any
      ++ \spad{puiseux(f,n)} returns a Puiseux expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable and terms will be computed
      ++ up to order at least n.
    puiseux: (FE,EQ FE) -> Any
      ++ \spad{puiseux(f,x = a)} expands the expression f as a Puiseux series
      ++ in powers of \spad{(x - a)}.
    puiseux: (FE,EQ FE,RN) -> Any
      ++ \spad{puiseux(f,x = a,n)} expands the expression f as a Puiseux
      ++ series in powers of \spad{(x - a)}; terms will be computed up to order
      ++ at least n.

    series: SY -> Any
      ++ \spad{series(x)} returns x viewed as a series.
    series: FE -> Any
      ++ \spad{series(f)} returns a series expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable.
    series: (FE,RN) -> Any
      ++ \spad{series(f,n)} returns a series expansion of the expression f.
      ++ Note: f should have only one variable; the series will be
      ++ expanded in powers of that variable and terms will be computed
      ++ up to order at least n.
    series: (FE,EQ FE) -> Any
      ++ \spad{series(f,x = a)} expands the expression f as a series
      ++ in powers of (x - a).
    series: (FE,EQ FE,RN) -> Any
      ++ \spad{series(f,x = a,n)} expands the expression f as a series
      ++ in powers of (x - a); terms will be computed up to order
      ++ at least n.

  Implementation ==> add
    performSubstitution: (FE,SY,FE) -> FE
    performSubstitution(fcn,x,a) ==
      zero? a => fcn
      xFE := x :: FE
      eval(fcn,xFE = xFE + a)

    iTaylor: (FE,SY,FE) -> Any
    iTaylor(fcn,x,a) ==
      pack := FS2UPS(R,FE,I,ULS(FE,x,a),_
                     EFULS(FE,UTS(FE,x,a),ULS(FE,x,a)),x)
      ans := exprToUPS(fcn,false,"just do it")$pack
      ans case %problem =>
        ans.%problem.prob = "essential singularity" =>
          error "No Taylor expansion: essential singularity"
        ans.%problem.func = "log" =>
          error "No Taylor expansion: logarithmic singularity"
        ans.%problem.func = "nth root" =>
          error "No Taylor expansion: fractional powers in expansion"
        error "No Taylor expansion"
      uls := ans.%series
      (uts := taylorIfCan uls) case "failed" =>
        error "No Taylor expansion: pole"
      any1 := ANY1(UTS(FE,x,a))
      coerce(uts :: UTS(FE,x,a))$any1

    taylor(x:SY) ==
      uts := UTS(FE,x,0$FE); any1 := ANY1(uts)
      coerce(monomial(1,1)$uts)$any1

    taylor(fcn:FE) ==
      null(vars := variables fcn) =>
        error "taylor: expression has no variables"
      not null rest vars =>
        error "taylor: expression has more than one variable"
      taylor(fcn,(first(vars) :: FE) = 0)

    taylor(fcn:FE,n:NNI) ==
      null(vars := variables fcn) =>
        error "taylor: expression has no variables"
      not null rest vars =>
        error "taylor: expression has more than one variable"
      x := first vars
      uts := UTS(FE,x,0$FE); any1 := ANY1(uts)
      series := retract(taylor(fcn,(x :: FE) = 0))$any1
      coerce(extend(series,n))$any1

    taylor(fcn:FE,eq:EQ FE) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      iTaylor(performSubstitution(fcn,x,a),x,a)

    taylor(fcn,eq,n) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      any1 := ANY1(UTS(FE,x,a))
      series := retract(iTaylor(performSubstitution(fcn,x,a),x,a))$any1
      coerce(extend(series,n))$any1

    iLaurent: (FE,SY,FE) -> Any
    iLaurent(fcn,x,a) ==
      pack := FS2UPS(R,FE,I,ULS(FE,x,a),_
                     EFULS(FE,UTS(FE,x,a),ULS(FE,x,a)),x)
      ans := exprToUPS(fcn,false,"just do it")$pack
      ans case %problem =>
        ans.%problem.prob = "essential singularity" =>
          error "No Laurent expansion: essential singularity"
        ans.%problem.func = "log" =>
          error "No Laurent expansion: logarithmic singularity"
        ans.%problem.func = "nth root" =>
          error "No Laurent expansion: fractional powers in expansion"
        error "No Laurent expansion"
      any1 := ANY1(ULS(FE,x,a))
      coerce(ans.%series)$any1

    laurent(x:SY) ==
      uls := ULS(FE,x,0$FE); any1 := ANY1(uls)
      coerce(monomial(1,1)$uls)$any1

    laurent(fcn:FE) ==
      null(vars := variables fcn) =>
        error "laurent: expression has no variables"
      not null rest vars =>
        error "laurent: expression has more than one variable"
      laurent(fcn,(first(vars) :: FE) = 0)

    laurent(fcn:FE,n:I) ==
      null(vars := variables fcn) =>
        error "laurent: expression has no variables"
      not null rest vars =>
        error "laurent: expression has more than one variable"
      x := first vars
      uls := ULS(FE,x,0$FE); any1 := ANY1(uls)
      series := retract(laurent(fcn,(x :: FE) = 0))$any1
      coerce(extend(series,n))$any1

    laurent(fcn:FE,eq:EQ FE) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      iLaurent(performSubstitution(fcn,x,a),x,a)

    laurent(fcn,eq,n) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      any1 := ANY1(ULS(FE,x,a))
      series := retract(iLaurent(performSubstitution(fcn,x,a),x,a))$any1
      coerce(extend(series,n))$any1

    iPuiseux: (FE,SY,FE) -> Any
    iPuiseux(fcn,x,a) ==
      pack := FS2UPS(R,FE,RN,UPXS(FE,x,a),_
                     EFUPXS(FE,ULS(FE,x,a),UPXS(FE,x,a),_
                     EFULS(FE,UTS(FE,x,a),ULS(FE,x,a))),x)
      ans := exprToUPS(fcn,false,"just do it")$pack
      ans case %problem =>
        ans.%problem.prob = "essential singularity" =>
          error "No Puiseux expansion: essential singularity"
        ans.%problem.func = "log" =>
          error "No Puiseux expansion: logarithmic singularity"
        error "No Puiseux expansion"
      any1 := ANY1(UPXS(FE,x,a))
      coerce(ans.%series)$any1

    puiseux(x:SY) ==
      upxs := UPXS(FE,x,0$FE); any1 := ANY1(upxs)
      coerce(monomial(1,1)$upxs)$any1

    puiseux(fcn:FE) ==
      null(vars := variables fcn) =>
        error "puiseux: expression has no variables"
      not null rest vars =>
        error "puiseux: expression has more than one variable"
      puiseux(fcn,(first(vars) :: FE) = 0)

    puiseux(fcn:FE,n:RN) ==
      null(vars := variables fcn) =>
        error "puiseux: expression has no variables"
      not null rest vars =>
        error "puiseux: expression has more than one variable"
      x := first vars
      upxs := UPXS(FE,x,0$FE); any1 := ANY1(upxs)
      series := retract(puiseux(fcn,(x :: FE) = 0))$any1
      coerce(extend(series,n))$any1

    puiseux(fcn:FE,eq:EQ FE) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      iPuiseux(performSubstitution(fcn,x,a),x,a)

    puiseux(fcn,eq,n) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      any1 := ANY1(UPXS(FE,x,a))
      series := retract(iPuiseux(performSubstitution(fcn,x,a),x,a))$any1
      coerce(extend(series,n))$any1

    iSeries: (FE,SY,FE) -> Any
    iSeries(fcn,x,a) ==
      pack := FS2UPS(R,FE,RN,UPXS(FE,x,a), _
                     EFUPXS(FE,ULS(FE,x,a),UPXS(FE,x,a), _
                     EFULS(FE,UTS(FE,x,a),ULS(FE,x,a))),x)
      ans := exprToUPS(fcn,false,"just do it")$pack
      ans case %problem =>
        ansG := exprToGenUPS(fcn,false,"just do it")$pack
        ansG case %problem =>
          ansG.%problem.prob = "essential singularity" =>
            error "No series expansion: essential singularity"
          error "No series expansion"
        anyone := ANY1(GSER(FE,x,a))
        coerce((ansG.%series) :: GSER(FE,x,a))$anyone
      any1 := ANY1(UPXS(FE,x,a))
      coerce(ans.%series)$any1

    series(x:SY) ==
      upxs := UPXS(FE,x,0$FE); any1 := ANY1(upxs)
      coerce(monomial(1,1)$upxs)$any1

    series(fcn:FE) ==
      null(vars := variables fcn) =>
        error "series: expression has no variables"
      not null rest vars =>
        error "series: expression has more than one variable"
      series(fcn,(first(vars) :: FE) = 0)

    series(fcn:FE,n:RN) ==
      null(vars := variables fcn) =>
        error "series: expression has no variables"
      not null rest vars =>
        error "series: expression has more than one variable"
      x := first vars
      upxs := UPXS(FE,x,0$FE); any1 := ANY1(upxs)
      series := retract(series(fcn,(x :: FE) = 0))$any1
      coerce(extend(series,n))$any1

    series(fcn:FE,eq:EQ FE) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      iSeries(performSubstitution(fcn,x,a),x,a)

    series(fcn,eq,n) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      any1 := ANY1(UPXS(FE,x,a))
      series := retract(iSeries(performSubstitution(fcn,x,a),x,a))$any1
      coerce(extend(series,n))$any1

