DiffeqToRecurrence.ml

(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)
(* DDMF is distributed under CeCILL-B license, *)
(* however this file is distributed under the LGPL version 2.1. *)

INCLUDE "preamble.ml"

let title (_, _, _, _, notation) =
  <:text<Finding a Recurrence for the Power Series Coefficients
         of <:imath< $(str: notation) >>>>

(* The following exception is used from goodinitvalues_rec to diffeqtorec. *)
exception No_solution

(* The following exception is used internally in goodinitvalues_rec only. *)
(* FC: This is how I understood the former code when porting it to ocaml. *)
(* But it is really strange. *)
exception Return_without_display of any maple

(*
# `goodinitvalues/rec`
# Input:  a recurrence in the format returned by gfun:-formatrec
#    the unknown sequence and its variable
#    some initial conditions
#    a boolean flag cleanup
#    (optional) an integer p
#    (optional) maxsing = nonnegint
# Output: a set of equalities u(k)=v_k, from which all the other values
#    (smaller than maxsing) can be deduced by solving the recurrence for its
#     maximal index.
#    This set is continued up to the pth term if p is given.
#    New variables _C[i] are introduced when a value is arbitrary.
#    When the flag is true, the unnecessary entries u(k)=_C[i] are removed
#    The result is an ERROR when no initial condition can be found.
*)

let goodinitvalues_rec recurrence u n ini flag argsix =

   (* :=proc (recurrence, u, n, ini, flag, arg6)
   option `Copyright (c) 1992-2009 by Algorithms Project, INRIA France. \
All rights reserved.`;
   local n0, order, i, inds, minind, maxind, sys, r, sol, b, a, j, k,
   rej, maxsingularity, gb, termorder, dorej, inds2; *)

  <:unit<
    # TODO: remove the following, when these are local vars
    unassign('n0', 'i', 'inds', 'minind', 'maxind', 'sys', 'r', 'sol', 'b',
             'a', 'j', 'k',
     'rej', 'maxsingularity', 'gb', 'termorder', 'dorej', 'inds2');

    recOrder:=nops($(recurrence))-2;  # SG: replaced order by recOrder
    maxind:=recOrder-1;
    if type($(ini),'set') then
      inds:=map(op,indets($(ini),$(u)('integer')))
    else
      inds:={}; maxind:=max(maxind,$(ini))
    end if;
    maxsingularity:=NULL;

    if $(argsix)<>NULL then    # SG: replaced "nargs=6"
      if type($(argsix),integer) then
        maxind:=max(maxind,$(argsix))
      elif type($(argsix),identical('maxsing')=nonnegint) then
        maxsingularity:=op(2,$(argsix))
      else error "invalid argument",$(argsix)
      fi
    end if;

    # The initial conditions are u(minind)..$(u)(maxind).

    if maxsingularity<>NULL then
      n0 := gfun:-firstnonzero(
        subs($(n) = $(n)-recOrder,
             $(recurrence)[nops($(recurrence))]),
        $(n),
        maxsingularity);
    else
      n0 := gfun:-firstnonzero(
        subs($(n) = $(n)-recOrder,
             $(recurrence)[nops($(recurrence))]),
        $(n));
    end if;

    # u(n0-1) cannot be deduced from the previous ones.
    maxind:=max(maxind,op(inds),n0-1);
    minind:=min(op(inds),0);
    r:=gfun:-makerec($(recurrence),$(u),$(n));
  >> ;

  (* SG: breakFlag instead of "break" statement in the following loop *)
  <:unit< breakFlag := false: >> ;

  (* Some cases will return a computed value but no do any display. *)
  (* An exception Return_without_display is used in this case. *)
  (* The 'catch' is at the very bottom of the function. *)
  try

    let sol_found = ref false and j = ref 1 in
    while not !sol_found && !j < 3 do
