--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.
 
import Field
import NonNegativeInteger
)abbrev package WEIER WeierstrassPreparation
++ Author:William H. Burge
++ Date Created:Sept 1988
++ Date Last Updated:Feb 15 1992
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description: This package implements the Weierstrass preparation
++ theorem f or multivariate power series.
++ weierstrass(v,p) where v is a variable, and p is a
++ TaylorSeries(R) in which the terms
++ of lowest degree s must include c*v**s where c is a constant,s>0,
++ is a list of TaylorSeries coefficients A[i] of the
++ equivalent polynomial
++ A = A[0] + A[1]*v + A[2]*v**2 + ... + A[s-1]*v**(s-1) + v**s
++ such that p=A*B , B being a TaylorSeries of minimum degree 0
WeierstrassPreparation(R): Defn == Impl where
    R : Field
    VarSet==>Symbol
    SMP ==> Polynomial R
    PS  ==> InnerTaylorSeries SMP
    NNI ==> NonNegativeInteger
    ST  ==> Stream
    StS ==> Stream SMP
    STPS==>StreamTaylorSeriesOperations
    STTAYLOR==>StreamTaylorSeriesOperations
    SUP==> SparseUnivariatePolynomial(SMP)
    ST2==>StreamFunctions2
    SMPS==>  TaylorSeries(R)
    L==>List
    null ==> empty?
    likeUniv ==> univariate
    coef ==> coefficient$SUP
    nil ==> empty
 
 
    Defn ==>  with
 
        crest:(NNI->( StS-> StS))
          ++\spad{crest n} is used internally.
        cfirst:(NNI->( StS-> StS))
          ++\spad{cfirst n} is used internally.
        sts2stst:(VarSet,StS)->ST StS
          ++\spad{sts2stst(v,s)} is used internally.
        clikeUniv:VarSet->(SMP->SUP)
          ++\spad{clikeUniv(v)} is used internally.
        weierstrass:(VarSet,SMPS)->L SMPS
          ++\spad{weierstrass(v,ts)} where v is a variable and ts is
          ++ a TaylorSeries, impements the Weierstrass Preparation
          ++ Theorem. The result is a list of TaylorSeries that
          ++ are the coefficients of the equivalent series.
        qqq:(NNI,SMPS,ST SMPS)->((ST SMPS)->ST SMPS)
          ++\spad{qqq(n,s,st)} is used internally.
 
    Impl ==>  add
        import TaylorSeries(R)
        import StreamTaylorSeriesOperations SMP
        import StreamTaylorSeriesOperations SMPS
 
 
        map1==>map$(ST2(SMP,SUP))
        map2==>map$(ST2(StS,SMP))
        map3==>map$(ST2(StS,StS))
        transback:ST SMPS->L SMPS
        transback smps==
            if null smps
            then nil()$(L SMPS)
            else
              if null first (smps:(ST StS))
              then nil()$(L SMPS)
              else
                cons(map2(first,smps:ST StS):SMPS,
                   transback(map3(rest,smps:ST StS):(ST SMPS)))$(L SMPS)
 
 
        clikeUniv(var)==likeUniv(#1,var)
        mind:(NNI,StS)->NNI
        mind(n, sts)==
           if null sts
           then error "no mindegree"
           else if first sts=0
                then mind(n+1,rest sts)
                else n
        mindegree (sts:StS):NNI== mind(0,sts)
 
 
        streamlikeUniv:(SUP,NNI)->StS
        streamlikeUniv(p:SUP,n:NNI): StS ==
          if n=0
          then cons(coef (p,0),nil()$StS)
          else cons(coef (p,n),streamlikeUniv(p,(n-1):NNI))
 
        transpose:ST StS->ST StS
        transpose(s:ST StS)==delay(
           if null s
           then nil()$(ST StS)
           else cons(map2(first,s),transpose(map3(rest,rst s))))
 
        zp==>map$StreamFunctions3(SUP,NNI,StS)
 
        sts2stst(var, sts)==
           zp(streamlikeUniv(#1,#2),
             map1(clikeUniv var, sts),(integers 0):(ST NNI))
 
        tp:(VarSet,StS)->ST StS
        tp(v,sts)==transpose sts2stst(v,sts)
        map4==>map$(ST2 (StS,StS))
        maptake:(NNI,ST StS)->ST SMPS
        maptake(n,p)== map4(cfirst n,p) pretend ST SMPS
        mapdrop:(NNI,ST StS)->ST SMPS
        mapdrop(n,p)== map4(crest n,p) pretend ST SMPS
        YSS==>Y$ParadoxicalCombinatorsForStreams(SMPS)
        weier:(VarSet,StS)->ST SMPS
        weier(v,sts)==
             a:=mindegree sts
             if a=0
             then error "has constant term"
             else
               p:=tp(v,sts) pretend (ST SMPS)
               b:StS:=rest(((first p pretend StS)),a::NNI)
               c:=retractIfCan first b
               c case "failed"=>_
 error "the coefficient of the lowest degree of the variable should _
 be a constant"
               e:=recip b
               f:= if e case "failed"
                   then error "no reciprocal"
                   else e::StS
               q:=(YSS qqq(a,f:SMPS,rest p))
               maptake(a,(p*q) pretend ST StS)
 
        cfirst n== first(#1,n)$StS
        crest n== rest(#1,n)$StS
        qq:(NNI,SMPS,ST SMPS,ST SMPS)->ST SMPS
        qq(a,e,p,c)==
            cons(e,(-e)*mapdrop(a,(p*c)pretend(ST StS)))
        qqq(a,e,p)==  qq(a,e,p,#1)
        wei:(VarSet,SMPS)->ST SMPS
        wei(v:VarSet,s:SMPS)==weier(v,s:StS)
        weierstrass(v,smps)== transback wei (v,smps)

