find_case(p,A,B,C,D,E,F) = {

f = (A*X_0^2 + B*X_0*X_1 + C*X_0*X_2 + D*X_1^2 + E*X_1*X_2 + F*X_2^2)*Mod(1,p);
M = [0];

for(i=0, p-1,
  for(l=0, p-1,
    if(i*l %p == A,
      for(k=0, p-1,
        for(n=0, p-1, 
          if(k*n %p == F,
            for(m=0, p-1,
              for(j=0, p-1,
                if(m*j %p == D, 
                  if( (i*X_0 + j*X_1 + k*X_2)*(l*X_0 + m*X_1 + n*X_2)*Mod(1,p) == f,
                    M = [i,j,k;l,m,n]
                  )
                )
              )
            )
          )
        )
      )
    )
  )
);

if(length(M) == 1,
  return("Intersects every line in 0, 1, or 2 points"),
    for(a = 1, p - 1,
      if(a*M[1,] == M[2,], return("Is a double-line"))
    );
  return("Is the union of two lines")
)
}

find_point(p,A,B,C,D,E,F) = {

if( find_case(p,A,B,C,D,E,F) != "Intersects every line in 0, 1, or 2 points",
  return("Reducible"),
  (f(X_0,X_1,X_2) = (A*X_0^2+B*X_0*X_1+C*X_0*X_2+D*X_1^2+E*X_1*X_2+F*X_2^2)*Mod(1,p));
  for(a=0, p-1,
    for(b=0, p-1,
      for(c=0, p-1,
        if( f(a,b,c)%p == 0 && (a != Mod(0,p) || b != Mod(0,p) || c != Mod(0,p)),
          return([a,b,c])
        )
      )
    )
  )
)
}

   
parametrize(p,A,B,C,D,E,F) = {

if( find_case(p,A,B,C,D,E,F) != "Intersects every line in 0, 1, or 2 points",
  return("Reducible"),
  point = find_point(p,A,B,C,D,E,F);
  (f(X_0,X_1,X_2) = (A*X_0^2+B*X_0*X_1+C*X_0*X_2+D*X_1^2+E*X_1*X_2+F*X_2^2)*Mod(1,p));
  for(i=0, p-1,
    for(j=0, p-1,
      for(k=0, p-1,
        for(l=0, p-1,
          for(m=0, p-1,
            for(n=0, p-1,
              for(r=0, p-1,
                for(s=0, p-1,
                  for(t=0, p-1,
                    if( (i != Mod(0,p) || j != Mod(0,p) || k != Mod(0,p)) && (l != Mod(0,p) || m != Mod(0,p) || n != Mod(0,p)) && (r != Mod(0,p) || s != Mod(0,p) || t != Mod(0,p)) 
                    && f(i*S^2 + j*S*T + k*T^2, l*S^2 + m*S*T + n*T^2, r*S^2 + s*S*T + t*T^2) == 0
                    && gcd(i*S^2 + j*S*T + k*T^2, gcd(l*S^2 + m*S*T + n*T^2, r*S^2 + s*S*T + t*T^2)) == 1,
                      return([point, i*S^2 + j*S*T + k*T^2, l*S^2 + m*S*T + n*T^2, r*S^2 + s*S*T + t*T^2])
                    )
                  )
                )
              )
            )
          )
        )
      )
    )
  )
)
}