(*
    <* j := $(int:j): *> (* added by SG; j is used below *)
*)

      if !j = 1 then
        <:unit<
          sys := { op($(ini)),
                   seq(subs($(n)=i,r), i=minind..maxind-recOrder) } :
        >>
      else
        <:unit<
          sys := { op($(ini)),
                   seq(subs($(n)=i-recOrder,r),
                       i={$minind+recOrder..maxind} minus inds) } :
        >> ;

      if <:bool< sys = {} >> then
        raise
          (Return_without_display
             (if not flag then
               << { seq($(u)(i) = _C[i], i=minind..maxind) } >>
             else
               <<
                 { seq($(u)(i) = _C[i],
                       i={$minind..maxind} minus {$0..recOrder-1}) }
               >>)) ;

      <:unit<
        a := gfun:-systomatrix(sys,[seq($(u)(i),i=minind..maxind)],'b');
      >> ;
      sol_found :=
        (* FIXME: We need a DynaMoW:-Raise badly enough! *)
        (* But will extraction of Maple code be easy? *)
        <:bool<
          try
            sol := LinearAlgebra:-LinearSolve(a,Vector['column'](b))
          catch "inconsistent system" :
            sol := NULL
          end try ;
          evalb(sol <> NULL)
        >> ;

      if not !sol_found && !j = 2 then
        <:unit< error "no valid initial conditions": >> ;

      incr j
    done ;

    (* SG: added the following if *)
    if not !sol_found then
      raise No_solution ;

    (* added by SG; TODO should become local var *)
    <:unit< sysOld :=sys : >> ;

    (* TODO: replaced k by kk in the following; should become local var *)
    <:unit<
      sol:=convert(sol,list);
  #    inds:=indets(sol,_t[anything]);
  #    inds:=indets(sol,'typeindex'( 'anything', 'suffixed'('_t', 'integer') ));
  #  SG: added the Or; before, e.g., _t[1] would not be replaced
      inds := indets(sol,
                     'typeindex'('anything',
                                 Or(identical('_t'),
                                    'suffixed'('_t', 'integer')) ));

      # replace the _t[anything] by _C[anything] depending on flag
      dorej:=$(bool:flag) and evalb(j=1); # SG: added evalb

      j:=max(op(map(op,indets([$(recurrence),$(ini)],
                              _C['anything']))));

      if j=-infinity then j:=-1 end if;
      for i in inds do
        if member(i,sol,'kk') and (not dorej or nops(select(has,sol,i))>1
                                     or kk<minind-1 or kk+minind>recOrder) then
          j:=j+1;
          sol:=subs(i=_C[minind+j],sol);
          rej[i]:=NULL
        else rej[i]:=kk
        end if
      end do;

      sys:={seq($(u)(i+minind-1)=sol[i],i={$1..nops(sol)} minus
                  {seq(rej[i],i=inds)})};
    >> ;

  (* The end of this function is structured as: *)
  (*   begin looong_display end, begin looong_symbolic end *)
  (* catch Return_without_display symbolic_system -> *)
  (*  <:par<>>, symbolic_system *)
  (* The unusual begin/end are used to mitigate the overload of nested *)
  (* ifs and structures. *)

  (* Start of long text computation, to be returned. *)
  begin
    if <:bool< n0 > recOrder >> then
      <:par<The first index from which on the leading coefficient is non-zero
        is <:isymb< $(n) = n0-recOrder >>.
        The value <:isymb< $(u)(n0-1) >>
        cannot be deduced from the previous ones by the recurrence. Hence >> @:@
      (let bound = <:int< n0-recOrder-1 >> in
      if <:bool< n0-recOrder-1 = 0 >> then
         <:par<the instance <:imath< $(symb:n) = 0 >> of the recurrence becomes
           an initial condition.>>
      else
        <:par<the instances
          <:imath< $(symb:n) = 0, \dots,~$(symb:n) = $(int:bound) >>
          of the recurrence become inital conditions.>>) @:@
      (if <:bool< maxind > minind >> then
        <:par<The initial sequence elements
          <:imath< $(symb:u)(minind), \dots,~$(symb:u)(maxind) >>
          thus satisfy the >> @:@
        (if <:bool< nops(sysOld) > 1 >> then
          <:par<system of linear equations>>
        else
          <:par<linear equation>>) @:@
        <:par<<:dsymb<
          op({ op(map(x->x=0, select(x-> not type(x,'equation'), sysOld))),
               op(select(x->type(x,'equation'), sysOld)) })
        >>>>
      else
        (if <:bool< maxind = minind >> then
          <:par<The sequence element <:imath< $(symb:u)(minind) >> thus
            satisfies
            <:dsymb<
              op({ op(map(x->x=0, select(x-> not type(x,'equation'), sysOld))),
                   op(select(x->type(x,'equation'), sysOld)) })
            >>
          >>
        else
          <:par<>>))
    else
      (if <:bool< sys <> sysOld >> then
        <:par<Because of the constraints on the initial conditions, >> @:@
        if <:bool< maxind > minind >> then
          <:par<the sequence elements
            <:imath< $(symb:u)(<< minind >>), \dots,~$(symb:u)(<< maxind >>) >>
            satisfies the >> @:@
          (if <:bool< nops(sysOld) > 1 >> then
            <:par<system of linear equations>>
          else
            <:par<linear equation>>) @:@
          <:par<<:dmath<
            << op({ op(map(x->x=0, select(x-> not type(x,'equation'), sysOld))),
                    op(select(x->type(x,'equation'), sysOld)) }) >> .
          >>>>
        else
          (if <:bool< maxind = minind >> then
            <:par<the sequence element <:imath< $(symb:u)(minind) >> satisfies
              <:dsymb<
                op({ op(map(x->x=0,
                            select(x-> not type(x,'equation'), sysOld))),
                     op(select(x->type(x,'equation'), sysOld)) })
              >>
            >>
          else
            <:par<>>)
      else
        <:par<>>) @:@

    (if <:bool< sys <> sysOld >> then
      (if <:bool< nops(sysOld) > 1 >> then
        <:par<The solution of this system is >>
      else
        <:par<The solution of this equation is >>) @:@
      <:par<<:dmath< << op(sys) >> . >>>>
    else
      <:par<>>)
  end,

    (* Clean the _C[] of initial conditions of the type _C[1] +/- _C[2]. *)
  begin
    if <:bool<
      hastype(remove(type,sys,$(u)('anything')='name'),_C['anything'])
      # limitations due to the cost of Groebner basis computation:
      and max(op(map(degree,[seq(op(2,i),i=sys)],
                    indets(sys,_C['anything']))))<3 and not has(sys,'RootOf')
      # also it has to be a system of polynomials
      and type([seq(op(2,i),i=sys)],
                 'list'('polynom'('rational',indets(sys,_C['anything']))))
    >> then
      let _ =
        <:unit<
          inds:=indets(sys,_C['anything']);
          inds2:=sort([op(map(op,indets(sys,$(u)('anything'))))]);
          sys:=subs([seq($(u)(i)=$(u)[i],i=inds2)],
                    {seq(op(1,i)-op(2,i),i=sys)}):
          # find algebraic relations between the u[i]
          termorder:=lexdeg([op(inds)],[seq($(u)[i],i=inds2)]);
          gb:=remove(hastype,Groebner:-Basis(sys,termorder),_C['anything']);
          sol:=subs(solve({op(gb)},{seq($(u)[i],
                    i=inds2[Groebner:-HilbertDimension(
                            gb,termorder)-nops(inds)..-1])}),
                    [seq($(u)[i],i=inds2)]);
          j:=-1;
          for i to nops(sol) do
            if sol[i]=$(u)[inds2[i]] then
              if dorej and nops(select(has,sol,sol[i]))=1 then
                rej[i]:=i
              else
                rej[i]:=NULL;
                j:=j+1;
                sol:=subs($(u)[inds2[i]]=_C[j],sol);
              end if
            end if
          end do;
        >> in
      <<
        {seq($(u)(inds2[i])=sol[i],i={$1..nops(sol)} minus
                {seq(rej[i],i=1..nops(sol))})}
      >>
    else
      << sys >>
  end

  (* Deals with the cases that require/permit no display. *)
  with Return_without_display symbolic_system ->
    <:par<>>, symbolic_system



