--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 STINPROD StreamInfiniteProduct
++ Author: Clifton J. Williamson
++ Date Created: 23 February 1990
++ Date Last Updated: 23 February 1990
++ Basic Operations: infiniteProduct, evenInfiniteProduct, oddInfiniteProduct,
++   generalInfiniteProduct
++ Related Domains: UnivariateTaylorSeriesCategory
++ Also See:
++ AMS Classifications:
++ Keywords: Taylor series, infinite product
++ Examples:
++ References:
++ Description: 
++   This package computes infinite products of Taylor series over an
++   integral domain of characteristic 0.  Here Taylor series are
++   represented by streams of Taylor coefficients.
StreamInfiniteProduct(Coef): Exports == Implementation where
  Coef: Join(IntegralDomain,CharacteristicZero)
  I  ==> Integer
  QF ==> Fraction
  ST ==> Stream
 
  Exports ==> with
 
    infiniteProduct: ST Coef -> ST Coef
      ++ infiniteProduct(f(x)) computes \spad{product(n=1,2,3...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    evenInfiniteProduct: ST Coef -> ST Coef
      ++ evenInfiniteProduct(f(x)) computes \spad{product(n=2,4,6...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    oddInfiniteProduct: ST Coef -> ST Coef
      ++ oddInfiniteProduct(f(x)) computes \spad{product(n=1,3,5...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    generalInfiniteProduct: (ST Coef,I,I) -> ST Coef
      ++ generalInfiniteProduct(f(x),a,d) computes
      ++ \spad{product(n=a,a+d,a+2*d,...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
 
  Implementation ==> add
 
    if Coef has Field then
 
      import StreamTaylorSeriesOperations(Coef)
      import StreamTranscendentalFunctions(Coef)
 
      infiniteProduct st             == exp lambert log st
      evenInfiniteProduct st         == exp evenlambert log st
      oddInfiniteProduct st          == exp oddlambert log st
      generalInfiniteProduct(st,a,d) == exp generalLambert(log st,a,d)
 
    else
 
      import StreamTaylorSeriesOperations(QF Coef)
      import StreamTranscendentalFunctions(QF Coef)
 
      applyOverQF:(ST QF Coef -> ST QF Coef,ST Coef) -> ST Coef
      applyOverQF(f,st) ==
        stQF := map(#1 :: QF(Coef),st)$StreamFunctions2(Coef,QF Coef)
        map(retract(#1)@Coef,f stQF)$StreamFunctions2(QF Coef,Coef)
 
      infiniteProduct st     == applyOverQF(exp lambert log #1,st)
      evenInfiniteProduct st == applyOverQF(exp evenlambert log #1,st)
      oddInfiniteProduct st  == applyOverQF(exp oddlambert log #1,st)
      generalInfiniteProduct(st,a,d) ==
        applyOverQF(exp generalLambert(log #1,a,d),st)

)abbrev package INFPROD0 InfiniteProductCharacteristicZero
++ Author: Clifton J. Williamson
++ Date Created: 22 February 1990
++ Date Last Updated: 23 February 1990
++ Basic Operations: infiniteProduct, evenInfiniteProduct, oddInfiniteProduct,
++   generalInfiniteProduct
++ Related Domains: UnivariateTaylorSeriesCategory
++ Also See:
++ AMS Classifications:
++ Keywords: Taylor series, infinite product
++ Examples:
++ References:
++ Description: 
++   This package computes infinite products of univariate Taylor series
++   over an integral domain of characteristic 0.
InfiniteProductCharacteristicZero(Coef,UTS):_
 Exports == Implementation where
  Coef : Join(IntegralDomain,CharacteristicZero)
  UTS  : UnivariateTaylorSeriesCategory Coef
  I  ==> Integer
 
  Exports ==> with
 
    infiniteProduct: UTS -> UTS
      ++ infiniteProduct(f(x)) computes \spad{product(n=1,2,3...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    evenInfiniteProduct: UTS -> UTS
      ++ evenInfiniteProduct(f(x)) computes \spad{product(n=2,4,6...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    oddInfiniteProduct: UTS -> UTS
      ++ oddInfiniteProduct(f(x)) computes \spad{product(n=1,3,5...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    generalInfiniteProduct: (UTS,I,I) -> UTS
      ++ generalInfiniteProduct(f(x),a,d) computes
      ++ \spad{product(n=a,a+d,a+2*d,...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
 
  Implementation ==> add
 
    import StreamInfiniteProduct Coef
 
    infiniteProduct x     == series infiniteProduct coefficients x
    evenInfiniteProduct x == series evenInfiniteProduct coefficients x
    oddInfiniteProduct x  == series oddInfiniteProduct coefficients x
 
    generalInfiniteProduct(x,a,d) ==
      series generalInfiniteProduct(coefficients x,a,d)

)abbrev package INPRODPF InfiniteProductPrimeField
++ Author: Clifton J. Williamson
++ Date Created: 22 February 1990
++ Date Last Updated: 23 February 1990
++ Basic Operations: infiniteProduct, evenInfiniteProduct, oddInfiniteProduct,
++   generalInfiniteProduct
++ Related Domains: UnivariateTaylorSeriesCategory
++ Also See:
++ AMS Classifications:
++ Keywords: Taylor series, infinite product
++ Examples:
++ References:
++ Description: 
++    This package computes infinite products of univariate Taylor series
++    over a field of prime order.
InfiniteProductPrimeField(Coef,UTS): Exports == Implementation where
  Coef : Join(Field,Finite,ConvertibleTo Integer)
  UTS  : UnivariateTaylorSeriesCategory Coef
  I  ==> Integer
  ST ==> Stream
 
  Exports ==> with
 
    infiniteProduct: UTS -> UTS
      ++ infiniteProduct(f(x)) computes \spad{product(n=1,2,3...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    evenInfiniteProduct: UTS -> UTS
      ++ evenInfiniteProduct(f(x)) computes \spad{product(n=2,4,6...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    oddInfiniteProduct: UTS -> UTS
      ++ oddInfiniteProduct(f(x)) computes \spad{product(n=1,3,5...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    generalInfiniteProduct: (UTS,I,I) -> UTS
      ++ generalInfiniteProduct(f(x),a,d) computes
      ++ \spad{product(n=a,a+d,a+2*d,...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
 
  Implementation ==> add
 
    import StreamInfiniteProduct Integer
 
    applyOverZ:(ST I -> ST I,ST Coef) -> ST Coef
    applyOverZ(f,st) ==
      stZ := map(convert(#1)@Integer,st)$StreamFunctions2(Coef,I)
      map(#1 :: Coef,f stZ)$StreamFunctions2(I,Coef)
 
    infiniteProduct x ==
      series applyOverZ(infiniteProduct,coefficients x)
    evenInfiniteProduct x ==
      series applyOverZ(evenInfiniteProduct,coefficients x)
    oddInfiniteProduct x ==
      series applyOverZ(oddInfiniteProduct,coefficients x)
    generalInfiniteProduct(x,a,d) ==
      series applyOverZ(generalInfiniteProduct(#1,a,d),coefficients x)

)abbrev package INPRODFF InfiniteProductFiniteField
++ Author: Clifton J. Williamson
++ Date Created: 22 February 1990
++ Date Last Updated: 23 February 1990
++ Basic Operations: infiniteProduct, evenInfiniteProduct, oddInfiniteProduct,
++   generalInfiniteProduct
++ Related Domains: UnivariateTaylorSeriesCategory
++ Also See:
++ AMS Classifications:
++ Keywords: Taylor series, infinite product
++ Examples:
++ References:
++ Description: 
++   This package computes infinite products of univariate Taylor series
++   over an arbitrary finite field.
InfiniteProductFiniteField(K,UP,Coef,UTS):_
 Exports == Implementation where
  K    :  Join(Field,Finite,ConvertibleTo Integer)
  UP   :  UnivariatePolynomialCategory K
  Coef :  MonogenicAlgebra(K,UP)
  UTS  :  UnivariateTaylorSeriesCategory Coef
  I   ==> Integer
  RN  ==> Fraction Integer
  SAE ==> SimpleAlgebraicExtension
  ST  ==> Stream
  STF ==> StreamTranscendentalFunctions
  STT ==> StreamTaylorSeriesOperations
  ST2 ==> StreamFunctions2
  SUP ==> SparseUnivariatePolynomial
 
  Exports ==> with
 
    infiniteProduct: UTS -> UTS
      ++ infiniteProduct(f(x)) computes \spad{product(n=1,2,3...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    evenInfiniteProduct: UTS -> UTS
      ++ evenInfiniteProduct(f(x)) computes \spad{product(n=2,4,6...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    oddInfiniteProduct: UTS -> UTS
      ++ oddInfiniteProduct(f(x)) computes \spad{product(n=1,3,5...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
    generalInfiniteProduct: (UTS,I,I) -> UTS
      ++ generalInfiniteProduct(f(x),a,d) computes
      ++ \spad{product(n=a,a+d,a+2*d,...,f(x**n))}.
      ++ The series \spad{f(x)} should have constant coefficient 1.
 
  Implementation ==> add
 
    liftPoly: UP -> SUP RN
    liftPoly poly ==
      -- lift coefficients of 'poly' to integers
      ans : SUP RN := 0
      while not zero? poly repeat
        coef := convert(leadingCoefficient poly)@I :: RN
        ans := ans + monomial(coef,degree poly)
        poly := reductum poly
      ans
 
    reducePoly: SUP RN -> UP
    reducePoly poly ==
      -- reduce coefficients of 'poly' to elements of K
      ans : UP := 0
      while not zero? poly repeat
        coef := numer(leadingCoefficient(poly)) :: K
        ans := ans + monomial(coef,degree poly)
        poly := reductum poly
      ans
 
    POLY := liftPoly definingPolynomial()$Coef
    ALG  := SAE(RN,SUP RN,POLY)
 
    infiniteProduct x ==
      stUP := map(lift,coefficients x)$ST2(Coef,UP)
      stSUP := map(liftPoly,stUP)$ST2(UP,SUP RN)
      stALG := map(reduce,stSUP)$ST2(SUP RN,ALG)
      stALG := exp(lambert(log(stALG)$STF(ALG))$STT(ALG))$STF(ALG)
      stSUP := map(lift,stALG)$ST2(ALG,SUP RN)
      stUP := map(reducePoly,stSUP)$ST2(SUP RN,UP)
      series map(reduce,stUP)$ST2(UP,Coef)
 
    evenInfiniteProduct x ==
      stUP := map(lift,coefficients x)$ST2(Coef,UP)
      stSUP := map(liftPoly,stUP)$ST2(UP,SUP RN)
      stALG := map(reduce,stSUP)$ST2(SUP RN,ALG)
      stALG := exp(evenlambert(log(stALG)$STF(ALG))$STT(ALG))$STF(ALG)
      stSUP := map(lift,stALG)$ST2(ALG,SUP RN)
      stUP := map(reducePoly,stSUP)$ST2(SUP RN,UP)
      series map(reduce,stUP)$ST2(UP,Coef)
 
    oddInfiniteProduct x ==
      stUP := map(lift,coefficients x)$ST2(Coef,UP)
      stSUP := map(liftPoly,stUP)$ST2(UP,SUP RN)
      stALG := map(reduce,stSUP)$ST2(SUP RN,ALG)
      stALG := exp(oddlambert(log(stALG)$STF(ALG))$STT(ALG))$STF(ALG)
      stSUP := map(lift,stALG)$ST2(ALG,SUP RN)
      stUP := map(reducePoly,stSUP)$ST2(SUP RN,UP)
      series map(reduce,stUP)$ST2(UP,Coef)
 
    generalInfiniteProduct(x,a,d) ==
      stUP := map(lift,coefficients x)$ST2(Coef,UP)
      stSUP := map(liftPoly,stUP)$ST2(UP,SUP RN)
      stALG := map(reduce,stSUP)$ST2(SUP RN,ALG)
      stALG := generalLambert(log(stALG)$STF(ALG),a,d)$STT(ALG)
      stALG := exp(stALG)$STF(ALG)
      stSUP := map(lift,stALG)$ST2(ALG,SUP RN)
      stUP := map(reducePoly,stSUP)$ST2(SUP RN,UP)
      series map(reduce,stUP)$ST2(UP,Coef)

