%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Computer program written in Mathematica
to calculate the spherical growth series 
of amalgamated free products of infinite cyclic groups

2023-01-25 by M. Fujii and T. Sakasai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


growth[x__] := Module[{p, n, q, r},
  p = Sort[{x}];
  n = Length[p];
  
  q = {};
  r = {};
  Do[If[EvenQ[p[[i]] - p[[1]]], q = Append[q, p[[i]]], 
    r = Append[r, p[[i]]]], {i, 1, n}];
  
  If[Length[q] < n,
    2*h[Flatten[
         Table[{addminus[p[[1]], p[[i]]], 
           addplus[p[[1]], p[[i]]]}, {i, 1, n}]]]/(1 - t^(p[[1]])) - 
     h[Flatten[
       Table[{addminus[p[[1]], p[[i]]], addminus[p[[1]], p[[i]]]}, {i,
          1, n}]]] + growth3zero[q, r], 
    2*h[Flatten[
         Table[{addminus[p[[1]], p[[i]]], 
           addplus[p[[1]], p[[i]]]}, {i, 1, n}]]]/(1 - t^(p[[1]])) - 
     h[Flatten[
       Table[{addminus[p[[1]], p[[i]]], addminus[p[[1]], p[[i]]]}, {i,
          1, n}]]] + growth3zeroeven[q]] // Factor
  ];

f[u_, v_] := Sum[t^i, {i, 1, u}] + Sum[t^i, {i, 1, v}];

g[list_] := Module[{i, j},
  1/(1 - Sum[(i - 1)*
       SymmetricPolynomial[i, 
        Table[f[list[[2 j - 1]], list[[2 j]]], {j, 1, 
          Length[list]/2}]], {i, 2, Length[list]/2}])
  ];

h[list_] := Module[{j},
  Product[
    1 + f[list[[2 j - 1]], list[[2 j]]], {j, 1, Length[list]/2}]*
   g[list]
  ];

addminus[pmin_, x_] := Floor[(x - pmin)/2];

addplus[pmin_, x_] := Floor[(pmin + x - 1)/2];

growth3zero[qlist_, rlist_] := Module[{q, r, en, em},
  q = qlist;
  r = rlist;
  
  Sum[Sum[growth30part1[q, r, en, em], {em, 0, Length[r] - 1}], {en, 
     0, q[[1]] - 1}] + 
   Sum[Sum[growth30part23[q, r, en, em], {em, 0, Length[r] - 2}], {en,
      0, q[[1]] - 1}] + 
   Sum[growth30part23[q, r, en, Length[r] - 1], {en, 0, q[[1]] - 2}] +
    Sum[Sum[
     growth30part4[q, r, en, em], {em, Length[r], 
      Length[q] + Length[r] - 1}], {en, 0, q[[1]] - 2}] +
   Sum[Sum[
     growth30part5[q, r, en, em], {em, Length[r], 
      Length[q] + Length[r] - 1}], {en, 0, q[[1]] - 2}]
  
  ];