let diffeqtorec_doit big_r y z u k iniconds yDisplay zDisplay =
   (* := proc(big_r,y,z,u,k,iniconds)
      option `Copyright (c) 1992-2009 by Algorithms Project, INRIA France. \
All rights reserved.`;
     local l, ini, i, rec, j, minordrec, maxordrec, m, r, dr1, inhdeg, inhpart,
       p, rr, cont; *)

  (* SG: added parameters yDisplay, zDisplay. They are used in DynaMoW
   instead of y and z, as diffeqtorec always sets the latter to Y and Z. *)

  (* TODO local variables *)
(*  annotate Maple <* local nn, firstk ; *> ; *)

  <:unit<
    if has($(big_r),$(k)) then
      error $(k),`cannot appear in the differential equation`
    elif has($(big_r),$(u)) then
      error $(u),`cannot appear in the differential equation` end if;

    # initial conditions
    # iniconds is there. SG: originally the test was "nargs=6"
    if $(iniconds)<>NULL then
      l:={seq(op(2,op(0,op(0,i))),
              i=indets([$(iniconds),$(big_r)],
                       `gfun/initeq`($(y)))
           minus {$(y)(0),'D'($(y))(0)})};
      ini:=[$(y)(0)=$(u)(0),
             'D($(y))(0)'=$(u)(1),
             seq(`@@`(D,i)($(y))(0)=$(u)(i)*i!,i=l)];
      r:=subs(ini,$(big_r)); ini:=subs(ini,$(iniconds))
    else r:=$(big_r); ini:={} end if;
  >> ;

  let par1 =
    if <:bool< ini <> {} >> then
      (if <:bool< nops(iniconds) = 1 >> then
        <:par<The initial condition of the differential equation implies>>
      else
        <:par<The initial conditions of the differential equation imply>>) @:@
      <:par<<:dmath< << op(ini) >> . >>>>
    else
      <:par<>>
  in

  <:unit<
    # In very special cases, this loop makes it possible to return
    # an inhomogeneous equation of lower order.
    # Ex: z*(-1+z)^3*(D@@2)(y)(z)+(-1+z)^3*D(y)(z)-(-1+z)^3*y(z)-z*(z-3)
    inhdeg := 0;
    if r[1]<>0 then
      for inhdeg from 0 do
        p:=1-$(z);
        for i from 2 to nops(r) while degree(p,$(z))=1 do
          p:=gcd(p,r[i])
        end do;
        if not has(p,$(z)) then break end if;
        r:=[r[1],seq(quo(r[i],1-$(z),$(z)),i=2..nops(r))]
      end do
    end if;
    # main loop
    ## this fixes a problem with the change in the definition of degree(0)
    ## in V.5
    rr:=map(proc(x) if x=0 then 1 else x end if end,r);
  >> ;

  <:unit<
    eqnWithPowerSeries :=
      add(subs($(z)=$(zDisplay), op(i+2,r)) *
            DDMF:-PowerSeries($(k), i, mul(k-j,j=0..i-1)*$(u)(k),
                  $(zDisplay), -i),
          i=0..nops(r)-2) =
      subs($(z)=$(zDisplay),
           -r[1])/(1-$(zDisplay))^inhdeg
  >> ;

  let inhom_deg = <:int< inhdeg >> in

  let par2 =
    (if inhom_deg > 0 then
      <:par<Dividing by the factor
        <:imath< (1-$(symb:zDisplay))^$(int:inhom_deg) >> and replacing
      >>
    else
      <:par<Replacing >>) @:@
    <:par<<:imath< $(symb:yDisplay)($(symb:zDisplay)) >> by the power series
      <:isymb< DDMF:-PowerSeries($(k),0,$(u)($(k)),$(zDisplay),0) >>
      in the differential equation results in
      <:dmath< << eqnWithPowerSeries >> . >>
    >>
  in

  <:unit<
    minordrec:=min(seq(i-degree(op(i+2,rr),$(z)),i=0..nops(r)-2));
    maxordrec:=max(seq(i-ldegree(op(i+2,rr),$(z)),i=0..nops(r)-2));
    rec:=Array(minordrec..maxordrec,'storage'='sparse');
    for i from 2 to nops(r) do
      for j from ldegree(op(i,r),$(z)) to degree(op(i,r),$(z)) do
        rec[i-2-j]:=rec[i-2-j]+coeff(op(i,r),$(z),j)*
         mul(k+m,m=1-j..i-2-j) end do end do;

    # SG: "expand" removed from main loop, put below
    # SG: added kInit
    # SG: added maxDeg
    # SG: added remPoly (polynomial built from initial power series terms,
    #      which do net reflect the full recurrence)

    maxDeg := max(seq(degree(op(i,r)), i=2..nops(r)));
    kInit := max(0,max(seq(degree(op(i,r))-(i-2), i=2..nops(r))));
    remPoly := collect(add(add(coeff(op(i+2,r),$(z),j) *
                                 add($(u)(kk+i-j) *
     mul(kk+i-j-nn, nn=0..i-1)*$(z)^kk, kk=j..kInit-1), j=0..maxDeg-1),
     i=0..nops(r)-2),$(z));

    # inhomogeneous part of the differential equation
    if r[1]=0 then dr1:=-1; inhpart:=0
    else
      dr1:=degree(r[1],$(z));
      if inhdeg<>0 then # r[1] stands for r[1]/(1-z)^inhdeg
        inhpart := expand(add(coeff(r[1],$(z),i)*mul($(k)+j,
                                                          j=1-i..inhdeg-i-1),
         i=0..dr1))/(inhdeg-1)!;
        r := subsop(1=series(r[1]/(1-$(z))^inhdeg,$(z),
                             max(dr1,maxordrec-minordrec)+1),r)
      else
        inhpart:=0
      end if
    end if;
  >> ;

  (* We now collect powers and shift summation indices. An equation is
     displayed only if it differs from the preceding one. *)

  let par3 =
    if inhom_deg > 0 then
      (let _ =
        (* TODO intermediate step on the rhs *)
        <:unit<
          eqnCollectPowers :=
            add(DDMF:-PowerSeries($(k), kInit+i,
                      subs($(k)=$(k)-i,rec[i])*$(u)(k),
                      $(zDisplay), -i),
                i=minordrec..maxordrec) +
              subs($(z)=$(zDisplay), remPoly) =
            - DDMF:-PowerSeries($(k), 0, inhpart, $(zDisplay), 0) :
        >> in
      if <:bool< eqnCollectPowers <> eqnWithPowerSeries >> then
        <:par<We expand the right-hand side, using the binomial theorem.
          After collecting powers of <:isymb< zDisplay >>
          on the left-hand side, the equation becomes
          <:dmath< << eqnCollectPowers >> . >>
        >>
      else
        <:par<>>) @:@
      (let _ =
        <:unit<
          eqnShift :=
            add(DDMF:-PowerSeries($(k), kInit,
                      rec[i]*$(u)(k+i), $(zDisplay), 0),
                i=minordrec..maxordrec) +
              subs($(z)=$(zDisplay), remPoly) =
            -DDMF:-PowerSeries($(k), 0, inhpart, $(zDisplay), 0):
        >> in
      if <:bool< eqnShift <> eqnCollectPowers >> then
        <:par<Now we shift summation indices, so that the same
          power <:isymb< $(zDisplay)^$(k) >>
          appears in all sums, and we obtain
          <:dsymb< eqnShift >>
        >>
      else
        <:par<>>)
    else
      (let _ =
        <:unit<
          eqnCollectPowers :=
            add(DDMF:-PowerSeries($(k), kInit+i,
                      subs($(k)=$(k)-i, rec[i]) * $(u)(k),
                      $(zDisplay), -i), i=minordrec..maxordrec) +
              subs($(z)=$(zDisplay), remPoly) =
            subs($(z)=$(zDisplay), -r[1]):
        >> in
      if <:bool< eqnCollectPowers <> eqnWithPowerSeries >> then
        <:par<After collecting powers of <:imath< $(symb:zDisplay) >>,
          the equation becomes
          <:dmath< << eqnCollectPowers >> . >>
        >>
      else
        <:par<>>) @:@
      (let _ =
        <:unit<
          eqnShift :=
            add(DDMF:-PowerSeries($(k), kInit,
                      rec[i]*$(u)(k+i), $(zDisplay), 0),
                i=minordrec..maxordrec) +
              subs($(z)=$(zDisplay), remPoly) =
            subs($(z)=$(zDisplay), -r[1]):
        >> in
      if <:bool< eqnShift <> eqnCollectPowers >> then
        <:par<Now we shift summation indices, so that the same
          power <:isymb< $(zDisplay)^$(k) >>
          appears in all sums, and we obtain
          <:dmath< << eqnShift >> . >>
        >>
      else
        <:par<>>)
  in

  <:unit<
    for i from minordrec to maxordrec do rec[i] := expand(rec[i]) end do;
    addini := {};  # addini added by SG

    # SG: Changed "nargs=6" to "iniconds<>NULL".
    if ini<>{} or $(iniconds)<>NULL or dr1<>-1 then
      addini := {op(map(convert,[seq(subs($(k)=i,[coeff(r[1],$(z),i),
        seq(rec[j]*$(u)(i+j),j=max(minordrec,-i)..maxordrec)]),
                                     i=0..dr1)],`+`))};
      for i from dr1+1 while i<-minordrec or subs(k=i,{seq(rec[maxordrec-j],
        j=0..min(i-dr1-1,maxordrec-minordrec))})={0} do
        addini := {op(addini),convert(subs(k=i,[seq(rec[j]*$(u)(i+j),
        j=max(minordrec,-i)..maxordrec),inhpart]),`+`)}
      end do;
    end if;
  >> ;

  let par4 =
    let sOrNo = Wording.ending_of_seq << addini >> in
    <:par<Comparing coefficients of <:isymb< $(zDisplay)^$(k) >>
      gives the >> @:@
    (if <:bool< addini <> {} >> then
      <:par<initial value constraint$(str:sOrNo)
        <:dsymb< op(map(x -> x=0, addini)) >>
        and the >>
    else
      <:par<>>) @:@
    (let _ =
      (* TODO: recurrence holds for smaller values, if inhdeg>0 *)
      (* firstk added by SG; the first k, for which the recurrence holds *)
      <:unit< firstk := max(kInit, dr1 + 1) >> in
    <:par<following recurrence, valid for <:isymb< $(k) >= firstk >>:
      <:dmath<
        << add(rec[i]*$(u)($(k)+i), i=minordrec..maxordrec) + inhpart = 0 >> .
      >>
    >>)
  in

  <:unit<
    ini := {op(ini), op(addini)};
    rec := gfun:-listprimpart(
      subs(k=k-minordrec, [inhpart,seq(rec[i],i=minordrec..maxordrec)]),
      $(k),
      'cont');
  >> ;

  let par5 =
    if <:bool<
         evalb(minordrec <> 0 or (cont <> 1 and not type(cont,'symbol')))
       >>
    then
      (if <:bool<
           evalb(minordrec <> 0 and cont <> 1 and not type(cont,'symbol'))
         >>
      then
        <:par<Cancelling the common
          factor <:isymb< subs($(k)=$(k)+minordrec, cont) >>
          and shifting <:isymb<$(k)>> by <:isymb< -minordrec >>,>>
      else
        if <:bool< cont <> 1 and not type(cont,'symbol') >> then
          <:par<Cancelling the common
            factor <:isymb< subs($(k)=$(k)+minordrec,cont) >>,>>
        else
          <:par<Shifting <:isymb<$(k)>> by <:isymb< -minordrec >>,>>) @:@
      (let _ =
        <:unit< firstk := max(kInit, dr1+1) + minordrec: >> in
      <:par< we find that, for <:isymb< $(k) >= firstk >>,
        <:dmath<
          << add(rec[nn]*$(u)($(k)+nn-2), nn=2..maxordrec-minordrec+2) = 0 >> .
        >>
      >>)
    else
      <:par<>>
  in

  <:unit<
    addini := {};  # added by SG
    # SG: changed "nargs=6" to "iniconds<>NULL"
    if has(cont,$(k)) and (ini<>{} or $(iniconds)<>NULL) then
      l:=-1;
      for j in select(type,gfun:-myisolve(cont,$(k)),nonnegint) do
        if j>=i then # this i is the last one in the previous loop
          l:=l+1;
          addini := addini union {u(j)=_C[l]}
        fi
      od
    fi;
  >> ;

  let par6 =
    if <:bool< addini <> {} >> then
      <:par<Due to the zeros of the canceled factor,
        the following sequence values are arbitrary:
        <:dmath< << addini >> . >>
      >>
    else
      <:par<>>
  in

  <:unit<
    ini := ini union addini;
    while nops(rec)>2 and rec[nops(rec)]=0 do
      rec:=subsop(nops(rec)=NULL,rec) end do;
  >> ;

  (par1 @@@ par2 @@@ par3 @@@ par4 @@@ par5 @@@ par6,
   (* firstk added by SG *)
   << map(collect, rec, $(k)), ini, cont, firstk >>)



