#
# LazyUnivariatePowerSeries(R)
#
# Create a univariate power series domain over the ring R
#
# Author: Dominik Gruntz 1991
#

LazyUnivariatePowerSeries := proc() local P,R, x:

    R := args[1]:
    x := args[2]:
    if not type(x, name) then ERROR(`2nd arguement must be a name`) fi:

    P := UnivariatePowerSeries(R):

    P[DomainName] := LazyUnivariatePowerSeries:

    env := ['_R' = R, '_P' = P]: 

    P[MakeUPS] := subs(env, proc(p) local coeffproc, Options, i, newproc:
       # p is a procedure!!
       # arg 1 : procedure
       # arg 2 : ldegree   : Startpoint of Laurent Series 
       # arg 3 : order     : Order of Laurent Series (finite or infinite
       # arg 4 : lorder    : first nonzero coefficient
       if type(p,procedure) then
          Options := op(3,p), `LUPS key `.`LUPS/UPSkey`:
          `LUPS/UPSkey` := `LUPS/UPSkey` + 1:
          coeffproc := subsop(3 = Options, p): newproc := true:
       else coeffproc := p: newproc := false: fi:

       coeffproc(Type)    := _P:
       coeffproc(ldegree) := simplify(args[2]):
       coeffproc(order)   := simplify(args[3]):

       if   args[2]=args[3] then # monomial
          coeffproc(lorder) := coeffproc(ldegree):
       elif nargs=4 and not type(args[4], list) then 
          if has(args[4], notDefined) then
             coeffproc(lorder) := notDefined
           else
             coeffproc(lorder) := simplify(args[4])
          fi:
       else coeffproc(lorder) := notDefined 
       fi:

       if type(args[nargs],list) then
          for i to nops(args[nargs]) do
             coeffproc(op(1,args[nargs][i])) := op(2,args[nargs][i]):
          od
       fi:
       if newproc then _ups(op(coeffproc)) else _ups(p) fi: 
    end):

    P[FixedPoint] := subs(env, proc(F) local e,p: # F : procedure
       if not type(F,procedure) then
	  ERROR(`first argument must be of type procedure`) fi:
          e := _P[MakeUPS](p, seq(args[i], i = 2..nargs)):
          e := F(e):
          p := op(1,e):
       op(1,e)(order) := simplify(args[3]):
       if op(1,e)(ldegree) <> args[2] then lprint(`error in FixedPoint`) fi:
       if nargs=4 then op(1,e)(lorder) := args[4] fi:
       e:
    end):

    P[Variable] := x:
    
    P[Coeff] := subs(env, proc(p, i):
       if i < _P[Ldegree](p) then _R[0] else eval(op(1,p)(i)) fi: 
       #op(1,p)(i)
    end):

    P[order]   := subs(env, proc(p) op(1,p)(order) end):

    P[Ldegree] := subs(env, proc(p) op(1,p)(ldegree) end):

    P[Lorder] := subs(env, proc(p) local i,next_proc:
       if op(1,p)(lorder) <> notDefined then RETURN(op(1,p)(lorder)) fi:
       for i from op(1,p)(ldegree) while _R[`=`](_P[Coeff](p,i), _R[0]) do od:
       op(1,p)(lorder) := i:
       op(1,p)(ldegree) := i:
       RETURN(i)
    end):

    P[Monomial] :=     subs(env, proc(n) option remember: 
       if nargs=0 then RETURN(_P[Monomial](1)) fi:
       _P[MakeUPS](proc(k) _R[0] end, n, n, n, [(n)=_R[1]]):
    end):

#   defOperation( Constant, R &-> P, P ):
    P[Constant] := subs(env, proc(r):
       #if not _R[Type](r) then ERROR(`defined as R &-> P`) fi:  
       _P[MakeUPS](proc(k) _R[0] end, 0, 0, 0,[(0)=r]): 
    end):

    P[Shift] := subs(env, proc(p, n) 
       if n=0 then p else
          _P[MakeUPS](subs(['_p'=p, '_n'=n],
             proc(k) option remember: _P[Coeff](_p,k-_n) end),
          _P[Ldegree](p)+n,
          _P[order](p)+n,
          op(1,p)(lorder)+n)
       fi:
    end):