growth30part1[qlist_, rlist_, numN_, numM_] := 
 Module[{q, r, n, m, en, em, p1, qminus, qplus, rminus, rplus, k, j, 
   i},
  q = qlist;
  r = rlist;
  m = Length[q];
  n = m + Length[r];
  en = numN;
  em = numM;
  p1 = q[[1]];
  
  qminus = Map[addminus[p1, #] &, q];
  qplus = Map[addplus[p1, #] &, q];
  rminus = Map[addminus[p1, #] &, r];
  rplus = Map[addplus[p1, #] &, r];
  
  (t^(r[[n - em - m]]))*
   Product[1 + f[qminus[[k]] + en, qplus[[k]] - en], {k, 1, m}]*
   Product[1 + f[rminus[[k - m]] + en, rplus[[k - m]] - en], {k, 
     m + 1, n - em - 1}]*
   Product[1 + 
     f[rminus[[k - m]] + en + 1, rplus[[k - m]] - en - 1], {k, 
     n - em + 1, n}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, m}], 
      Table[{rminus[[j - m]] + en, rplus[[j - m]] - en}, {j, m + 1, 
        n - em - 1}],
      {rminus[[n - em - m]] + en, rplus[[n - em - m]] - en - 1}, 
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        n - em + 1, n}]]]]*
   Sum[i*SymmetricPolynomial[i,
      Join[Table[f[qminus[[j]] + en, qplus[[j]] - en], {j, 1, m}], 
       Table[f[rminus[[j - m]] + en, rplus[[j - m]] - en], {j, m + 1, 
         n - em - 1}], {0}, 
       Table[f[rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1], {j,
          n - em + 1, n}]]], {i, 1, n - 1}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, m}], 
      Table[{rminus[[j - m]] + en, rplus[[j - m]] - en}, {j, m + 1, 
        n - em - 1}], 
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        n - em, n}]]]]*
   Product[1 + f[qminus[[k]] + en, qplus[[k]] - en], {k, 1, m}]*
   Product[1 + f[rminus[[k - m]] + en, rplus[[k - m]] - en], {k, 
     m + 1, n - em - 1}]*
   Product[1 + 
     f[rminus[[k - m]] + en + 1, rplus[[k - m]] - en - 1], {k, 
     n - em + 1, n}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, m}], 
      Table[{rminus[[j - m]] + en, rplus[[j - m]] - en}, {j, m + 1, 
        n - em}], 
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        n - em + 1, n}]]]]
  ];

growth30part23[qlist_, rlist_, numN_, numM_] := 
 Module[{q, r, n, m, en, em, p1, qminus, qplus, rminus, rplus, k, j, 
   i},
  q = qlist;
  r = rlist;
  m = Length[q];
  n = m + Length[r];
  en = numN;
  em = numM;
  p1 = q[[1]];
  
  qminus = Map[addminus[p1, #] &, q];
  qplus = Map[addplus[p1, #] &, q];
  rminus = Map[addminus[p1, #] &, r];
  rplus = Map[addplus[p1, #] &, r];
  
  h[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, m}], 
      Table[{rminus[[j - m]] + en, rplus[[j - m]] - en}, {j, m + 1, 
        n - em - 1}],
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        n - em, n}]]]] - 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, m}], 
      Table[{rminus[[j - m]] + en, rplus[[j - m]] - en}, {j, m + 1, 
        n - em - 1}],
      {rminus[[n - em - m]] + en, rplus[[n - em - m]] - en - 1},
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        n - em + 1, n}]]]] - 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qminus[[j]]}, {j, 1, m}], 
      Table[{rminus[[j - m]] + en, rminus[[j - m]]}, {j, m + 1, 
        n - em - 1}],
      Table[{rminus[[j - m]] + en + 1, rminus[[j - m]]}, {j, n - em, 
        n}]]]] + 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qminus[[j]]}, {j, 1, m}], 
      Table[{rminus[[j - m]] + en, rminus[[j - m]]}, {j, m + 1, 
        n - em}],
      Table[{rminus[[j - m]] + en + 1, rminus[[j - m]]}, {j, 
        n - em + 1, n}]]]]
  ];

