(*                                                                 *)
(* fie-w   Mathematica program for computing Fiedler's             *)
(*         (long) knot invariant W                                 *)
(*                                                                 *)
(* Author: S.Orevkov   orevkov(at)math.ups-tlse.fr                 *)
(*                                                                 *)
(* Usage:  << fie-w                                                *)
(*         W[knot]                                                 *)
(*                                                                 *)
(* where:  fie-w       is this file                                *)
(*         knot        is a list {Wr,reg}                          *)
(*         Wr          is list of writhes of crossings (+1 or -1)  *)
(*         reg         is the list of region-vectors of crossings  *)
(*                                                                 *)
(* Known bug:  The program W does not work correctly when the knot *)
(*             diagram is non-minimal in the sense that there is a *)
(*             component of its complement whose closure is not    *)
(*             simply connected (for example, diagrams obtained by *)
(*             the first Reidemeister move).                       *)
(*                                                                 *)
(* We suppose that all crossings are indexed by  1,2,...,n  and    *)
(* the regions of the complement of the knot diagram are indexed   *)
(* indexed by  1,2,...,m  so that  1  and  2  correspond to the  *)
(* unbounded regions (recall that the knot is long, hence there    *)
(* are 2 unbounded regions).                                       *)
(*                                                                 *)
(* A region-vector of a crossing is the vector {j_1,j_2,j_3,j_4}   *)
(* where  j_k  is the number of the region corresponding to the    *)
(* k-th quarter at the crossing. The quarters are numerated in the *)
(* standard way: 1,2,3,4  counterclockwise starting with the       *)
(* quarter between the to arrows:                                  *)
(*      A                                                          *)
(*   2  |  1                                                       *)
(*  ----+---->                                                     *)
(*   3  |  4                                                       *)
(*      |                                                          *)
(*                                                                 *)
(* Examples:                                                       *)
(*           Wr3={1,1,1}; reg3={{4,3,1,2},{5,3,4,2},{1,3,5,2}};    *)
(*           trefoil={Wr3,reg3};                                   *)
(*           W[trefoil]                                            *)
(*                                                                 *)
(*           Wr8 = {1, 1, -1, -1};                                 *)
(*           reg8 = {{3,1,2,4},{2,5,3,4},{1,6,5,2},{5,6,1,3}};     *)
(*           figure8 = {Wr8,reg8};                                 *)
(*                                                                 *)
(*           W[figure8]                                            *)
(*                                                                 *)
(*  See figures in  fie-w.ps                                       *)
(*                                                                 *)
W = Function[knot,
      Module[{Wr, reg, n, w, contr, TabAS, AS, 
        rr, nr, j0, r, i, j, P, JS, col, j1, j2, c1, c2, f, s},
        Wr = knot[[1]]; reg = knot[[2]];
        n = Length[Wr];
        w = (Wr + 3)/2;
        
(* contr[[(w + 1)/2, js + 1, as]]  is the contribution of a crossing *)
(* if Alexander state is 'as', Jones state is 'js', and writhe is  w *)
        
        contr = {
            {{x^2, -I*x*y, y^2, -I*x*y}, {I*x/y, 1, I*y/x, 1}},
            {{1/x^2, I/x/y, 1/y^2, I/x/y}, {-I/x*y, 1, -I*x/y, 1}}};
        
        TabAS = {};
        AS = Table[1, {n}]; rr = {}; nr = 0;
        While[True,
          (*Print["========================="];
                            Print[AS, " nr=", nr, " rr=", rr];*)
          Do[j0 = 0;
          Do[If[(r = reg[[i, j]]) < 3, Continue[]];
            If[MemberQ[rr, r], Continue[]];
            j0 = j; Break[], {j, AS[[i]], 4}];
          (* Print["  i=", i, " -> j0=", j0]; *)
            If[j0 > 0, AS[[i]] = j0; AppendTo[rr, r]; nr++, Break[]],
            {i, nr + 1, n}
          ];
          If[nr == n, (*Print[AS];*)  AppendTo[TabAS, AS];
            AS[[n - 1]]++; AS[[n]] = 1; 
            rr = Delete[rr, {{n - 1}, {n}}]; nr = n - 2,
 
            If[nr == 0, Break[]];
            AS[[nr]]++; Do[AS[[i]] = 1, {i, nr + 1, n}];
            rr = Delete[rr, nr]; nr--
          ]
        ];
        P = 0;
        JS = Table[0, {n}];
        While[True,
          
          (* Compute the number of circles *)
          col = Table[i, {i, n + 2}];
          Do[
            If[JS[[i]] == 1, j1 = 1; j2 = 3, j1 = 2; j2 = 4];
            c1 = col[[reg[[i, j1]]]]; c2 = col[[reg[[i, j2]]]];
            Do[If[col[[j]] == c2, col[[j]] = c1], {j, n + 2}],
            {i, n}
          ];
          f = (I*(x*y - 1/x/y))^(Length[Union[col]] - 2);
          s = 0;
          Do[
            s += Product[contr[[w[[j]],JS[[j]]+1, TabAS[[k,j]]]], {j,n}],
            {k, Length[TabAS]}
          ];
          P += s*f;
          
          (*pass to the next JS*)
          i0 = 0;
          Do[If[JS[[i]] == 0, i0 = i; Break[]], {i, n, 1, -1}];
          If[i0 == 0, Break[]];
          JS[[i0]] = 1;
          Do[JS[[i]] = 0, {i, i0 + 1, n}]
        ];
        Factor[P]
]];