(* ::Package:: *)

zmult[v_, l_, s_, prec_] := 
  Block[{z1, q, t, ze, b, h, r, j, z, m, n}, (
    t = Length[v];
    ze = Table[1, {s}];
    b = Table[0, {l}];
    
    (* prec = Floor[Log[s]*Sqrt[s]]; *)
    
    b[[1]] = SetPrecision[ 1.0, prec];
    h = -1;
    
    For[r = 1, r <= t, r++,
     j = v[[t + 1 - r]];
     h += j - 1;
     z = SetPrecision[0.0, prec];
     q = SetPrecision[1.0, prec] / s^h;
     
     For[m = 1, m <= l, m++,
      	q /= s;
      	b[[
        m]] = (b[[m]] + 
          Sum[-Binomial[h + m, m - n + 1] * b[[n]], {n, 1, 
            m - 1}])/(h + m);
      	z += b[[m]] * q;
      	];
     
     For[m = s, m >= 1, m--,
      	z1 = z + ze[[m]] / m^j;
      	ze[[m]] = z;
      	z = z1;
      	];
     ];
    Return[z];
    )];

ztmult[v_, t_, M_, l_, prec_] := 
  Block[{k, z, a, b, w, h, u, n , ze , ze0, i, j, b0}, (
    n = Length[v];
    a = Table[0, {n + 1}, {l}];
    b = Table[0, {l + 1}];
    
    For[i = 1, i <= n, i++,
     h = v[[n + 1 - i]];
     If[i == 1,
      b[[h]] = 1;,
      For[j = 1 + l, j >= 1, j--,
        b[[j]] = 
          If[j >= h + 1, a[[i, j - h]] + t*b[[j - h]], 0] // Expand //
            SetPrecision[#, prec] &;
        ];
      ];
     For[j = 1, j <= l , j++,
      a[[i + 1, 
          j]] = ((b[[j + 1]] + 
              Sum[-Binomial[j, u - 1] * a[[i + 1, u]], {u, 1, 
                j - 1}]) / j) // Expand // SetPrecision[#, prec] &;
      ];
     ];
    
    ze = Table[SetPrecision[0.0, prec], {n + 1}];
    ze[[n + 1]]  =  SetPrecision[1.0, prec];
    ze0 = ze;
    
    For[i = 1, i <= n, i++,
     ze[[i]] = 
       Sum[a[[n + 2 - i, j]]/M^j, {j, 1, l}] // Expand // 
        SetPrecision[#, prec] &;
     ];
    
    For[i = M, i >= 1, i--,
     b0 =  SetPrecision[0.0, prec];
     For[j = n, j >= 1, j--,
      b0 = (t *b0 + ze[[j + 1]]) / i^v[[j]] // Expand // 
        SetPrecision[#, prec] &;
      ze0[[j]] = ze[[j]] + b0 // Expand // SetPrecision[#, prec] &;
      ];
     ze = ze0 // Expand // SetPrecision[#, prec] &;
     ];
    
    Return[ze[[1]]];
    )];
    
(* --------------------- *)

(* Memoisation *)
Clear[mzv];
Clear[mzvt];
mzv[v_, s_, prec_] := 
  mzv[v, s, prec] = zmult[v, Floor[Sqrt[s*Log[s]/2]], s, prec];
mzvt[v_, t_, s_, prec_] := 
  mzvt[v, t, s, prec] = 
   ztmult[v, t, s, Floor[Sqrt[s*Log[s]/2]], prec];
   
(* (* No memoisation *)
Clear[mzv];
Clear[mzvt];
mzv[v_, s_, prec_] := zmult[v,Floor[Sqrt[s*Log[s]/2]],s, prec];
mzvt[v_, t_, s_, prec_] := ztmult[v,t,s, \
Floor[Sqrt[s*Log[s]/2]],prec]; *)

ztnaive[v_, t_] := Block[{zTemp, i}, (
    Plus @@ (Composition[
         t^(1 + Length[v] - Length[#]) * 
           zTemp @@ 
            Table[Plus @@ v[[#[[i]] + 1 ;; #[[i + 1]]]], {i, 1, 
              Length[#] - 1}] &, 
         FoldList[Plus, 0, #] &] /@ (Flatten[
          Permutations /@ IntegerPartitions[Length[v]], 
          1])) /. {zTemp[expr___] :> z[{expr}]}
    )];

z[v_] := mzv[v, 10000, 500];
zt[v_,t_] := mzvt[v, t, 10000, 500];
