BeginPackage["SzegedIndex`", "DiscreteMath`Combinatorica`"];

SzegedIndex::usage = "Let G be a graph and X a binary relation on the
vertex set of G. Let F denote the fasciagraph composed of n copies of G
where adjacencies between vertices of two neighbouring copies of G are
determined by X. Then SzegedIndex[G, X, n] gives the Szeged index of F
for large n, provided that each copy of G is an isometric subgraph of F."

Unprotect[Ceiling];
  Ceiling[a_?IntegerQ] := a;
  Ceiling[a_?IntegerQ + b_] := a + Ceiling[b];
  Ceiling[Literal[Times[a__?IntegerQ] + b_.]] := Times[a] + Ceiling[b];
Protect[Ceiling]; 

Unprotect[Floor];
  Floor[a_?IntegerQ] := a;
  Floor[a_?IntegerQ + b_] := a + Floor[b];
  Floor[Literal[Times[a__?IntegerQ] + b_.]] := Times[a] + Floor[b];
Protect[Floor]; 

Begin["`Private`"];

TropicalProduct[a_, b_] := Inner[Plus, a, b, Min];
 
ConstantMatrixQ[a_] := a - a[[1,1]] === 0 * a;

PolySum[t_, {i_, a_Integer, b_Integer}] :=
  Plus @@ Table[t, {i, a, b}];
PolySum[t_, {i_, a_, b_}] /; a =!= 0 := 
  PolySum[Expand[t /. i :> i+a], {i, 0, b-a}];