#   defOperation(`R*`, [R,P] &-> P, P)
    P[`R*`] := subs(env, proc(a,u):
       if not ( _R[Type](a) and _P[Type](u) ) then
          ERROR(` defined as [R,P] &-> P, P)`):
       fi:
       _P[MakeUPS](subs(['_a'=a,'_u'=u], 
          proc(k) option remember: _R[`*`](_a,_P[Coeff](_u,k)) end),
          _P[Ldegree](u),
          _P[order](u)
       )
    end):

#   defOperation( Series, List(R,Integer) &-> P, P ):
    P[Series] := subs(env, proc(L) local exponents, coefficients:
       if not type(L,listlist) or op({op(map(nops,L))}) <> 2 then 
          ERROR(`defined as List(R,Integer) &-> P`) fi:
       exponents := map(x -> op(2,x), L): 
       if not type(exponents, list(integer)) then 
ERROR(`defined as List(R,Integer) &-> P: exponents must be integers.`): fi:
       coefficients := map(x -> op(1,x), L):
       for i to nops(coefficients) do
          if not _R[Type](coefficients[i]) then
             ERROR(`defined as List(R,Integer) &-> P: `,
                   `coefficients must belong to the coefficient ring.`): 
          fi:
       od:

       _P[MakeUPS](
           proc(k) _R[0] end, 
           min(op(exponents)), 
           max(op(exponents)),     
           map(x -> (op(2,x))=op(1,x) , L) 
       )
    end):
                