(* diffeqtorec
 Input:  eqn: differential equation (for example output of algeqtodiffeq)
    y(z): its unknown function
    u(n): the name of the sequence of Taylor coefficients at the origin
    ini: (optional) boolean indicating whether initial conditions should
     be computed
    content_rec: (optional) will be assigned the content that was removed from
      the recurrence. This is useful when the diff. eqn. is only the homogenous
      part of the actual one.
 Output: the linear recurrence satisfied by u(n) *)
 (* Added by SG: the second member of the output ist the first n for which *)
 (* the recurrence holds. *)

let diffeqtorec eqn yofz uofk ini =

  (* TODO: {ini::boolean:=false,return_content::boolean:=false,
     content_rec::name:=dummy} *)
  (*  diffeqtorec:=proc (eqn,yofz,uofk, {ini::boolean:=false,
     return_content::boolean:=false,content_rec::name:=dummy}) *)
  (*   option `Copyright (c) 1992-2009 by Algorithms Project, INRIA France. \
All rights reserved.`; *)
  (*   local iniconds, f, y, z, u, k, Y, Z, rec, inirec, cont; *)

  (* retVal added by SG *)
(*  symbolic retVal ; *)

  (* TODO remove, as soon as iniconds is local *)
  <:unit< unassign('iniconds'): >> ;

  let inirec = ref << "empty" >> in

  <:unit<
    # This also checks that >=3 args are passed.
    gfun:-getname($(uofk),'u','k'):
  >> ;

  if <:bool< $(bool:ini) or type($(eqn),'set') >> then
    <:unit< f := gfun:-formatdiffeq([$(eqn),$(yofz)],'y','z','iniconds'): >>
  else
    <:unit< f := gfun:-formatdiffeq([$(eqn),$(yofz)],'y','z'): >> ;

  let par1 =

    <:par<The function under consideration satisfies the differential equation
      <:dsymb< add(op(i+2,f)*diff(y(z),[z$i]), i=0..nops(f)-2) = -op(1,f) >>
    >> @:@

    (let sOrNo = Wording.ending_of_seq << iniconds >> in
    if <:bool< type(iniconds, 'set') and not iniconds = {} >> then
      <:par<with initial condition$(str:sOrNo)
        <:dmath< << op(iniconds) >> . >>
      >>
    else
      <:par<>>) @:@
    <:par<Our goal is to find a recurrence for the Taylor
      coefficients <:isymb< $(uofk) >>
      of <:isymb< $(yofz) >>.>>

  in

  (* Avoid problems when u=y or u=z, ... *)
  <:unit< f := subs([y=Y,z=Z], f): >> ;

  if <:bool< $(bool:ini) or type($(eqn),'set') >> then
    <:unit< iniconds := subs(y=Y, iniconds): >>
  else
    <:unit< iniconds := NULL: >> ;

  (* SG: added d2recReturn; TODO should become local *)
  let par2, d2recReturn =
    diffeqtorec_doit
      << f >>
      << Y >>
      << Z >>
      << u >>
      << k >>
      << iniconds >>
      << y >>
      << z >> in

  <:unit< rec := op(1,[$(d2recReturn)]): >> ;
  inirec := << op(2,[$(d2recReturn)]) >> ;
  <:unit<
    cont := op(3,[$(d2recReturn)]):
    firstk := op(4,[$(d2recReturn)]):
  >> ;

  let par3 =
    if <:bool< iniconds<>{} or iniconds<>NULL: >> then
      let p, ir =
        goodinitvalues_rec
          << rec >>
          << u >>
          << k >>
          !inirec
          true
          << NULL >> in
      begin
        inirec := ir ;
        p
      end
    else
      <:par<>>
  in

  (* SG: Reassign u and k; this is a work-around. So far there is no way *)
  (* to have local variables, hence k may have been erroneously modified *)
  (* by a procedure call. *)
  <:unit< gfun:-getname($(uofk),'u','k'): >> ;

  let result, par4 =
    (* if added by SG *)
    (* FIXME: FC, I changed the coding-style so that not a false is returned *)
    (* but an exception is raised. *)
    let the_inirec = !inirec in
    if <:bool< $(the_inirec) = false >> then
      << false, infinity >>,
      <:par<The initial conditions are contradictory, hence
        the differential equation above has no power series >> @:@
      (if <:bool< type(iniconds,'set') and iniconds<>{} >> then
        <:par<solution with the given initial values.>>
      else
        <:par<solution.>>)
    else
      let r = << gfun:-makerec(rec,u,k,$(the_inirec)), firstk >> in
      r,
      if <:bool< type(op(1, [$(r)]), 'set') and nops(op(1,[$(r)])) > 1 >>
      then
        let the_inirec = << firstk >> in
        <:par<Summarizing, the coefficients <:imath< u(k) >>
          of any power series solution of the differential equation above
          satisfy, for <:imath< k \geq $(symb:the_inirec) >>,
          the recurrence relation
          <:dsymb< op(1,op(1,[$(r)])) = 0 >>
        >> @:@
        (* TODO: Give references to the steps where the initial conditions *)
        (* occur. *)
        (let sOrNo =
          Wording.ending_of_int (<:int< nops(op(1,[$(r)])) >> - 1) in
        <:par<with initial condition$(str:sOrNo)
          <:dmath< << op(2..nops(op(1,[$(r)])), op(1,[$(r)])) >> . >>
        >>)
      else <:par<>> in