PolySum[t__Plus, it_] := Expand[PolySum[#, it]& /@ t];
PolySum[c_ t_, {i_, a_, b_}] /; FreeQ[c, i] := 
  Expand[c PolySum[t, {i, a, b}]];
PolySum[c_, {i_, a_, b_}] /; FreeQ[c, i] := Expand[c (b - a + 1)];
PolySum[i_, {i_, 0, n_}] := Expand[n(n + 1)/2];
PolySum[i_^2, {i_, 0, n_}] := Expand[n(n + 1)(2n + 1)/6];
PolySum[i_^m_, {i_, 0, n_}] := 
  Expand[(BernoulliB[m + 1, n + 1] - BernoulliB[m + 1])/(m + 1)];  

valPlus[p_, d_, displ_, q_, da_] := 
  Module[{sp, sd, m, mat, mat2, nuP = {}, nvP = {}}, 
    Do[mat = da[[m]]; mat2 = da[[m + displ]]; 
       sp = Count[mat2[[p]] - mat[[d]], _?Negative]; 
       sd = Count[mat2[[p]] - mat[[d]], _?Positive]; 
       AppendTo[nuP, sp]; AppendTo[nvP, sd], {m, q}]; 
    Return[{nuP, nvP}]];
 
valMinus[p_, d_, displ_, q_, da_] := 
  Module[{sp, sd, m, mat, mat2, nuM = {}, nvM = {}}, 
    Do[mat = Transpose[da[[m]]]; mat2 = Transpose[da[[m + displ]]]; 
       sp = Count[mat[[p]] - mat2[[d]], _?Negative]; 
       sd = Count[mat[[p]] - mat2[[d]], _?Positive]; 
       AppendTo[nuM, sp]; AppendTo[nvM, sd], {m, q}]; 
    Return[{nuM, nvM}]];
 
sze[displ_, p_, q_, pp_, da_, m_][e_, remainder_] := 
  Module[{i, j, ee, nuP, nvP, nuM, nvM, a, n = pp * m + remainder}, 
    {nuP, nvP} = valPlus[e[[1]], e[[2]], displ, q, da]; 
    {nuM, nvM} = valMinus[e[[1]], e[[2]], displ, q, da]; 
    i/: IntegerQ[i] = True; 
    j/: IntegerQ[j] = True; 
    a = Evaluate //@ 
      (PolySum[	(PolySum[Ceiling[(j-q-i+1)/pp] * nuM[[p+i]], {i, 1, pp}] 
        + PolySum[nuM[[i+1]], {i, 0, q-1}] 
        + PolySum[nuP[[i+1]], {i, 1, q-1}] 
        + PolySum[Ceiling[(n-q-j-i+2)/pp] * nuP[[p+i]], {i, 1, pp}]) 
      * (PolySum[Ceiling[(j-q-i+1)/pp] * nvM[[p+i]], {i, 1, pp}] 
        + PolySum[nvM[[i+1]], {i, 0, q-1}] 
        + PolySum[nvP[[i+1]], {i, 1, q-1}] 
        + PolySum[Ceiling[(n-q-j-i+2)/pp] * nvP[[p+i]], {i, 1, pp}]), 
           {j, q, n-q+1-displ}] 
     + PolySum[ (PolySum[nuM[[i]], {i, 1, j}] 
        + PolySum[nuP[[i+1]], {i, 1, q-1}] 
        + PolySum[Ceiling[(n-q-j-i+2)/pp] * nuP[[p+i]], {i, 1, pp}])
      * (PolySum[nvM[[i]], {i, 1, j}] 
        + PolySum[nvP[[i+1]], {i, 1, q-1}] 
        + PolySum[Ceiling[(n-q-j-i+2)/pp] * nvP[[p+i]], {i, 1, pp}]),
           {j, 1, q-1}] 
     + PolySum[ (PolySum[Ceiling[(j-q-i+1)/pp] * nuM[[p+i]], {i, 1, pp}] 
        + PolySum[nuM[[i+1]], {i, 0, q-1}] 
        + PolySum[nuP[[i+1]], {i, 1, n-j}])
      * (PolySum[Ceiling[(j-q-i+1)/pp] * nvM[[p+i]], {i, 1, pp}] 
        + PolySum[nvM[[i+1]], {i, 0, q-1}] 
        + PolySum[nvP[[i+1]], {i, 1, n-j}]), 
           {j, n-q+2-displ, n-displ}] 
      /. Ceiling[x_] :> Ceiling[Expand[x]]); 
    If[pp > 1,  
      ee /: IntegerQ[ee] = True; 
      a = a /. Literal[PolySum[x_, {ind_, u_, v_}]] :> 
            Plus @@ Table[PolySum[x /. ind :> pp * ee + Mod[r, pp], 
              {ee, Ceiling[(u-r)/pp], Floor[(v-r)/pp]}], {r, 0, pp-1}]
        /. {Ceiling[x_] :> Ceiling[Expand[x]], 
            Floor[x_] :> Floor[Expand[x]]}; 
      a = a /. Literal[PolySum[x_, {ind_, u_, v_}]] :> 
            PolySum[Expand[x], {ind, u, v}] ];
    Return[a] ];
 
SzegedIndex[g_Graph, x_List, n_Symbol] := 
  Module[{continue, d, v, t, edges, di, d1, i, j, da, p, q, 
          pp, ee, remainder, xx, result, printout, nn}, 
    Off[Part::pspec]; 
    nn/: IntegerQ[nn] = True; 
    d = AllPairsShortestPath[g]; v = Length[Vertices[g]]; 
    t = ReplacePart[Table[Infinity, {i, v}, {j, v}], 1, x]; 
    edges = ToUnorderedPairs[g]; di = TropicalProduct[d, t]; 
    d1 = TropicalProduct[di, d]; di = d1; i = 1; da = {d1}; 
    continue = True; 
    While[continue, 
      di = TropicalProduct[di, d1]; AppendTo[da, di]; 
      While[i > 0 && continue, 
        If[ConstantMatrixQ[di - da[[i]]],
           continue = False; p = i; q = Length[da],
             i--] ]; 
      If[continue, i = Length[da]] ]; 
    da = Prepend[da, d]; pp = q - p; 
    Do[ee = {edges, Table[remainder, {Length[edges]}]}; 
       xx = {x, Table[remainder, {Length[x]}]}; 
       result = Factor[
         Plus @@ MapThread[sze[0, p, q, pp, da, nn], ee] + 
         Plus @@ MapThread[sze[1, p, q, pp, da, nn], xx]] /. nn -> n; 
       If[pp > 1, 
          printout = {n >= Ceiling[(2p - 1 - remainder) / pp], 
                      ":   Sz(", pp * n}; 
          If[remainder > 0, 
             printout = Join[printout, {"+", remainder}] ]; 
          printout = Join[printout, {") = ", result}], 
            printout = {n >= 2p - 1, ":   Sz(n) = ", result}]; 
       Print @@ printout, 
    {remainder, 0, pp - 1}];
    On[Part::pspec] ];

End[ ];

EndPackage[ ];