#Implemented by Yongjae Cha and M.van Hoeij, Florida State Univ. at 2008.10#
#Implemented in Maple code#

sing:=proc(L,ext) 
        local F, i;
        F := lcoeff(L, _Env_LRE_tau)*tcoeff(L, _Env_LRE_tau); 
        F := {seq(i[1], i = `LREtools/g_p/factors`(F, _Env_LRE_x, ext))};       
        F := {seq(`LREtools/g_p/shorten_int`(i[1], _Env_LRE_x), i = `LREtools/g_p/int_dif`(F, _Env_LRE_x))}; 
        PolynomialTools:-Sort([op(F)], _Env_LRE_x)      
end proc;

plist:=proc(sing,L)
        local res,s,minmax;
        res := {};
        for s in sing do
                minmax:=`LREtools/g_p`(primpart(L,_Env_LRE_tau), RootOf(s,_Env_LRE_x))[1];
                minmax := {min(op(minmax)), max(op(minmax))};
                if minmax<>{0} then
                res := res union { [s, minmax] }
                fi;
        od;
        res
end;

combinations:=proc(pl)  
        local s,t1,t2,i,t;
        if pl = {} then return {1} fi;
        s := pl[1];
        t := procname(pl minus {s});
        t1 := s[1]^min(op(s[2])) * eval(s[1],_Env_LRE_x=_Env_LRE_x - 1)^max(op(s[2]));  
        t2 := s[1]^max(op(s[2])) * eval(s[1],_Env_LRE_x=_Env_LRE_x - 1)^min(op(s[2]));
        {seq(t1*i, i=t), seq(t2*i, i=t)}
end;

sqsols := proc (L) 
        local M, v, u;  
        M := y(_Env_LRE_x+1)+subs(_Env_LRE_x = 2*_Env_LRE_x, tcoeff(L, _Env_LRE_tau))*y(_Env_LRE_x); 
        v := LREtools[hypergeomsols](M, y(_Env_LRE_x), {},output=basis); 
        u := subs(_Env_LRE_x = (1/2)*_Env_LRE_x, v)
end proc;

checkgauge:=proc(comb,L)
        local pL,c,const;
        const:=(lcoeff(coeff(L,_Env_LRE_tau,0),_Env_LRE_x))/(lcoeff(coeff(L,_Env_LRE_tau,2),_Env_LRE_x));
        for c in comb do
          if findgauge(_Env_LRE_tau^2+const*c,L)<>FAIL then
                return [_Env_LRE_tau^2+const*c,findgauge(_Env_LRE_tau^2+const*c,L)]
          fi
        od;     
        FAIL
end;

check_3_curvature := proc(L1)
        local a,b,c, pr, pc_trace,L;
        a,b,c := coeff(L,_Env_LRE_tau,2), coeff(L,_Env_LRE_tau,1), coeff(L,_Env_LRE_tau,0);
        pr := Expand(a * subs(_Env_LRE_x=_Env_LRE_x+1, c) * subs(_Env_LRE_x=_Env_LRE_x+2, b)) mod 3;
        pc_trace := pr + subs(_Env_LRE_x=_Env_LRE_x+1,pr) + subs(_Env_LRE_x=_Env_LRE_x+2, pr) - b * subs(_Env_LRE_x=_Env_LRE_x+1,b) * subs(_Env_LRE_x=_Env_LRE_x+2,b);
        evalb( 0 = Expand(pc_trace) mod 3 )
end:

