//////////////
// SETUP:
//////////////
_<l> := PolynomialRing(Rationals());
K<l> := NumberField(l^2-3);
P<x0,x1,x2,x3> := PolynomialRing(K,4);
f := x0^2+x1^2+x2^2-25*x3^2;
g := (x0+l*x3)*(x0-l*x1-3*x3)*(x0+l*x1-3*x3)-2*x3^3;

CC := ComplexField(250);
_<l> := PolynomialRing(CC);
ll := Roots(l^2-3)[2][1];
roundingMap1 := hom<K->CC|Evaluate(l,ll)>;

CCout:= ComplexField(10);
roundingMap2 := hom<CC->CCout|>;



//////////////
// Computation
//////////////
print "computing multiple root loci";
S<a0,a1,a2,a3,a4,a5,a6> := PolynomialRing(Rationals(),7);
R<b0,b1,b2,b3,t> := PolynomialRing(Rationals(),5);
phi := hom<S->R|Coefficients((b0+b1*t+b2*t^2+b3*t^3)^2,t)>;
ker := MinimalBasis(PolyMapKernel(phi));


print "forming intersection of curve with parametric plane ax+by+cz+d=0";
Ru<u0,u1,u2,u3> := PolynomialRing(K,4);
Ku<u0,u1,u2,u3> := FieldOfFractions(Ru);
Pu<x0,x1,x2,x3> := PolynomialRing(Ku,4);
plane := u0*x0+u1*x1+u2*x2+u3*x3;
iota := hom<P->Pu|[x0,x1,x2,x3]>;
f := iota(f);
g := iota(g);


print "projecting onto all axes and forcing projection to be 3 double points";
f3 := Evaluate(f,[x0,x1,x2,(-u0*x0-u1*x1-u2*x2)/u3]);
g3 := Evaluate(g,[x0,x1,x2,(-u0*x0-u1*x1-u2*x2)/u3]);
r32,den := ClearDenominators(Resultant(f3,g3,x2));
r32 := r32*den;
C32 := Coefficients(Evaluate(r32,[x0,1,0,0]));
r31,den := ClearDenominators(Resultant(f3,g3,x1));
r31 := r31*den;
C31 := Coefficients(Evaluate(r31,[x0,0,1,0]));
r30,den := ClearDenominators(Resultant(f3,g3,x0));
r30 := r30*den;
C30 := Coefficients(Evaluate(r30,[0,x1,1,0]));

f2 := Evaluate(f,[x0,x1,(-u0*x0-u1*x1-u3*x3)/u2,x3]);
g2 := Evaluate(g,[x0,x1,(-u0*x0-u1*x1-u3*x3)/u2,x3]);
r23,den := ClearDenominators(Resultant(f2,g2,x3));
r23 := r23*den;
C23 := Coefficients(Evaluate(r23,[x0,1,0,0]));
r21,den := ClearDenominators(Resultant(f2,g2,x1));
r21 := r21*den;
C21 := Coefficients(Evaluate(r21,[x0,0,0,1]));
r20,den := ClearDenominators(Resultant(f2,g2,x0));
r20 := r20*den;
C20 := Coefficients(Evaluate(r20,[0,x1,0,1]));

f1 := Evaluate(f,[x0,(-u0*x0-u2*x2-u3*x3)/u1,x2,x3]);
g1 := Evaluate(g,[x0,(-u0*x0-u2*x2-u3*x3)/u1,x2,x3]);
r13,den := ClearDenominators(Resultant(f1,g1,x3));
r13 := r13*den;
C13 := Coefficients(Evaluate(r13,[x0,0,1,0]));
r12,den := ClearDenominators(Resultant(f1,g1,x2));
r12 := r12*den;
C12 := Coefficients(Evaluate(r12,[x0,0,0,1]));
r10,den := ClearDenominators(Resultant(f1,g1,x0));
r10 := r10*den;
C10 := Coefficients(Evaluate(r10,[0,0,x2,1]));

