(* :Title: Differential Dynamical Systems      *)
(* :Name: Iterate`                             *)
(* :Author: Daniel Geisler, Dec 2003.          *)
(* :Summary:                                   *)
(* :Context: Iterate`                          *)
(* :Package Version: 0.3                       *)
(* :Copyright: Copyright 2003, Daniel Geisler. *)
(* :Mathematica Version: 4.1                   *)


BeginPackage["Iterate`"]
Off[General::"spell1"];

(* Bad sign; the Mathematica code is likely implemented in an inefficient manner. *)  
$RecursionLimit=Infinity;  

(* User functions *)
Iterate::usage = "Iterate[function, time variable, space variable, fixed point, derivatives computed]; expample Iterate[f, n, z, p, size]. Computers the continuous iteration of f at fixed point p."
Bell::Usage = "Bell[n] is the nth Bell polynomial."
dyne::Usage = "dyne[n] the internal hybrid analytic/combinatoric representation of the nth derivative of iterated function."
Dyne::Usage = "The nth derivative of iterated function."
hierarchies::Usage = "hierarchies[n] computes the instances of Schroeder's Fourth Problem for n items."

Hyperbolic::Usage = ""
Parabolic::Usage = ""
Neutral::Usage = ""
Superattracting::Usage = ""


(* The following are internal functions exported to properly display internal structures. *)
dyn::Usage = "Instance of the unlabeled hierarchies combinatoric structure."
k::Usage = "Used to manage summation iterators."
d::Usage = "d[k] is used to represent the kth derivative of the current iterated function."
hier::Usage = "Instance of the hierarchies combinatoric structure."


Begin["`Private`"]

Options[Iterate] =  {Verbose -> False};

Format[k[_,i_],TeXForm]:=Subscript[k,i]; 
Format[d[i_],TeXForm]:=Subscript[$function,i];
Format[dyn[i_],TeXForm]:="{dyn("<>ToString[TeXForm[i]]<>")}";
Format[Derivative[i_][f_][_],TeXForm]:=Subscript[f,i];

Format[k[_,i_]]:=Subscript[k,i]; 
Format[d[i_]]:=Subscript[$function,i]; 


(* Bell polynomials -----------------------------------------------------*) 
  Bell[0]=f[g[z]];
  Bell[i_Integer/;i>0] := Bell[i]=D[Bell[i-1],z];
(* Bell polynomials - End *)