growth30part4[qlist_, rlist_, numN_, numM_] := 
 Module[{q, r, n, m, en, em, p1, qminus, qplus, rminus, rplus, k, j, 
   i},
  q = qlist;
  r = rlist;
  m = Length[q];
  n = m + Length[r];
  en = numN;
  em = numM;
  p1 = q[[1]];
  
  qminus = Map[addminus[p1, #] &, q];
  qplus = Map[addplus[p1, #] &, q];
  rminus = Map[addminus[p1, #] &, r];
  rplus = Map[addplus[p1, #] &, r];
  
  (t^(q[[n - em]]))*
   Product[
    1 + f[qminus[[k]] + en, qplus[[k]] - en], {k, 1, n - em - 1}]*
   Product[1 + f[qminus[[k]] + en + 1, qplus[[k]] - en - 1], {k, 
     n - em + 1, m}]*
   Product[1 + 
     f[rminus[[k - m]] + en + 1, rplus[[k - m]] - en - 1], {k, m + 1, 
     n}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}],
      {qminus[[n - em]] + en, qplus[[n - em]] - en - 1},
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, 
        n - em + 1, m}],
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        m + 1, n}]]]]*
   Sum[i*SymmetricPolynomial[i,
      Join[
       Table[f[qminus[[j]] + en, qplus[[j]] - en], {j, 1, 
         n - em - 1}], {0}, 
       Table[f[qminus[[j]] + en + 1, qplus[[j]] - en - 1], {j, 
         n - em + 1, m}], 
       Table[f[rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1], {j,
          m + 1, n}]]], {i, 1, n - 1}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}], 
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, n - em, 
        m}], Table[{rminus[[j - m]] + en + 1, 
        rplus[[j - m]] - en - 1}, {j, m + 1, n}]]]]*
   Product[
    1 + f[qminus[[k]] + en, qplus[[k]] - en], {k, 1, n - em - 1}]*
   Product[1 + f[qminus[[k]] + en + 1, qplus[[k]] - en - 1], {k, 
     n - em + 1, m}]*
   Product[1 + 
     f[rminus[[k - m]] + en + 1, rplus[[k - m]] - en - 1], {k, m + 1, 
     n}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, n - em}], 
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, 
        n - em + 1, m}],
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        m + 1, n}]]]]
  ];

growth30part5[qlist_, rlist_, numN_, numM_] := 
 Module[{q, r, n, m, en, em, p1, qminus, qplus, rminus, rplus, k, j, 
   i},
  q = qlist;
  r = rlist;
  m = Length[q];
  n = m + Length[r];
  en = numN;
  em = numM;
  p1 = q[[1]];
  
  qminus = Map[addminus[p1, #] &, q];
  qplus = Map[addplus[p1, #] &, q];
  rminus = Map[addminus[p1, #] &, r];
  rplus = Map[addplus[p1, #] &, r];
  
  h[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}], 
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, n - em, 
        m}],
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        m + 1, n}]]]] - 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}],
      {qminus[[n - em]] + en, qplus[[n - em]] - en - 1},
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, 
        n - em + 1, m}],
      Table[{rminus[[j - m]] + en + 1, rplus[[j - m]] - en - 1}, {j, 
        m + 1, n}]]]] - 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qminus[[j]]}, {j, 1, n - em - 1}],
      Table[{qminus[[j]] + en + 1, qminus[[j]]}, {j, n - em, m}],
      Table[{rminus[[j - m]] + en + 1, rminus[[j - m]]}, {j, m + 1, 
        n}]]]] + 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qminus[[j]]}, {j, 1, n - em}], 
      Table[{qminus[[j]] + en + 1, qminus[[j]]}, {j, n - em + 1, m}],
      Table[{rminus[[j - m]] + en + 1, rminus[[j - m]]}, {j, m + 1, 
        n}]]]]
  ];

growth3zeroeven[qlist_] := Module[{q, en, em},
  q = qlist;
  
  Sum[Sum[growth30evenpart1[q, en, em], {em, 0, Length[q] - 1}], {en, 
     0, q[[1]] - 2}] + 
   Sum[Sum[growth30evenpart23[q, en, em], {em, 0, 
      Length[q] - 2}], {en, 0, q[[1]] - 2}] + 
   Sum[growth30evenpart23[q, en, Length[q] - 1], {en, 0, q[[1]] - 3}]
  ];