(result, par1 @@@ par2 @@@ par3 @@@ par4)


let_service DiffeqToRecurrence
  (eqn : diffeq maple)
  (yofz : fcall maple)
  (uofk : fcall maple)
  (ini : bool)
  (notation : string) :
  DC.sec_entities * unit with { title = title } =

  let title = title _req_params
  and rec_eq, pars = diffeqtorec eqn yofz uofk ini in

(* Not used because of following FIXME.
  (* The second member of diffeqtorec's return value ist firstk,
    the index from which on the recurrence holds. *)
  let rec_eq = << op(1,[$(symb:rec_eq)]) >> in
*)

  (* FIXME: This has to become a test file. *)
  (* FIXME: Such comments, <!-- ... --> cannot be expressed. *)
  (*
  if <:bool< $(symb:rec_eq) = $(symb:desired_result) >> then
    [t: <!-- diffeqtorec passed the test for sf_id=$(symb:sf_id)
      and eqn=$(symb:eqn). --> :t]
  else
    [t: <p class='warning'>Proof failed:
           The return value of diffeqtorec_ddmf_wrapper<br/>$(symb:rec)<br/>
      does not match the desired result<br/>$(symb:desired_result). :t]
  endif ;
  *)

  (DC.section title pars, ())

Generated by GNU Enscript 1.6.5.90.