(* Iterate -----------------------------------------------------   *)
(*                                                                   *)
  (* d[k] is the kth derivative of f(z_0) where z_0 is a fixed point *)
  (* Dyne[k] is the kth derivative of dynamical system f^n(z_0) while*)
  (* dyne[k] is the quasi-combinatorial internal version of Dyne[k].  *)  
  (* The following is classical identity of Lyapunov exponent as     *)
  (* first derivative of dynamical system.                           *)
  dyne[1]=d[1]^$n; 

  (* Done manually because of limitation of RSolve function in Mathematica 4.1         *)
  (* Solves recurrence equation for kth derivative of dynamical system f^n(z_0)        *)
  (* dyne[2]=d[2]*dyn[{2}]                                                             *)  
  (* dyne[3]=3*d[2]^2*dyn[{3, {2}}] + d[3]*dyn[{3}]                                    *)      
  (* The Recurrence[] function computes the characteristic equation used to solve the  *)
  (* recurrence equation for kth derivative of dynamical system f^n(z_0). The solve[]  *)
  (* function actually finds the solution of the characteristic equation.              *) 
  (* The dyn[] function is a representation of series-reduced planted trees with n     *)
  (* leaves and unlabeled connected cographs on n nodes - EIS A000669; also unlabeled  *)   
  (* hierarchies - Flajolet. The series begins 1, 1, 2, 5, 12, 33, 90, 261, 766, 2312. *)  
  (* The connection with Schroeder's Fourth Problem EIS A000311 can be seen by summing *)     
  (* the coefficients. The first number in the bracket is the total number of nodes    *)
  (* within the brackets with partitions possibly following. The expression 3*d[2]^2   *)
  (* from dyne[3] could be derived from dyn[{3, {2}}] but the information is already   *)
  (* computed. An earlier version of this program directly calculated the nested       *)
  (* summations but it became easier to reduce the nested summations to their          *)
  (* combinatoric basis, and derived the rules that transformed the multiplication and *) 
  (* powers of the nested summations into the generation of a combinatoric structure   *)
  (* from simpler structures. The Decode[] function translates the internal            *)
  (* quasi-combinatorial structure into the associated nested summation.               *)
  
  Recurrence::usage ="Recurrence[n] computes the recurrence equation."
  Recurrence[i_/;i>0]:= 
    Recurrence[i]=
      ExpandAll[
        Bell[i] /. D[g[z],{z,i}] \[Rule] 0  
                /. g[z] \[Rule] z 
                /. Derivative[m_][g][z] \[RuleDelayed] dyne[m] /. $n \[Rule] $n-1 
                /. Derivative[m_][f][z] \[Rule] d[m]
    ];

  (* These rules are probably why $RecursionLimit=Infinity must be used *)
  (* They are the combinatoric analogs of multiplying nested summations  *)
  (* together and raising them to an integer power. *)     
  dyn /: dyn[{a_,b__}]^q_ := dyn[Join[{a},Table[b,{j,q}]]] /; a \[Equal] $derivative; 
  dyn /: dyn[{a_,b__}]^q_ := dyn[Table[{a,b},{j,q}]] /; a \[NotEqual] $derivative;   
  dyn /: dyn[{a_,b__}]*dyn[{a_,c__}] := dyn[{a,b,c}] /; a \[Equal] $derivative;
  dyn /: dyn[{a_,b__}]*dyn[{a_,c__}] := dyn[{{a,b},{a,c}}] /; a \[NotEqual] $derivative;
  dyn[lst:{a_,{b__}}] := dyn[{a,b}] /; ¬IntegerQ[lst[[2]][[1]]];

  (* The solve[] function finds the solution of the characteristic equation. *)
  solve[x_,i_]:= 
      If[FreeQ[x,dyn[__]], 
        (x /. d[1] \[Rule] 1)*dyn[{i}],
         x /. d[1] \[Rule] 1 /. dyn[a_] \[Rule] dyn[{i,a}]
      ];

  solve[a_+b_,i_] := solve[a,i]+solve[b,i];


  (* Computes the exponential of the nested summation *)
  exp[x_]:=d[1]^(First[First[x]]*$n-Plus @@ (First /@ Flatten[x])
                 -Plus @@((First[#]-1)*# &)/@ Flatten[x]);

  (* Recursively computes the range constraints of the nested summation's iterator *)
  prop[{{k[a__],0,b_}},d_] := {{k[a],0,b-d}};
  prop[{{k[a__],0,b_},c__},d_] := {{k[a],0,b-d},(prop[#,1+d+k[a]]&)/@{c}};

  (* Computes the nested summation's iterator *)
  iter[x_] := Partition[Flatten[prop[x /. k[a__] \[Rule] {k[a],0,$n-1},0]],3];

  (* This is an analytic functor transforming the combinatoric structure represented *)
  (* by dyn[] into the associated nested summation.                                  *)
  Dyn[x_] := Module[{position=1,y},
          y = x /. a_Integer \[RuleDelayed] k[a,position++];
          sum[exp[y],iter[y]] 
        ];

  (* dyne[k] is the quasi-combinatorial internal version of the kth derivative of *)
  (* the dynamical system f^n(z_0). *)  
  dyne[i_] := dyne[i]=solve[Recurrence[i],i];

  (* Dyne[k] is the kth derivative of dynamical system f^n(z_0) in the form of    *)
  (* nested summations. *)
  Dyne[i_] := Dyne[i]=ReleaseHold[Decode[dyne[i]]];

  (* Decodes from dyne[k] to Dyne[k]. *)
  Decode::usage ="Decode[expression] converts expression from an internal format."
  Decode[x_] :=
      Module[{},
          x (* /. Times \[Rule] NonCommutativeMultiply *)
            /. dyn \[Rule] Dyn
            /. d[m_] \[RuleDelayed] derv[m] 
            /. sum[a_,{b__}] \[Rule] HoldForm[Sum[a,b]]
      ];             


   Hyperbolic[i_] := Hyperbolic[i] = ExpandAll[Dyne[i] /. Derivative[1][$function][$p] \[Rule] \[Lambda] /. \[Lambda] \[Rule] Derivative[1][$function][$p]];
   Parabolic[i_] := Parabolic[i] = ExpandAll[Dyne[i] /. Derivative[1][$function][$p] \[Rule] 1];
   Neutral[i_] := Neutral[i] = Hyperbolic[i] /. n \[Rule] Mod[n,j];
   sup[x_] := Module[{position=1,y},
         y = x /. a_Integer \[RuleDelayed] k[a,position++];
         KroneckerDelta[(First[First[y]]*$n-Plus @@ (First /@ Flatten[y])-
           Plus @@((First[#]-1)*# &)/@ Flatten[y])]/.k[__]\[Rule]0
       ];
   Superattracting[i_] := Superattracting[i]=dyne[i] /. dyn \[Rule] sup /. d[m_] \[Rule] Derivative[m][$function][$p] ;


  (* Main function of package *)
  (* The scope of $derivative and $n is the package so that the dyn rules will work *)  
  (* The scope of $function, $p, and $z is the package to support different formats *)  
  Iterate[f_, n_, z_, p_, max_Integer:4 ,opts___] := 
    Module[{verbose},
      {verbose} = {Verbose} /.{opts} /.Options[Iterate];
      If[verbose, 
        Print["----------------------------------------------------"];
      ];
      $n=n;
      $function=f;
      $p=p;
      $z=z;
      derv[1]=D[f[z],{z,1}] /. z \[Rule] p;
      Do[
        derv[$derivative]=D[f[z],{z,$derivative}] /. z \[Rule] p;
        If[verbose, 
          Print[Superscript[D,$derivative] Superscript[f,n][p]//TraditionalForm];   
        ];        
        dyne[$derivative];
        If[verbose, 
          Print[dyne[$derivative] /. d[k_] \[Rule] Subscript[f,k] //TraditionalForm];
          Print[Dyne[$derivative]//TraditionalForm];
          Print["----------------------------------------------------"];
        ];          
      ,{$derivative,2,max}
      ];
      p+Sum[1/$derivative!*Dyne[$derivative]*(z-p)^$derivative,{$derivative,1,max}]
    ];

(* Iterate - End*)


(* Combinatorics -----------------------------------------------------*)

  (* Following use pointing operator to create a hierarchies version of *)
  (* use of pointing operator in set partitions *)

  (* lead partition functor *)
  lead[_, n_ /; n < 3] := 0

  lead[h_, n_] := Module[{p,i},
      p = Position[h, {___}];
      Sum[MapAt[{#, n} &, h, p[[i]]], {i, Length[p]}]
      ]

  (* append partition functor *)
  follow[h_, n_] := Module[{r,i},
      r = Replace[Position[h, {___}], {a__} \[Rule] {a, -1}, 1];
      Sum[Insert[h, n, r[[i]]], {i, Length[r]}]
      ]

  (* join partition functor *)
  marry[_, n_ /; n < 3] := 0

  marry[h_, n_] := Module[{p,i},
      p = Position[h, _Integer];
      Sum[MapAt[{#, n} &, h, p[[i]]], {i, Length[p]}]
      ]

  (* Combinatoric differentiation *)
  extend[a_ + b_, n_] := extend[a, n] + extend[b, n]

  extend[a_, n_] := lead[a, n] + follow[a, n] + marry[a, n]

  (* Seed for hierarchies *)
  hierarchies[1] := hierarchies[1] = extend[hier[{}], 1]

  (* Recursive generation of hierarchies structure *)
  hierarchies[n_] := hierarchies[n] = extend[hierarchies[n - 1], n]

(* Combinatorics - End*)

End[]

EndPackage[]