f0 := Evaluate(f,[(-u1*x1-u2*x2-u3*x3)/u0,x1,x2,x3]);
g0 := Evaluate(g,[(-u1*x1-u2*x2-u3*x3)/u0,x1,x2,x3]);
r03,den := ClearDenominators(Resultant(f0,g0,x3));
r03 := r03*den;
C03 := Coefficients(Evaluate(r03,[0,x1,1,0]));
r02,den := ClearDenominators(Resultant(f0,g0,x2));
r02 := r02*den;
C02 := Coefficients(Evaluate(r02,[0,x1,0,1]));
r01,den := ClearDenominators(Resultant(f0,g0,x1));
r01 := r01*den;
C01 := Coefficients(Evaluate(r01,[0,0,x2,1]));

iota := hom<Pu->Ru|[0,0,0,0]>;
C32 := iota(C32);
C31 := iota(C31);
C30 := iota(C30);
psi32 := hom<S->Ru|C32>;
C32 := psi32(ker);
psi31 := hom<S->Ru|C31>;
C31 := psi31(ker);
psi30 := hom<S->Ru|C30>;
C30 := psi30(ker);

C23 := iota(C23);
C21 := iota(C21);
C20 := iota(C20);
psi23 := hom<S->Ru|C23>;
C23 := psi23(ker);
psi21 := hom<S->Ru|C21>;
C21 := psi21(ker);
psi20 := hom<S->Ru|C20>;
C20 := psi20(ker);

C13 := iota(C13);
C12 := iota(C12);
C10 := iota(C10);
psi13 := hom<S->Ru|C13>;
C13 := psi13(ker);
psi12 := hom<S->Ru|C12>;
C12 := psi12(ker);
psi10 := hom<S->Ru|C10>;
C10 := psi10(ker);

C03 := iota(C03);
C02 := iota(C02);
C01 := iota(C01);
psi03 := hom<S->Ru|C03>;
C03 := psi03(ker);
psi02 := hom<S->Ru|C02>;
C02 := psi02(ker);
psi01 := hom<S->Ru|C01>;
C01 := psi01(ker);


print "gathering equations, restricting to affine chart d=1";
C := C32 cat C31 cat C30 cat C23 cat C21 cat C20 cat C13 cat C12 cat C10 cat C03 cat C02 cat C01;
Su<u0,u1,u2> := PolynomialRing(K,3);
pi := hom<Ru->Su|[u0,u1,u2,1]>;
I := ideal<Su|pi(C)>;
SetVerbose("Groebner",1);
// VarietySizeOverAlgebraicClosure(I); // omit this step if too costly


print "computing radical decomposition";
RD := RadicalDecomposition(I);


print "computing tritangents numerically";
CCu<u0,u1,u2> := PolynomialRing(CC,3);
phi := hom<Su->CCu|roundingMap1,u0,u1,u2>;
RRD := [phi(Generators(rd)) : rd in RD];

realRoots := [];
complexRoots := [];
for H in RRD do
  _,h := IsUnivariate(H[3]);
  R2 := Roots(h);
  for rr2 in R2 do
    r2 := rr2[1];
    _, h := IsUnivariate(Evaluate(H[2],[u0,u1,r2]));
    R1 := Roots(h);
    for rr1 in R1 do
      r1 := rr1[1];
      _, h := IsUnivariate(Evaluate(H[1],[u0,r1,r2]));
      R0 := Roots(h);
      for rr0 in R0 do
        r0 := rr0[1];
        if (Im(r0) eq 0 and Im(r1) eq 0 and Im(r2) eq 0) then
          realRoots := Append(realRoots,[r0,r1,r2,1]);
        else
          complexRoot := Append(complexRoots,[r0,r1,r2,1]);
        end if;
      end for;
    end for;
  end for;
end for;


print "trimming precision of output";
realRoots_out := [[roundingMap2(r[1]),roundingMap2(r[2]),roundingMap2(r[3])]: r in realRoots];
complexRoots_out := [[roundingMap2(r[1]),roundingMap2(r[2]),roundingMap2(r[3])]: r in complexRoots];
