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[ ];