growth30evenpart1[qlist_, numN_, numM_] := 
 Module[{q, n, en, em, p1, qminus, qplus, k, j, i},
  q = qlist;
  n = Length[q];
  en = numN;
  em = numM;
  p1 = q[[1]];
  
  qminus = Map[addminus[p1, #] &, q];
  qplus = Map[addplus[p1, #] &, q];
  
  
  (t^(q[[n - em]]))*
   Product[
    1 + f[qminus[[k]] + en, qplus[[k]] - en], {k, 1, n - em - 1}]*
   Product[1 + f[qminus[[k]] + en + 1, qplus[[k]] - en - 1], {k, 
     n - em + 1, n}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}],
      {qminus[[n - em]] + en, qplus[[n - em]] - en - 1}, 
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, 
        n - em + 1, n}]]]]*
   Sum[i*SymmetricPolynomial[i,
      Join[
       Table[f[qminus[[j]] + en, qplus[[j]] - en], {j, 1, 
         n - em - 1}], {0}, 
       Table[f[qminus[[j]] + en + 1, qplus[[j]] - en - 1], {j, 
         n - em + 1, n}]]], {i, 1, n - 1}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}], 
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, n - em, 
        n}]]]]*
   Product[
    1 + f[qminus[[k]] + en, qplus[[k]] - en], {k, 1, n - em - 1}]*
   Product[1 + f[qminus[[k]] + en + 1, qplus[[k]] - en - 1], {k, 
     n - em + 1, n}]*
   g[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, n - em}], 
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, 
        n - em + 1, n}]]]]
  ];

growth30evenpart23[qlist_, numN_, numM_] := 
 Module[{q, n, en, em, p1, qminus, qplus, k, j, i},
  q = qlist;
  n = Length[q];
  
  en = numN;
  em = numM;
  p1 = q[[1]];
  
  qminus = Map[addminus[p1, #] &, q];
  qplus = Map[addplus[p1, #] &, q];
  
  
  h[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}], 
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, n - em, 
        n}]]]] - 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qplus[[j]] - en}, {j, 1, 
        n - em - 1}],
      {qminus[[n - em]] + en, qplus[[n - em]] - en - 1},
      Table[{qminus[[j]] + en + 1, qplus[[j]] - en - 1}, {j, 
        n - em + 1, n}]]]] - 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qminus[[j]]}, {j, 1, n - em - 1}],
      Table[{qminus[[j]] + en + 1, qminus[[j]]}, {j, n - em, n}]]]] + 
   h[Flatten[
     Join[Table[{qminus[[j]] + en, qminus[[j]]}, {j, 1, n - em}], 
      Table[{qminus[[j]] + en + 1, qminus[[j]]}, {j, n - em + 1, n}]]]]
  ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Example
< x_1,x_2,x_3,x_4 | (x_1)^2=(x_2)^3=(x_3)^7=(x_4)^8 >

[In]
growth[2, 3, 7, 8]
[Out]
(1 + 4 t - 16 t^2 - 163 t^3 - 516 t^4 - 320 t^5 + 4296 t^6 + 
   24213 t^7 + 81073 t^8 + 206772 t^9 + 434218 t^10 + 778907 t^11 + 
   1218441 t^12 + 1683284 t^13 + 2070095 t^14 + 2277157 t^15 + 
   2246173 t^16 + 1987792 t^17 + 1576248 t^18 + 1116383 t^19 + 
   702382 t^20 + 389314 t^21 + 187754 t^22 + 77324 t^23 + 
   26408 t^24 + 7112 t^25 + 1368 t^26 + 
   144 t^27)/((-1 + t) (-1 + t + 17 t^2 + 43 t^3 + 69 t^4 + 80 t^5 + 
     62 t^6 + 36 t^7 + 12 t^8)^2 (-1 + t + 12 t^2 + 33 t^3 + 55 t^4 + 
     68 t^5 + 64 t^6 + 47 t^7 + 26 t^8 + 11 t^9 + 3 t^10))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%