# Implementation of Set() - Functions:
# ------------------------------------

    P[Random] := subs(env, proc():
       _P[MakeUPS](proc(k) option remember: _R[Random]() end,
          - rand(10)(), 
          infinity
       ):
    end):

    P[Input] := subs(env, proc(expr) local i,t,p: # a is a procedure or a list
       if not has(expr, _P[Variable]) then
          _P[Constant](_R[Input](expr))
       elif expr = _P[Variable] then _P[Monomial](1)
       elif type(expr, series) then _P[Input](convert(expr, polynom))
       elif type(expr, `^`) then
          if (op(1,expr) = _P[Variable]) and type(op(2,expr), integer) then
             _P[Monomial](op(2,expr))
          else
             t := _P[Input](op(1,expr)):
             if type(op(2,expr), rational) then
                _P[`^`](t, op(2,expr))
             else
                _P[`^`](t, _P[Input](op(2,expr)))
             fi:
          fi
       elif type(expr, `*`) then
          t := _P[Input](op(1,expr)):
          for i from 2 to nops(expr) do
             t := _P[`*`](t, _P[Input](op(i,expr)))
          od: t
       elif type(expr, polynom(anything, _P[Variable])) then
          p := [coeffs(expr, _P[Variable], t)]:
          p := map(x -> _R[Input](x), p):
          seq([p[i], degree(t[i])], i=1..nops(p)):
          _P[Series](["]):
       elif type(expr, `+`) then
          t := _P[Input](op(1,expr)):
          for i from 2 to nops(expr) do
             t := _P[`+`](t, _P[Input](op(i,expr)))
          od: t
       elif type(expr, function) then
          if op(0,expr) = O then 
             i := _P[Lorder](_P[Input](op(1,expr))):
             _P[MakeUPS](subs(['_i'=i],
                proc(k) 
                  if k<_i then _R[0] 
                  else ERROR(`this term is not available.`) 
                  fi 
                end),
                i, infinity, i
             ):
          elif op(0,expr) = exp  then _P[Exp](_P[Input](op(1,expr)))
          elif op(0,expr) = ln   then _P[Ln ](_P[Input](op(1,expr)))
          elif op(0,expr) = log  then _P[Log](_P[Input](op(1,expr)))
          elif op(0,expr) = sin  then _P[Sin](_P[Input](op(1,expr)))
          elif op(0,expr) = cos  then _P[Cos](_P[Input](op(1,expr)))
          elif op(0,expr) = tan  then _P[Tan](_P[Input](op(1,expr)))
          elif op(0,expr) = sinh then _P[Sinh](_P[Input](op(1,expr)))
          elif op(0,expr) = cosh then _P[Cosh](_P[Input](op(1,expr)))
          elif op(0,expr) = tanh then _P[Tanh](_P[Input](op(1,expr)))
          else FAIL
          fi:
       else FAIL
       fi:
    end):

    P[Output] := subs(env, proc(p) local d, v, t, c, k, ldeg:
       if has(p,FAIL) then RETURN(FAIL) fi:
       if nargs=2 then d := args[2] else d := Order-1 fi:
       v := _P[Variable]:
       ldeg := _P[Ldegree](p): if ldeg = infinity then RETURN(0) fi:
       t := 0: # Maple Variable
       for k from ldeg to ldeg+d do
          c := traperror(_P[Coeff](p,k)):
          if lasterror = `this term is not available.` then
             RETURN(t + O(v^k))
          fi:
          t := t + _R[Output](c)*v^k
       od:
       if   _P[order](p)=infinity or _P[order](p) > d+ldeg then 
          t + O(v^(d+1+ldeg))
       else
          t
       fi:
    end):

    P[Type] := subs(env, proc(p) 
       type(p,function) and type(op(1,p),procedure) and op(1,p)(Type)=_P 
    end):

    P[`=`] := subs(env, proc(u,v) local i,n,lu,lv:
       if _P[Ldegree](u)=infinity and _P[Ldegree](v)=infinity then 
          RETURN(true) # both are 0
       fi:
       lu := _P[Lorder](u): lv := _P[Lorder](v):
       if lu <> lv then 
          RETURN(false)
       else
          if nargs=2 then n := Order else n := args[3] fi:
          for i from lu to lu+n do
             if _R[`<>`](_P[Coeff](u,i),_P[Coeff](v,i)) then RETURN(false) fi:
          od:
          lprint(`equal up to `.n.` Terms`):
          RETURN(true)
       fi:
    end):

# Implementation of SemiGroup() / Monoid() / Group() - Functions:
# ---------------------------------------------------------------

    P[`+`] := subs(env, proc(p1,p2):
       if   _P[Ldegree](p1)=infinity then RETURN(p2)
       elif _P[Ldegree](p2)=infinity then RETURN(p1)
       fi:
       _P[MakeUPS](subs(['_p1'=p1,'_p2'=p2],
          proc(k) option remember:
            _R[`+`](_P[Coeff](_p1,k),_P[Coeff](_p2,k)) 
          end),
          min(_P[Ldegree](p1), _P[Ldegree](p2)),
          max(_P[order](p1),   _P[order](p2))
       ): 
    end):

    P[`-`] := subs(env, proc(p1) local p2:
       if nargs=2 then p2 := args[2]:
          if   _P[Ldegree](p2)=infinity then RETURN(p1)
          elif _P[Ldegree](p1)=infinity then RETURN(_P[`-`](p2))
          fi:
          _P[MakeUPS](subs(['_p1'=p1,'_p2'=p2],
             proc(k) option remember:
                _R[`-`](_P[Coeff](_p1,k),_P[Coeff](_p2,k)) 
             end),
             min(_P[Ldegree](p1), _P[Ldegree](p2)),
             max(_P[order](p1),   _P[order](p2))
          )
       else            
          _P[MakeUPS](subs(['_p1'=p1],
             proc(k) option remember:
                _R[`-`](_P[Coeff](_p1,k)) 
            end),
            _P[Ldegree](p1),
            _P[order](p1)
          )
       fi:
    end):

P[0] := subs(env, _P[MakeUPS](proc(x) _R[0] end, infinity,infinity,infinity)):

# Implementation of Ring() - Functions:
# -------------------------------------

    P[1] := P[Monomial](0):

# defined in Ring()  as [Integer,R] &-> R, [R,R] &-> R 
    
    P[`*`] := subs(env, proc(u,v) local lu,lv,lw:
       if not ( (type(u,integer) or _P[Type](u)) and _P[Type](v) ) then
          ERROR(`defined as [Integer,R] &-> R or [R,R] &-> R`):
       fi:
       if type(u,integer) then
          _P[MakeUPS](subs(['_a'=u,'_v'=v], 
             proc(k) option remember: _R[`*`](_a,_P[Coeff](_v,k)) end),
             _P[Ldegree](v),
             _P[order](v)
          )
       elif _P[Ldegree](u)=infinity or _P[Ldegree](v)=infinity then _P[0]
       elif _P[Ldegree](u)=_P[order](u) then # u is a monomial
          _P[Shift](_P[`R*`](_P[Coeff](u,_P[Ldegree](u)),v), _P[Ldegree](u))
       elif _P[Ldegree](v)=_P[order](v) then # v is a monomial
          _P[Shift](_P[`R*`](_P[Coeff](v,_P[Ldegree](v)),u), _P[Ldegree](v))
       else
          lu := _P[Ldegree](u):
          lv := _P[Ldegree](v):
          lw := simplify(lu + lv):
          _P[MakeUPS](subs(['_u'=u, '_v'=v, '_lu'=lu, '_lv'=lv, '_lw'=lw],
             proc(k) local i,t: option remember:
                t := _R[0]:
                for i from 0 to k-_lw do
                   t := _R[`+`](t, _R[`*`](_P[Coeff](_u,i + _lu), 
                           _P[Coeff](_v,(k-_lw)-i + _lv)))
                od:
                t
             end),
             lw,
             _P[order](u)+_P[order](v)
          )
       fi
    end):

    P[Inv] := subs(env, proc(v) local w, lw, lv, inv_v0:
       if _P[Lorder](v)=infinity then RETURN(FAIL) fi: # Division by zero
       lv := op(1,v)(lorder): lw := -lv:
       inv_v0 := _R[Inv](_P[Coeff](v,lv)): if inv_v0=FAIL then RETURN(FAIL) fi:
       if lv=_P[order](v) then # inverse of a monomial
          _P[MakeUPS](proc(k) _R[0] end, lw, lw, lw, [(lw)=inv_v0])
       else
#            w := _P[Shift](v, -lv):
#            _P[FixedPoint]('e=_P[`-`](_P[Shift](_P[Constant](inv_v0),0), 
#                        _P[Integrate](_P[`*`](_P[`*`](e,e),
#                                      _P[Diff](w))))',0,infinity):
#            _P[Shift](", -lv)

w := _P[MakeUPS](subs(['_v'=v, '_w'=w, '_inv_v0'=inv_v0, '_lv'=lv, '_lw'=lw],
             proc(k) local i,t: option remember:
                t := _R[0]:
                for i from 0 to (k-_lw)-1 do
                   t := _R[`+`](t, _R[`*`](_P[Coeff](_w, i        + _lw),
                                   _P[Coeff](_v,(k-_lw)-i + _lv)))
                od:
                if k=_lw then
                   _R[`*`](_R[`-`](_R[Coerce](1),t) , _inv_v0)
                else
                   _R[`*`](_R[`-`](t) , _inv_v0)
                fi:
             end),
             lw,
             infinity
          ):
       fi:
    end):

#  defOperation( `^`, {[P,Integer]     &-> P,
#                      [P,P]           &-> P,
#                      [P,Rational]    &-> P},         P ):
    P[`^`] := subs(env, proc(u, alpha) local m, w0, u1, u2, w, o:
       if _P[Type](alpha) then
          _P[Exp](_P[`*`](alpha, _P[Ln](u)))
       else
          if _P[Ldegree](u)=infinity then # u = 0
             if alpha=0 then RETURN(FAIL) else RETURN(_P[0]) fi:
                     elif alpha=0 then RETURN(_P[1])
          elif alpha=1 then RETURN(u)
          else
             m := _P[Lorder](u):
             if not type(alpha*m,integer) then RETURN(FAIL) fi:
             #w0 := _R[`^`](_P[Coeff](u,m), alpha):
             if _P[order](u)=m then # monomial ^ alpha
                _P[`R*`](_R[`^`](_P[Coeff](u,m), alpha), _P[Monomial](alpha*m))
             else
                if type(alpha,integer) and (alpha>0) then 
                   o := (_P[order](u)-m)*alpha
                else 
                   o := infinity 
                fi:
             u1 := _P[Shift](u,-m):
             if type(alpha,integer) then
       u2 := _P[FixedPoint](subs(['_alpha'=alpha,'_m'=m,'_w0'=w0, '_u1'=u1],
                         w -> _P[`+`](
                         _P[Constant](_R[`^`](_P[Coeff](_u1,0), _alpha)),
                             _P[`*`](_alpha, 
                                _P[Integrate](
                                   _P[`*`](w,
                                   _P[`/`](_P[Diff](_u1),_u1)))))),
                                          0,o):
             else
                w0 := _R[`^`](_P[Coeff](u1,0), alpha):
                   u2 := _P[FixedPoint](subs(['_alpha'=alpha,'_u1'=u1,'_w0'=w0],
                            w -> _P[`+`](
                        _P[Constant](_R[`^`](_P[Coeff](_u1,0), _alpha)),
                        _P[`/`](
                               _P[`*`](numer(_alpha), 
                                  _P[Integrate](
                              _P[`*`](w,
                              _P[`/`](_P[Diff](_u1),_u1)))),denom(_alpha)))),
                                         0,o):
             fi:
            _P[Shift](u2,alpha*m)
             fi:
          fi
       fi
    end):

    P[Diff] := subs(env, proc(u) local ldeg:
       ldeg := _P[Ldegree](u)-1: if ldeg = -1 then ldeg := 0 fi:
       _P[MakeUPS](subs(['_u'=u],
          proc(k) option remember:
             _R[`*`](k+1,_P[Coeff](_u,k+1))
          end),
          ldeg,
          _P[order](u)-1
       )
    end):

    P[Integrate] := subs(env, proc(u) local ldeg:
       if _P[Ldegree](u)<0 then
          if _R[`<>`](_P[Coeff](u,-1),_R[0]) then RETURN(FAIL) fi
       fi:
       ldeg := _P[Ldegree](u)+1: # if ldeg = 1 then ldeg := 0 fi:
       _P[MakeUPS](subs(['_u'=u],
          proc(k) option remember:
             if k=0 then _R[0]
             else _R[`/`](_P[Coeff](_u,k-1),k)
             fi:
          end),
          ldeg,
          _P[order](u)+1
       )
    end):


    #    diff(exp(p(x)), x) = exp(p(x)) * diff(p(x), x)
    # => exp(p(x))          = int(exp(p(x)) * diff(p(x), x), x) + exp(p(0))
    P[Exp] := subs(env, proc(s):
       if   _P[Ldegree](s)=infinity then RETURN(_P[1])
       elif _P[Lorder](s) < 1       then RETURN(FAIL)
       else
          _P[FixedPoint](subs(['_s'=s],
             e -> _P[`+`](
               _P[Integrate](
                  _P[`*`](
                     _P[Diff](_s), 
                     e 
                  )
               ),     
               _P[1]
                )), 
           0, infinity)
       fi:
    end):

    P[Ln]  := subs(env, proc(s) 
       if   _P[Lorder](s)<>0               then RETURN(FAIL)
       elif _R[`<>`](_P[Coeff](s,0),_R[1]) then RETURN(FAIL)
       else _P[Integrate](_P[`/`](_P[Diff](s),s))
       fi
    end):

    P[Tan] := subs(env, proc(s) local e:
       if   _P[Ldegree](s)=infinity then RETURN(_P[0])
       elif _P[Lorder](s) < 1       then RETURN(FAIL)
       else
          _P[FixedPoint](subs(['_s'=s],
             e -> _P[`+`](
                _P[Integrate](
               _P[`*`](
                  _P[`+`](
                     _P[1],
                     _P[`*`](e,e)
                  ),
                  _P[Diff](_s)
                    )  ),         
                _P[0]
                    )),
             _P[Ldegree](s), infinity)
       fi
    end):

    P[Sin] := subs(env, proc(s) local e, ds: 
       if   _P[Ldegree](s)=infinity then RETURN(_P[0])
       elif _P[Lorder](s) < 1       then RETURN(FAIL)
       else
          ds := _P[Diff](s):
          _P[FixedPoint](subs(['_ds'=ds],
             e -> _P[Integrate](
                _P[`*`](
                                   _P[`-`](
                          _P[1],
                      _P[Integrate](_P[`*`](e,_ds))
                       ), 
                   _ds
                             )  )),
             _P[Ldegree](s), infinity)
       fi:
    end):

    P[Cos] := subs(env, proc(s) local e, ds: 
       if   _P[Ldegree](s)=infinity then RETURN(_P[1])
       elif _P[Lorder](s) < 1       then RETURN(FAIL)
       else
          ds := _P[Diff](s):
          _P[FixedPoint](subs(['_ds'=ds],
             e -> _P[`-`](
                _P[1], 
                _P[Integrate](
               _P[`*`](
                  _P[Integrate](_P[`*`](e,_ds)), 
                  _ds
                             )  )  )),
             0, infinity)
       fi 
    end):

    op(P):
end:

`print/_ups` := proc(p) local i,s:
   p(Type)[Output](_ups(p)):
end:

# the fillowing procedure is useful for debugging purposes, but not necessary
showups := proc(ups, arg): 
    if nargs=1 then print(op(1,ups)): print(op(4,op(1,ups)))
    elif arg=`proc` or arg=1 then lprint(op(1,ups))
    elif arg=table or arg=2  then print(op(4,op(1,ups)))
    else print(op(1,ups)): print(op(4,op(1,ups)))
    fi:
end:

`LUPS/UPSkey` := 0:

save LazyUnivariatePowerSeries,
    `print/_ups`,
    `LUPS/UPSkey`,
#   showups,
    
    `LUPS.m`:

quit

