--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 GENUPS GenerateUnivariatePowerSeries
++ Author: Clifton J. Williamson
++ Date Created: 29 April 1990
++ Date Last Updated: 31 May 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: series, Taylor, Laurent, Puiseux
++ Examples:
++ References:
++ Description:
++   \spadtype{GenerateUnivariatePowerSeries} provides functions that create
++   power series from explicit formulas for their \spad{n}th coefficient.
GenerateUnivariatePowerSeries(R,FE): Exports == Implementation where
  R  : Join(IntegralDomain,RetractableTo Integer,_
            LinearlyExplicitRingOver Integer)
  FE : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_
            FunctionSpace R)
  ANY1 ==> AnyFunctions1
  EQ   ==> Equation
  I    ==> Integer
  NNI  ==> NonNegativeInteger
  RN   ==> Fraction Integer
  SEG  ==> UniversalSegment
  ST   ==> Stream
  SY   ==> Symbol
  UTS  ==> UnivariateTaylorSeries
  ULS  ==> UnivariateLaurentSeries
  UPXS ==> UnivariatePuiseuxSeries
 
  Exports ==> with
    taylor: (I -> FE,EQ FE) -> Any
      ++ \spad{taylor(n +-> a(n),x = a)} returns
      ++ \spad{sum(n = 0..,a(n)*(x-a)**n)}.
    taylor: (FE,SY,EQ FE) -> Any
      ++ \spad{taylor(a(n),n,x = a)} returns \spad{sum(n = 0..,a(n)*(x-a)**n)}.
    taylor: (I -> FE,EQ FE,SEG NNI) -> Any
      ++ \spad{taylor(n +-> a(n),x = a,n0..)} returns
      ++ \spad{sum(n=n0..,a(n)*(x-a)**n)};
      ++ \spad{taylor(n +-> a(n),x = a,n0..n1)} returns
      ++ \spad{sum(n = n0..,a(n)*(x-a)**n)}.
    taylor: (FE,SY,EQ FE,SEG NNI) -> Any
      ++ \spad{taylor(a(n),n,x = a,n0..)} returns
      ++ \spad{sum(n = n0..,a(n)*(x-a)**n)};
      ++ \spad{taylor(a(n),n,x = a,n0..n1)} returns
      ++ \spad{sum(n = n0..,a(n)*(x-a)**n)}.
 
    laurent: (I -> FE,EQ FE,SEG I) -> Any
      ++ \spad{laurent(n +-> a(n),x = a,n0..)} returns
      ++ \spad{sum(n = n0..,a(n) * (x - a)**n)};
      ++ \spad{laurent(n +-> a(n),x = a,n0..n1)} returns
      ++ \spad{sum(n = n0..n1,a(n) * (x - a)**n)}.
    laurent: (FE,SY,EQ FE,SEG I) -> Any
      ++ \spad{laurent(a(n),n,x=a,n0..)} returns
      ++ \spad{sum(n = n0..,a(n) * (x - a)**n)};
      ++ \spad{laurent(a(n),n,x=a,n0..n1)} returns
      ++ \spad{sum(n = n0..n1,a(n) * (x - a)**n)}.
 
    puiseux: (RN -> FE,EQ FE,SEG RN,RN) -> Any
      ++ \spad{puiseux(n +-> a(n),x = a,r0..,r)} returns
      ++ \spad{sum(n = r0,r0 + r,r0 + 2*r..., a(n) * (x - a)**n)};
      ++ \spad{puiseux(n +-> a(n),x = a,r0..r1,r)} returns
      ++ \spad{sum(n = r0 + k*r while n <= r1, a(n) * (x - a)**n)}.
    puiseux: (FE,SY,EQ FE,SEG RN,RN) -> Any
      ++ \spad{puiseux(a(n),n,x = a,r0..,r)} returns
      ++ \spad{sum(n = r0,r0 + r,r0 + 2*r..., a(n) * (x - a)**n)};
      ++ \spad{puiseux(a(n),n,x = a,r0..r1,r)} returns
      ++ \spad{sum(n = r0 + k*r while n <= r1, a(n) * (x - a)**n)}.
 
    series: (I -> FE,EQ FE) -> Any
      ++ \spad{series(n +-> a(n),x = a)} returns
      ++ \spad{sum(n = 0..,a(n)*(x-a)**n)}.
    series: (FE,SY,EQ FE) -> Any
      ++ \spad{series(a(n),n,x = a)} returns
      ++ \spad{sum(n = 0..,a(n)*(x-a)**n)}.
    series: (I -> FE,EQ FE,SEG I) -> Any
      ++ \spad{series(n +-> a(n),x = a,n0..)} returns
      ++ \spad{sum(n = n0..,a(n) * (x - a)**n)};
      ++ \spad{series(n +-> a(n),x = a,n0..n1)} returns
      ++ \spad{sum(n = n0..n1,a(n) * (x - a)**n)}.
    series: (FE,SY,EQ FE,SEG I) -> Any
      ++ \spad{series(a(n),n,x=a,n0..)} returns
      ++ \spad{sum(n = n0..,a(n) * (x - a)**n)};
      ++ \spad{series(a(n),n,x=a,n0..n1)} returns
      ++ \spad{sum(n = n0..n1,a(n) * (x - a)**n)}.
    series: (RN -> FE,EQ FE,SEG RN,RN) -> Any
      ++ \spad{series(n +-> a(n),x = a,r0..,r)} returns
      ++ \spad{sum(n = r0,r0 + r,r0 + 2*r..., a(n) * (x - a)**n)};
      ++ \spad{series(n +-> a(n),x = a,r0..r1,r)} returns
      ++ \spad{sum(n = r0 + k*r while n <= r1, a(n) * (x - a)**n)}.
    series: (FE,SY,EQ FE,SEG RN,RN) -> Any
      ++ \spad{series(a(n),n,x = a,r0..,r)} returns
      ++ \spad{sum(n = r0,r0 + r,r0 + 2*r..., a(n) * (x - a)**n)};
      ++ \spad{series(a(n),n,x = a,r0..r1,r)} returns
      ++ \spad{sum(n = r0 + k*r while n <= r1, a(n) * (x - a)**n)}.
 
  Implementation ==> add
 
    genStream: (I -> FE,I) -> ST FE
    genStream(f,n) == delay concat(f(n),genStream(f,n + 1))
 
    genFiniteStream: (I -> FE,I,I) -> ST FE
    genFiniteStream(f,n,m) == delay
      n > m => empty()
      concat(f(n),genFiniteStream(f,n + 1,m))
 
    taylor(f,eq) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      coerce(series(genStream(f,0))$UTS(FE,x,a))$ANY1(UTS(FE,x,a))
 
    taylor(an:FE,n:SY,eq:EQ FE) ==
      taylor(eval(an,(n :: FE) = (#1 :: FE)),eq)
 
    taylor(f:I -> FE,eq:EQ FE,seg:SEG NNI) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      hasHi seg =>
        n0 := lo seg; n1 := hi seg
        if n1 < n0 then (n0,n1) := (n1,n0)
        uts := series(genFiniteStream(f,n0,n1))$UTS(FE,x,a)
        uts := uts * monomial(1,n0)$UTS(FE,x,a)
        coerce(uts)$ANY1(UTS(FE,x,a))
      n0 := lo seg
      uts := series(genStream(f,n0))$UTS(FE,x,a)
      uts := uts * monomial(1,n0)$UTS(FE,x,a)
      coerce(uts)$ANY1(UTS(FE,x,a))
 
    taylor(an,n,eq,seg) ==
      taylor(eval(an,(n :: FE) = (#1 :: FE)),eq,seg)
 
    laurent(f,eq,seg) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "taylor: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      hasHi seg =>
        n0 := lo seg; n1 := hi seg
        if n1 < n0 then (n0,n1) := (n1,n0)
        uts := series(genFiniteStream(f,n0,n1))$UTS(FE,x,a)
        coerce(laurent(n0,uts)$ULS(FE,x,a))$ANY1(ULS(FE,x,a))
      n0 := lo seg
      uts := series(genStream(f,n0))$UTS(FE,x,a)
      coerce(laurent(n0,uts)$ULS(FE,x,a))$ANY1(ULS(FE,x,a))
 
    laurent(an,n,eq,seg) ==
      laurent(eval(an,(n :: FE) = (#1 :: FE)),eq,seg)
 
    modifyFcn:(RN -> FE,I,I,I,I) -> FE
    modifyFcn(f,n0,nn,q,m) == (zero?((m - n0) rem nn) => f(m/q); 0)
 
    puiseux(f,eq,seg,r) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "puiseux: left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      not positive? r => error "puiseux: last argument must be positive"
      hasHi seg =>
        r0 := lo seg; r1 := hi seg
        if r1 < r0 then (r0,r1) := (r1,r0)
        p0 := numer r0; q0 := denom r0
        p1 := numer r1; q1 := denom r1
        p2 := numer r; q2 := denom r
        q := lcm(lcm(q0,q1),q2)
        n0 := p0 * (q quo q0); n1 := p1 * (q quo q1)
        nn := p2 * (q quo q2)
        ulsUnion := laurent(modifyFcn(f,n0,nn,q,#1),eq,segment(n0,n1))
        uls := retract(ulsUnion)$ANY1(ULS(FE,x,a))
        coerce(puiseux(1/q,uls)$UPXS(FE,x,a))$ANY1(UPXS(FE,x,a))
      p0 := numer(r0 := lo seg); q0 := denom r0
      p2 := numer r; q2 := denom r
      q := lcm(q0,q2)
      n0 := p0 * (q quo q0); nn := p2 * (q quo q2)
      ulsUnion := laurent(modifyFcn(f,n0,nn,q,#1),eq,segment n0)
      uls := retract(ulsUnion)$ANY1(ULS(FE,x,a))
      coerce(puiseux(1/q,uls)$UPXS(FE,x,a))$ANY1(UPXS(FE,x,a))
 
    puiseux(an,n,eq,r0,m) ==
      puiseux(eval(an,(n :: FE) = (#1 :: FE)),eq,r0,m)
 
    series(f:I -> FE,eq:EQ FE) == puiseux(f(numer #1),eq,segment 0,1)
    series(an:FE,n:SY,eq:EQ FE) == puiseux(an,n,eq,segment 0,1)
    series(f:I -> FE,eq:EQ FE,seg:SEG I) ==
      ratSeg : SEG RN := map(#1::RN,seg)$UniversalSegmentFunctions2(I,RN)
      puiseux(f(numer #1),eq,ratSeg,1)
    series(an:FE,n:SY,eq:EQ FE,seg:SEG I) ==
      ratSeg : SEG RN := map(#1::RN,seg)$UniversalSegmentFunctions2(I,RN)
      puiseux(an,n,eq,ratSeg,1)
    series(f:RN -> FE,eq:EQ FE,seg:SEG RN,r:RN) == puiseux(f,eq,seg,r)
    series(an:FE,n:SY,eq:EQ FE,seg:SEG RN,r:RN) == puiseux(an,n,eq,seg,r)

