(** qp3.m : quasipositivity problem for 3-braids      **)
(** qp3[b] = True/False (True iff b is quasipositive) **)
(**                                                   **)
(** Author; S.Orevkov                                 **)  
(**            http://picard.ups-tlse.fr/~orevkov     **)
(** Last change: January 12, 2011                     **)
(**                                                   **)
(** Input: b list of generators of the braid group B3 **)
(**  1 and  2 stand for sigma_1 and sigma_2           **)
(** -i (i=1,2) stands for the inverse of sigma_i      **)
(**                                                   **)
qp3 = Function[b, Module[{deb=False,pw=0,a={},rc=False,
  i,bi,b0,aa,a0,a1,qp,br,e,e0,A,qptest,cycred},
  Do[bi=b[[i]];
    If[ a=={},
      If[ bi>0,
        a={1}; b0=bi,
      (*else*)
        pw-=1; a={1,1}; b0=bi+3
      ],
    (*else*)
      If[ bi>0,
        If[ bi==b0,
	  a[[-1]]++,
	(*else*)
	  If[ Length[a]==1,
	    AppendTo[a,1]; b0=bi,
	  (*else*)
	    If[ a[[-1]]>1,
	      AppendTo[a,1]; b0=bi,
	    (*else*)
	      a=Delete[a,-1]; a[[-1]]--; pw++;
	      If[a[[-1]]==0, a=Delete[a,-1]; b0=bi ]
	    ]
	  ]
	],
      (*else*)
        If[ bi==-b0,
	  a[[-1]]--;
	  If[ a[[-1]]==0, a=Delete[a,-1]; b0=3+bi ],
	(*else*)
	  pw--; a[[-1]]++; AppendTo[a,1]
	]
      ]
    ],
  {i,Length[b]}];
  If[ Length[a]>2 && Mod[Length[a]+pw,2]==1,
    a[[1]]+=a[[-1]]; a=Delete[a,-1]
  ];
(*                                                           *)
(* Now the braid is in Garside normal form                   *)
(* Delta^pw sigma_1^a[[1]] sigma_2^a[[2]] sigma_1^a[[3]] etc *)
(* where each a[[i]] is geater than 1                        *)
(*                                                           *)
  If[deb,Print[{pw,a}]];
  e=3pw+Plus@@a;
  cycred=Function[Null,   (* cyclic reduction *)
    While[Length[a]>1,
      If[ a[[1]]==1 && a[[-1]]==1,
        If[ Length[a]==2, Break[] ];
        pw++;
        If[ Length[a]==3, a={a[[2]]-1}; Break[] ];
        a[[2]] += (a[[-2]]-1);
	a=Delete[Delete[Delete[a,1],-1],-1];
        Break[]
      ];
      If[ a[[-1]]==1, a=Prepend[Delete[a,-1],1] ];
      If[ a[[1]]==1,
        pw++; a=Delete[a,1]; a[[1]]--; a[[-1]]--,
      (*else*)
        Break[]
      ]
    ];
    aa=a;           (* minimize under cycling and symetries *)
    Do[aa=Reverse[aa];
      Do[ a1=aa[[-1]]; aa=Prepend[Delete[aa,-1],a1];
        a=Sort[{a,aa}][[1]],
      {Length[a]-1}],
    {2}]
  ];
  If[deb,Print[{pw,a}]];
  cycred[];
  If[deb,Print[{pw,a}]];
  qptest=Function[Null, (* Return True if the problem is solved *)
    If[ pw>=0,      (* In this case qp is True if quasipositive *)
      qp=True; True,
    (*else*)
      If[  e <  2,
        qp=False; True,
      (*else*)
        If[ pw+Length[a] < 2,
          qp=False; True,
        (*else*)
          If[ pw >= -3, qp=True; True, False ]
        ]
      ]
    ]
  ];
  rc=br=False;
  If[ qptest[], rc=qp; br=True ];
  A[e]={{pw,a}};
  Do[A[i]={},{i,2,e-1}]; e0=e;
  While[ e0>1 && Not[br],  If[deb,Print["e0=",e0," ",A[e0]]];
    Do[a=A[e0][[iA]]; pw0=a[[1]]; a0=a[[2]];
      If[deb,Print["  ",a]];
      Do[a1=a0[[-1]]; a=a0=Prepend[Delete[a0,-1],a1]; pw=pw0;
        If[deb,Print["    e=",e0," ",{pw,a}]];
        e=e0-a[[1]]+1;
	If[e<0,Continue[]];
        a[[1]]=1; pw=pw0; cycred[];
	If[deb,Print["    e=",e, " ",{pw,a}]];
	If[ qptest[],
	  If[ qp,
	    br=rc=True; Break[],
	  (*else*)
	    Continue[]
	  ],
	(*else*)
	  A[e] = Union[A[e],{{pw,a}}]
	],
      {Length[a0]}];
      If[br,Break[]],
    {iA,Length[A[e0]]}];
    e0--    
  ];
  rc
]];