findgauge:=proc(L1,L2)
        local G,eq0,eq1,eq2,eq3,eq4,c0_x1,eq5,eq6,c0_x0,eqc1,c,d,K,cfs,e;
        cfs := indets([L1,L2],function);
        cfs := cfs union subs(_Env_LRE_x=_Env_LRE_x +1, cfs);
        cfs := cfs union subs(_Env_LRE_x =_Env_LRE_x +2, cfs);
        G := c1(_Env_LRE_x )*_Env_LRE_tau + c0(_Env_LRE_x );
        K:=`LREtools/rightdivision`( `LREtools/mult`(L2, G), L1)[2];
        eq0, eq1 := coeff(%,_Env_LRE_tau,0), coeff(%,_Env_LRE_tau,1);
        if coeff(L1,_Env_LRE_tau,1)=0 and coeff(L2,_Env_LRE_tau,1)<>0 then
                c:=solve(eq1, c0(_Env_LRE_x +1));
                eval(eq0, {c0(_Env_LRE_x ) = subs(_Env_LRE_x  = _Env_LRE_x -1, %), c0(_Env_LRE_x +2) = subs(_Env_LRE_x  = _Env_LRE_x +1, %)});
                collect(%, {c1(_Env_LRE_x ),c1(_Env_LRE_x +1),c1(_Env_LRE_x +2),c1(_Env_LRE_x +3),c1(_Env_LRE_x +4)}, factor);
                d:=LREtools[ratpolysols](%, c1(_Env_LRE_x ), {}, output = basis);
                if whattype(d)=exprseq then return FAIL fi;
                d:=d[1];
                eval(c, {c1(_Env_LRE_x ) = d, c1(_Env_LRE_x +2) = subs(_Env_LRE_x =_Env_LRE_x +2,d), c1(_Env_LRE_x +1) = (_Env_LRE_x =_Env_LRE_x +1,d)});
                return collect(simplify(d*_Env_LRE_tau+subs(_Env_LRE_x =_Env_LRE_x -1,%)),_Env_LRE_tau)
        elif coeff(L1,_Env_LRE_tau,1)<>0 and coeff(L2,_Env_LRE_tau,1)=0 then
                c:=solve(eq1, c0(_Env_LRE_x +2));
                eval(eq0, {c0(_Env_LRE_x +2) = %, c0(_Env_LRE_x ) = eval(%, _Env_LRE_x  = _Env_LRE_x -2)});
                d:=LREtools[ratpolysols](%, c1(_Env_LRE_x ), {}, output = basis);
                if whattype(d)=exprseq then return FAIL fi;
                d:=d[1];
                eval(c, {c1(_Env_LRE_x ) = %, c1(_Env_LRE_x +2) = subs(_Env_LRE_x =_Env_LRE_x +2,%), c1(_Env_LRE_x +1) = (_Env_LRE_x =_Env_LRE_x +1,%)});
                return collect(simplify(d*_Env_LRE_tau+subs(_Env_LRE_x =_Env_LRE_x -2,%)),_Env_LRE_tau)
        elif coeff(L1,_Env_LRE_tau,1)=0 and coeff(L2,_Env_LRE_tau,1)=0 then
                d:=LREtools[ratpolysols](eq1, c1(_Env_LRE_x ), {}, output = basis);
                e:=LREtools[ratpolysols](eq0, c0(_Env_LRE_x ), {}, output = basis);
                if whattype(d)=exprseq then d:=0
                else d:=d[1]
                fi;
                if whattype(e)=exprseq then e:=0
                else e:=simplify(e[1])
                fi;
                if d*tau+e=0 then return FAIL fi;
                return d*tau+e
        else
                eq2 := coeff(eq1, c0(_Env_LRE_x +2) ) * eq0 - coeff(eq0, c0(_Env_LRE_x +2) ) * eq1; # eq2 has no c0(_Env_LRE_x +2) term.
                eq3 := subs(_Env_LRE_x =_Env_LRE_x +1, eq2);
                eq4 := coeff(eq1, c0(_Env_LRE_x +2) ) * eq3 - coeff(eq3, c0(_Env_LRE_x +2) ) * eq1; # Has only c0(_Env_LRE_x +1) term
                c0_x1 := solve( expand(eq4), c0(_Env_LRE_x +1) );
                eq5 := coeff(eq2, c0(_Env_LRE_x +1) ) * eq4 - coeff(eq4, c0(_Env_LRE_x +1) ) * eq2;
                eq6 := coeff(eq0, c0(_Env_LRE_x +2) ) * eq3 - coeff(eq3, c0(_Env_LRE_x +2) ) * eq0;
                if not has(eq2, c0(_Env_LRE_x )) then
                        c0_x0 := solve(expand(eq6), c0(_Env_LRE_x ))
                else
                        c0_x0 := solve(expand(eq5), c0(_Env_LRE_x ) )
                fi;
                eqc1 := subs(_Env_LRE_x =_Env_LRE_x +1, c0_x0) - c0_x1;
                # vars := indets(eqc1,function) minus cfs;
                # collect(primpart(eqc1,vars), vars, factor);
                d:=LREtools[ratpolysols](%, c1(_Env_LRE_x ), {}, output = basis);
                d:=d[1];
                if whattype(d)=exprseq then return FAIL
                elif nops(d)>1 then
                        error "input must be reducible"
                fi;
                d := d[1];
                c := subs(c1(_Env_LRE_x )=d, c1(_Env_LRE_x +1)=subs(_Env_LRE_x =_Env_LRE_x +1,d), c2(_Env_LRE_x +2)=subs(_Env_LRE_x =_Env_LRE_x +2,d), c0_x0);
                c:=normal(c)*tau + d;
                collect(simplify(c),tau)
        fi;
