BeginPackage["Hosoya`"];

  Needs["DiscreteMath`RSolve`"];

  Off[General::spell1, General::spell];

  alpha::usage = "alpha[s, x] where s is a string over {1,2,3} gives the \
Hosoya polynomial of the hexagonal chain consisting of Length[s] hexagons \
where s defines the type (1, 2, or 3) of attachment of consecutive hexagons. \
The last symbol of s is immaterial."; 

  beta::usage = "beta[s, x] is an auxiliary function giving \
H(B(n), u(n), x) where B(n) is the hexagonal chain \
consisting of n hexagons defined by the string s over {1,2,3} of length n. \
The last symbol of s is immaterial."; 

  gamma::usage = "gamma[s, x] is an auxiliary function giving \
H(B(n), v(n), x) where B(n) is the hexagonal chain \
consisting of n hexagons defined by the string s over {1,2,3} of length n. \
The last symbol of s is immaterial."; 

  Hosoya::usage = "Hosoya[s,n,x] gives the n-th Hosoya polynomial \
H(B(n),x) where B(n) is the hexagonal chain \
consisting of n hexagons defined by the string s repeated indefinitely. \
Hosoya[s,n,x,k] gives H(B(dn+k),x) where d is the \
length of s and k is between 0 and d-1.";

  HosoyaGF::usage = "HosoyaGF[s,x,z] gives the ordinary generating \
function Sum[H(B(n),x)z^n, {n,0,\[Infinity]}] of Hosoya \
polynomials H(B(n),x) where B(n) is the hexagonal chain \
consisting of n hexagons defined by the string s repeated indefinitely. \
HosoyaGF[s,x,z,k] gives the generating function Sum[H(B(dn)+\
k,x)z^n, {n,0,\[Infinity]}] where d is the length of s and \
k is between 0 and d-1.";

  WienerIndex::usage = "WienerIndex[s,n] gives the Wiener index of \
B(n), the hexagonal chain consisting of n hexagons defined by the string \
s repeated indefinitely. WienerIndex[s,n,k] gives the Wiener index of H(B\
(dn)+k,x) where d is the length of s and k is between \
0 and d-1.";

  HyperWienerIndex::usage = "WienerIndex[s,n] gives the hyper-Wiener \
index of B(n), the hexagonal chain consisting of n hexagons defined by \
the string s repeated indefinitely. HyperWienerIndex[s,n,k] gives the \
hyper-Wiener index of H(B(dn+k),x) where d is the \
length of s and k is between 0 and d-1.";

  MakePretty::usage = "MakePretty[p,n,x] tries to \
simplify the exponential polynomial p.";

  MakePrettyGF::usage = "MakePrettyGF[f,z] tries to \
simplify the rational function f.";

  MakeTeX::usage = "MakeTeX[exp] converts the \
traditional form of expression exp into TeX.";\


  Begin["`Private`"];

  alpha[{}, x_] = 2 + x;

  beta[{}, x_] = gamma[{}, x_] = 1 + x;

  alpha[s_, x_] := 
      Expand[alpha[Drop[s, -1], x] + 
          x (1 + x) (beta[Drop[s, -1], x] + 
                gamma[Drop[s, -1], x]) + 4 + 3  x + 2  x^2 + x^3];
  beta[s_, x_] := 
      Expand[Switch[Last[s], 1, 
          x beta[Drop[s, -1], x] + 1 + x + x^2 + x^3, 2, 
          x^2 beta[Drop[s, -1], x] + (1 + x)^2, 3, 
          x^2 gamma[Drop[s, -1], x] + (1 + x)^2]];
  gamma[s_, x_] := 
      Expand[Switch[Last[s], 1, 
          x^2 beta[Drop[s, -1], x] + (1 + x)^2, 2, 
          x^2 gamma[Drop[s, -1], x] + (1 + x)^2, 3, 
          x gamma[Drop[s, -1], x] + 1 + x + x^2 + 
            x^3]];
  Hadamard[f_, g_, z_Symbol] := 
      Module[{f1 = Factor[f], g1 = Factor[g], p, q, r, d, c, a, y, 
          dq}, {p, q} = Denominator /@ {f1, g1}; dq = 
          Exponent[q, z]; r = 
          Resultant[p, Expand[z^dq (q /. z -> y/z)], 
            z]; c = 
          CoefficientList[r, y]; c = -Rest[c]/
            First[c]; d = 
          Length[c]; a = 
          Times @@ (Pad[CoefficientList[Series[#, {z, 0, d - 1}], z], 
                    d] & /@ {f, g}); Return[
          ToGfn[Rec[c, a], z]]];
  ToGfn[Rec[c_List, a_List], z_Symbol] := 
      Module[{i, j, d, zz, z1, r, den, fc}, 
        d = Length[c]; zz = 
          Table[z^i, {i, d}]; z1 = 
          If[zz =!= {}, Drop[zz, -1], zz]; r = 
          Factor[((a . 
                    Table[z^i - 
                        Take[c, d - i - 1] . 
                          Take[z1, -(d - i - 1)], {i, 0, 
                        d - 1}]))/(1 - 
                  c . zz)]; Return[
          r]];
  Pad[l_List, d_Integer?(# >= 0 &)] :=  (* 
        Pad l to length d with zeros on the right \
*) Module[{k = Length[l]}, 
        If[d <= k, l, Join[l, Table[0, {d - k}]]]];
  HosoyaGF[s_, x_, z_] := 
      Module[{m = Length[s], bvars, b, n, gvars, g, vars, 
          inh1 = 1 + x + x^2 + x^3, inh2 = (1 + x)^2, beqns, i, j, 
          lowerb, lowerg, geqns, inits, eqns, sol, gfb, gfg, 
          gfa}, bvars = b[#][n] & /@ 
            Range[m]; gvars = g[#][n] & /@ 
            Range[m]; vars = Join[bvars, gvars]; 
        beqns = Table[b[i][
                n] == (lowerb = 
                  If[i > 1, b[i - 1][n], b[m][n - 1]]; 
                lowerg = If[i > 1, g[i - 1][n], g[m][n - 1]]; 
                Switch[s[[i]], 1, x lowerb + inh1, 2, 
                  x^2 lowerb + inh2, 3, x^2 lowerg + inh2]), {i, m}]; 
        geqns = Table[g[i][
                n] == (lowerb = 
                  If[i > 1, b[i - 1][n], b[m][n - 1]]; 
                lowerg = If[i > 1, g[i - 1][n], g[m][n - 1]]; 
                Switch[s[[i]], 1, x^2 lowerb + inh2, 2, 
                  x^2 lowerg + inh2, 3, x lowerg + inh1]), {i, m}]; 
        inits = Join[{b[m][0] == 1 + x, g[m][0] == 1 + x}, 
            Join @@ Table[{b[i][0] == 0, g[i][0] == 0}, {i, m - 1}]]; 
        eqns = Join[inits, beqns, geqns]; sol = 
          GeneratingFunction[eqns, vars, n, z]; sol = 
          Factor[sol[[1]]]; gfb = 
          Factor[Sum[z^(j - 1) (sol[[j]] /. z -> z^m), {j, m}]/
              z^(m - 1)]; 
        gfg = Factor[
            Sum[z^(j - m - 1) (sol[[j]] /. z -> z^m), {j, m + 1, 
                  2  m}]/z^(m - 1)]; 
        gfa = Factor[(2 + x + x (1 + x) z (gfb + gfg))/(1 - 
                    z) + z (4 + 3  x + 2  x^2 + x^3)/(1 - z)^2]; 
        Return[gfa]];
  Hosoya[s_, nn_, x_] := 
      Module[{n}, SeriesTerm[
            HosoyaGF[s, x, z], {z, 0, n}] /. n -> nn];
  HosoyaGF[s_, x_, z_, k_] /; 0 <= k <= Length[s] - 1 := 
      Module[{ls = Length[s], gf = HosoyaGF[s, x, z], num, den, q, r, qlist, 
          full}, {num, den} = {Numerator[gf], 
            Denominator[gf]}; {q, 
            r} = {PolynomialQuotient[num, den, z], 
            PolynomialRemainder[num, den, z]}; qlist = 
          CoefficientList[q, z]; 
        full =  (qlist Pad[
                    CoefficientList[
                      Series[z^k/(1 - z^ls), {z, 0, Length[qlist] - 1}], 
                      z], Length[qlist]]) . z^Range[0, Length[qlist] - 1] + 
            Hadamard[r/den, z^k/(1 - z^ls), z]; full = 
          ExpandDenominator[ExpandNumerator[Cancel[full/z^k]]] /. 
            z^m_.  -> z^(m/ls); Return[
          Factor[full]]];
  Hosoya[s_, nn_, x_, k_] /; 0 <= k <= Length[s] - 1 := 
      Module[{n}, SeriesTerm[
            HosoyaGF[s, x, z, k], {z, 0, n}] /. 
          n -> nn];
  WienerIndex[s_, nn_] := 
      Module[{n, x, 
          z}, SeriesTerm[
            Factor[D[HosoyaGF[s, x, z], x] /. x -> 1], {z, 0, n}] /. 
          n -> nn];
  WienerIndex[s_, nn_,  k_] /; 0 <= k <= Length[s] - 1 := 
      Module[{n, x, 
          z}, SeriesTerm[
            Factor[D[HosoyaGF[s, x, z, k], x] /. x -> 1], {z, 0, n}] /. 
          n -> nn];
  HyperWienerIndex[s_, nn_] := 
      Module[{n, x, z, p}, p = 
          HosoyaGF[s, x, z]; SeriesTerm[
            Factor[(D[p, x] + D[p, {x, 2}]/2) /. x -> 1], {z, 0, n}] /. 
          n -> nn];
  HyperWienerIndex[s_, nn_,  k_] /; 0 <= k <= Length[s] - 1 := 
      Module[{n, x, z, p}, p = 
          HosoyaGF[s, x, z, k]; SeriesTerm[
            Factor[(D[p, x] + D[p, {x, 2}]/2) /. x -> 1], {z, 0, n}] /. 
          n -> nn];
  MakePretty[p_, n_, x_] := 
      Module[{pp = Factor[p], num, den, nfree, 
          nfull}, num = 
          Expand[Numerator[pp]]; den = 
          Factor[Expand[Denominator[pp]]]; nfree = 
          Select[num, FreeQ[#, n] &]; nfull = 
          Collect[num - nfree, {n, 
              x^n, (-1)^n}]; num = 
          Factor[nfree] + Factor /@ nfull; Return[
          TraditionalForm[num/den]]];
  MakePrettyGF[p_, z_] := 
      Module[{num = Numerator[p], den = Denominator[p]}, 
        Return[TraditionalForm[
            Collect[num, z]/den]]];
  MakeTeX[x_] := TeXForm[TraditionalForm[x]];

  On[General::spell1, 
      General::spell];

  End[];

  EndPackage[];