end:

tausqsols:=proc(LL)
        local s,c,csol,sqsol,t1,t2,L,ext,pl,pl2,F,i;
        L := evala(Primpart(LL,_Env_LRE_tau));
        if not check_3_curvature(L) then
                return FAIL
        fi;
        ext := indets([args], {'RootOf', 'radical', 'nonreal'}); 
        s:=sing(L,ext);
        pl:=plist(s,L);
        c:=combinations(pl);
        csol:=checkgauge(c,L/lcoeff(L,_Env_LRE_tau));
        if csol <> FAIL then
	        sqsol:=op(sqsols(csol[1]));
	        t1:=coeff(csol[2],_Env_LRE_tau,0)*sqsol;
	        t2:=coeff(csol[2],_Env_LRE_tau,1)*subs(_Env_LRE_x = _Env_LRE_x+1, sqsol);
	        return (C[1]+((-1)^_Env_LRE_x)*C[2])*t1+(C[1]-((-1)^_Env_LRE_x)*C[2])*t2
	fi;
	# We're done searching over the base field, now search only over
	# algebraic extensions
	pl2 := [seq(`if`(nops(i[2])>1, i[1], NULL), i=pl)];
	if pl2=[] then
		return FAIL
	fi;
	for i in pl2 do if type(degree(i, _Env_LRE_x), odd) then
		return FAIL
	fi od;
	pl2 := PolynomialTools:-Sort(pl2, _Env_LRE_x);
	for F in `LREtools/hsolsR/subfields`(2, [pl2[1]], _Env_LRE_x, ext) do
         c:=combinationsF(pl,F);
         csol:=checkgauge(c,L/lcoeff(L,_Env_LRE_tau));
         if csol <> FAIL then
                sqsol:=op(sqsols(csol[1]));
                t1:=coeff(csol[2],_Env_LRE_tau,0)*sqsol;
                t2:=coeff(csol[2],_Env_LRE_tau,1)*subs(_Env_LRE_x = _Env_LRE_x+1, sqsol);
                return (C[1]+((-1)^_Env_LRE_x)*C[2])*t1+(C[1]-((-1)^_Env_LRE_x)*C[2])*t2
         fi;
	od;
	FAIL
end;

combinationsF:=proc(pl,F)
        local s,t1,t2,i,t,v;
        if pl = {} then return {1} fi;
        s := pl[1];
        t := procname(pl minus {s},F);
	if t={} then
		{}
	elif nops(s[2])=1 then
		t1 := (s[1]*eval(s[1],_Env_LRE_x=_Env_LRE_x - 1))^s[2][1];
		{seq(t1*i, i=t)}
	else
		v := map(i -> i[1],`DEtools/kovacicsols/factors`(s[1],_Env_LRE_x,F));
		if nops(v)<>2 then return {} fi;
		t1 := (v[1]*eval(v[2],_Env_LRE_x=_Env_LRE_x - 1))^min(op(s[2]))
		*     (v[2]*eval(v[1],_Env_LRE_x=_Env_LRE_x - 1))^max(op(s[2]));
                t2 := (v[1]*eval(v[2],_Env_LRE_x=_Env_LRE_x - 1))^max(op(s[2]))
                *     (v[2]*eval(v[1],_Env_LRE_x=_Env_LRE_x - 1))^min(op(s[2]));
		{seq(t1*i, i=t), seq(t2*i, i=t)}
	fi
end;