text: # -*- coding: utf-8 -*-
  r"""
  Projective plane conics over a rational function field
  
  The class :class:`ProjectiveConic_rational_function_field` represents a
  projective plane conic over a rational function field `F(t)`, where `F`
  is any field. Instances can be created using :func:`Conic`.
  
  AUTHORS:
  
  - Lennart Ackermans (2016-02-07): initial version
      
  EXAMPLES:
      
  Create a conic::
  
      sage: K = FractionField(PolynomialRing(QQ, 't'))
      sage: P.<X, Y, Z> = K[]
      sage: Conic(X^2 + Y^2 - Z^2)
      Projective Conic Curve over Fraction Field of Univariate
      Polynomial Ring in t over Rational Field defined by
      X^2 + Y^2 - Z^2
  
  Points can be found using :meth:`has_rational_point`::
  
      sage: K.<t> = FractionField(QQ['t'])
      sage: C = Conic([1,-t,t])
      sage: C.has_rational_point(point = True)
      (True, (0 : 1 : 1))
  """
  
  #*****************************************************************************
  #       Copyright (C) 2016 Lennart Ackermans
  #
  # This program is free software: you can redistribute it and/or modify
  # it under the terms of the GNU General Public License as published by
  # the Free Software Foundation, either version 2 of the License, or
  # (at your option) any later version.
  #                  http://www.gnu.org/licenses/
  #*****************************************************************************
  from __future__ import division
  from __future__ import absolute_import
  
  from sage.rings.all import PolynomialRing
  from sage.matrix.constructor import diagonal_matrix, matrix, block_matrix
  from sage.schemes.plane_conics.con_field import ProjectiveConic_field
  from sage.arith.all import lcm, gcd
  from sage.modules.free_module_element import vector
  from sage.rings.fraction_field import is_FractionField
  
  class ProjectiveConic_rational_function_field(ProjectiveConic_field):
      r"""
      Create a projective plane conic curve over a rational function field
      `F(t)`, where `F` is any field.
  
      The algorithms used in this class come mostly from [HC2006]_.
  
      EXAMPLES::
  
          sage: K = FractionField(PolynomialRing(QQ, 't'))
          sage: P.<X, Y, Z> = K[]
          sage: Conic(X^2 + Y^2 - Z^2)
          Projective Conic Curve over Fraction Field of Univariate
          Polynomial Ring in t over Rational Field defined by
          X^2 + Y^2 - Z^2
  
      TESTS::
  
          sage: K = FractionField(PolynomialRing(QQ, 't'))
          sage: Conic([K(1), 1, -1])._test_pickling()
  
      REFERENCES:
  
      .. [HC2006] Mark van Hoeij and John Cremona, Solving Conics over
         function fields. J. Théor. Nombres Bordeaux, 2006.
      .. [ACKERMANS2016] Lennart Ackermans, Oplosbaarheid van Kegelsneden.
         http://www.math.leidenuniv.nl/nl/theses/Bachelor/.
      """
      def __init__(self, A, f):
          r"""
          See ``Conic`` for full documentation.
  
          EXAMPLES::
  
              sage: c = Conic([1, 1, 1]); c
              Projective Conic Curve over Rational Field defined by
              x^2 + y^2 + z^2
          """
          ProjectiveConic_field.__init__(self, A, f)
       
      def has_rational_point(self, point = False, algorithm = 'default',
          read_cache = True):
          r"""
          Returns True if and only if the conic ``self``
          has a point over its base field `F(t)`, which is a field of rational
          functions.
  
          If ``point`` is True, then returns a second output, which is
          a rational point if one exists.
  
          Points are cached whenever they are found. Cached information
          is used if and only if ``read_cache`` is True.
          
          The default algorithm does not (yet) work for all base fields `F`.
          In particular, sage is required to have:
          
          * an algorithm for finding the square root of elements in finite
            extensions of `F`;
          
          * a factorization and gcd algorithm for `F[t]`;
          
          * an algorithm for solving conics over `F`.
          
          ALGORITHM:
          
          The parameter ``algorithm`` specifies the algorithm
          to be used:
  
          * ``'default'`` -- use a native Sage implementation, based on the
            algorithm Conic in [HC2006]_.
  
          * ``'magma'`` (requires Magma to be installed) --
            delegates the task to the Magma computer algebra
            system.
          
          EXAMPLES:
          
          We can find points for function fields over (extensions of) `\QQ`
          and finite fields::
          
              sage: K.<t> = FractionField(PolynomialRing(QQ, 't'))
              sage: C = Conic(K, [t^2-2, 2*t^3, -2*t^3-13*t^2-2*t+18])
              sage: C.has_rational_point(point=True)
              (True, (-3 : (t + 1)/t : 1))
              sage: R.<t> = FiniteField(23)[]
              sage: C = Conic([2, t^2+1, t^2+5])
              sage: C.has_rational_point()
              True
              sage: C.has_rational_point(point=True)
              (True, (5*t : 8 : 1))
              sage: F.<i> = QuadraticField(-1)
              sage: R.<t> = F[]
              sage: C = Conic([1,i*t,-t^2+4])
              sage: C.has_rational_point(point = True)
              (True, (-t - 2*i : -2*i : 1))
  
          It works on non-diagonal conics as well::
  
              sage: K.<t> = QQ[]
              sage: C = Conic([4, -4, 8, 1, -4, t + 4])
              sage: C.has_rational_point(point=True)
              (True, (1/2 : 1 : 0))
  
          If no point exists output still depends on the argument ``point``::
  
              sage: K.<t> = QQ[]
              sage: C = Conic(K, [t^2, (t-1), -2*(t-1)])
              sage: C.has_rational_point()
              False
              sage: C.has_rational_point(point=True)
              (False, None)
          
          Due to limitations in Sage of algorithms we depend on, it is not
          yet possible to find points on conics over multivariate function fields
          (see the requirements above)::
          
              sage: F.<t1> = FractionField(QQ['t1'])
              sage: K.<t2> = FractionField(F['t2'])
              sage: a = K(1)
              sage: b = 2*t2^2+2*t1*t2-t1^2
              sage: c = -3*t2^4-4*t1*t2^3+8*t1^2*t2^2+16*t1^3-t2-48*t1^4
              sage: C = Conic([a,b,c])
              sage: C.has_rational_point()
              Traceback (most recent call last):
              ...
              NotImplementedError: is_square() not implemented for elements of
              Univariate Quotient Polynomial Ring in tbar over Fraction Field
              of Univariate Polynomial Ring in t1 over Rational Field with
              modulus tbar^2 + t1*tbar - 1/2*t1^2
          
          In some cases, the algorithm requires us to be
          able to solve conics over `F`. In particular, the following does not
          work::
  
              sage: P.<u> = QQ[]
              sage: E = P.fraction_field()
              sage: Q.<Y> = E[]
              sage: F.<v> = E.extension(Y^2 - u^3 - 1)
              sage: R.<t> = F[]
              sage: K = R.fraction_field()
              sage: C = Conic(K, [u, v, 1])
              sage: C.has_rational_point()
              Traceback (most recent call last):
              ...
              NotImplementedError: has_rational_point not implemented for conics
              over base field Univariate Quotient Polynomial Ring in v over
              Fraction Field of Univariate Polynomial Ring in u over Rational
              Field with modulus v^2 - u^3 - 1
  
          ``has_rational_point`` fails for some conics over function fields
          over finite fields, due to :trac:`20003`::
  
              sage: K.<t> = PolynomialRing(GF(7))
              sage: C = Conic([5*t^2+4, t^2+3*t+3, 6*t^2+3*t+2, 5*t^2+5, 4*t+3, 4*t^2+t+5])
              sage: C.has_rational_point()
              Traceback (most recent call last):
              ...
              TypeError: self (=Scheme morphism:
                From: Projective Conic Curve over Fraction Field of Univariate Polynomial Ring in t over Finite Field of size 7 defined by (5*t^2 + 4)*x^2 + ((6*t^3 + 3*t^2 + 5*t + 5)/(t + 3))*y^2 + ((6*t^6 + 3*t^5 + t^3 + 6*t^2 + 6*t + 2)/(t^4 + t^3 + 4*t^2 + 3*t + 1))*z^2
                To:   Projective Conic Curve over Fraction Field of Univariate Polynomial Ring in t over Finite Field of size 7 defined by (5*t^2 + 4)*x^2 + (t^2 + 3*t + 3)*x*y + (5*t^2 + 5)*y^2 + (6*t^2 + 3*t + 2)*x*z + (4*t + 3)*y*z + (4*t^2 + t + 5)*z^2
                Defn: Defined on coordinates by sending (x : y : z) to
                      (x + ((2*t + 5)/(t + 3))*y + ((3*t^4 + 2*t^3 + 5*t^2 + 5*t + 3)/(t^4 + t^3 + 4*t^2 + 3*t + 1))*z : y + ((6*t^3 + 6*t^2 + 3*t + 6)/(t^3 + 4*t^2 + 2*t + 2))*z : z)) domain must equal right (=Scheme morphism:
                From: Projective Conic Curve over Fraction Field of Univariate Polynomial Ring in t over Finite Field of size 7 defined by (5*t^3 + 6*t^2 + 3*t + 3)*x^2 + (t + 4)*y^2 + (6*t^7 + 2*t^5 + t^4 + 2*t^3 + 3*t^2 + 6*t + 6)*z^2
                To:   Projective Conic Curve over Fraction Field of Univariate Polynomial Ring in t over Finite Field of size 7 defined by (5/(t^3 + 4*t^2 + 2*t + 2))*x^2 + (1/(t^3 + 3*t^2 + 5*t + 1))*y^2 + ((6*t^6 + 3*t^5 + t^3 + 6*t^2 + 6*t + 2)/(t^9 + 5*t^8 + t^7 + 6*t^6 + 3*t^5 + 4*t^3 + t^2 + 5*t + 3))*z^2
                Defn: Defined on coordinates by sending (x : y : z) to
                      ((t^3 + 4*t^2 + 2*t + 2)*x : (t^2 + 5)*y : (t^5 + 4*t^4 + t^2 + 3*t + 3)*z)) codomain
  
              
          TESTS::
  
              sage: K.<t> = FractionField(PolynomialRing(QQ, 't'))
              sage: a = (2*t^2 - 3/2*t + 1)/(37/3*t^2 + t - 1/4)
              sage: b = (1/2*t^2 + 1/3)/(-73*t^2 - 2*t + 11/4)
              sage: c = (6934/3*t^6 + 8798/3*t^5 - 947/18*t^4 + 3949/9*t^3 + 20983/18*t^2 + 28/3*t - 131/3)/(-2701/3*t^4 - 293/3*t^3 + 301/6*t^2 + 13/4*t - 11/16)
              sage: C = Conic([a,b,c])
              sage: C.has_rational_point(point=True)
              (True, (4*t + 4 : 2*t + 2 : 1))
  
          A long time test::
  
              sage: K.<t> = FractionField(PolynomialRing(QQ, 't'))
              sage: a = (-1/3*t^6 - 14*t^5 - 1/4*t^4 + 7/2*t^2 - 1/2*t - 1)/(24/5*t^6 - t^5 - 1/4*t^4 + t^3 - 3*t^2 + 8/5*t + 5)
              sage: b = (-3*t^3 + 8*t + 1/2)/(-1/3*t^3 + 3/2*t^2 + 1/12*t + 1/2)
              sage: c = (1232009/225*t^25 - 1015925057/8100*t^24 + 1035477411553/1458000*t^23 + 7901338091/30375*t^22 - 1421379260447/729000*t^21 + 266121260843/972000*t^20 + 80808723191/486000*t^19 - 516656082523/972000*t^18 + 21521589529/40500*t^17 + 4654758997/21600*t^16 - 20064038625227/9720000*t^15 - 173054270347/324000*t^14 + 536200870559/540000*t^13 - 12710739349/50625*t^12 - 197968226971/135000*t^11 - 134122025657/810000*t^10 + 22685316301/120000*t^9 - 2230847689/21600*t^8 - 70624099679/270000*t^7 - 4298763061/270000*t^6 - 41239/216000*t^5 - 13523/36000*t^4 + 493/36000*t^3 + 83/2400*t^2 + 1/300*t + 1/200)/(-27378/125*t^17 + 504387/500*t^16 - 97911/2000*t^15 + 1023531/4000*t^14 + 1874841/8000*t^13 + 865381/12000*t^12 + 15287/375*t^11 + 6039821/6000*t^10 + 599437/1500*t^9 + 18659/250*t^8 + 1218059/6000*t^7 + 2025127/3000*t^6 + 1222759/6000*t^5 + 38573/200*t^4 + 8323/125*t^3 + 15453/125*t^2 + 17031/500*t + 441/10)
              sage: C = Conic([a,b,c])
              sage: C.has_rational_point(point = True) # long time (4 seconds)
              (True,
               ((-2/117*t^8 + 304/1053*t^7 + 40/117*t^6 - 1/27*t^5 - 110/351*t^4 - 2/195*t^3 + 11/351*t^2 + 1/117)/(t^4 + 2/39*t^3 + 4/117*t^2 + 2/39*t + 14/39) : -5/3*t^4 + 19*t^3 : 1))
          """
          from .constructor import Conic
          
          if read_cache:
              if self._rational_point is not None:
                  return (True, self._rational_point) if point else True
          
          if algorithm != 'default':
              return ProjectiveConic_field.has_rational_point(self, point,
                  algorithm, read_cache)
          
          # Default algorithm
          if self.base_ring().characteristic() == 2:
              raise NotImplementedError("has_rational_point not implemented \
  for function field of characteristic 2.")
          new_conic, transformation, inverse = self.diagonalization()
          coeff = new_conic.coefficients()
          if coeff[0] == 0:
              return (True, transformation([1,0,0])) if point else True
          elif coeff[3] == 0:
              return (True, transformation([0,1,0])) if point else True
          elif coeff[5] == 0:
              return (True, transformation([0,0,1])) if point else True
          
          # We save the coefficients of the reduced form in coeff
          # A zero of the reduced conic can be multiplied by multipliers
          # to get a zero of the old conic
          (coeff, multipliers) = new_conic._reduce_conic()
          new_conic = Conic(coeff)
          transformation = transformation \
              * new_conic.hom(diagonal_matrix(multipliers))
          if coeff[0].degree() % 2 == coeff[1].degree() % 2 and \
                  coeff[1].degree() % 2 == coeff[2].degree() % 2:
              case = 0
          else:
              case = 1
          
          t, = self.base_ring().base().gens() # t in F[t]
          supp = []
          roots = [[], [], []]
          remove = None
          # loop through the coefficients and find a root of f_i (as in
          # [HC2006]) modulo each element in the coefficients' support
          for i in (0,1,2):
              supp.append(list(coeff[i].factor()))
              for p in supp[i]:
                  if p[1] != 1:
                      raise ValueError("Expected factor of exponent 1.")
                  # Convert to monic factor
                  x = p[0]/list(p[0])[-1]
                  N = p[0].base_ring().extension(x, 'tbar')
                  R = PolynomialRing(N, 'u')
                  u, = R.gens()
                  # If p[0] has degree 1, sage might forget the "defining
                  # polynomial" of N, so we define our own modulo operation
                  if p[0].degree() == 1:
                      mod = t.parent().hom([-x[0]])
                  else:
                      mod = N
                  if i == 0:
                      x = -mod(coeff[2])/mod(coeff[1])
                  elif i == 1:
                      x = -mod(coeff[0])/mod(coeff[2])
                  else:
                      x = -mod(coeff[1])/mod(coeff[0])
                  if x.is_square():
                      root = N(x.sqrt())
                  else:
                      return (False, None) if point else False
                  # if case == 0 and p[0] has degree 1, we switch to case
                  # 1 and remove this factor out of the support. In [HC2006]
                  # this is done later, in FindPoint.
                  if case == 0 and p[0].degree() == 1:
                      case = 1
                      # remove later so the loop iterator stays in place.
                      remove = (i,p)
                  else:
                      roots[i].append(root)
          if remove:
              supp[remove[0]].remove(remove[1])
          supp = [[p[0] for p in supp[i]] for i in (0,1,2)]
  
          if case == 0:
          # Find a solution of (5) in [HC2006]
              leading_conic = Conic(self.base_ring().base_ring(),
                          [coeff[0].leading_coefficient(),
                          coeff[1].leading_coefficient(),
                          coeff[2].leading_coefficient()])
              has_point = leading_conic.has_rational_point(True)
              if has_point[0]:
                  if point:
                      pt = new_conic.find_point(supp, roots, case,
                          has_point[1])
                  else:
                      pt = True
                  return (True, transformation(pt)) if point else True
              else:
                  return (False, None) if point else False
          # case == 1:
          if point:
              pt = new_conic.find_point(supp, roots, case)
          else:
              pt = True
          return (True, transformation(pt)) if point else True
              
          
      
      def _reduce_conic(self):
          r"""
          Return the reduced form of the conic, i.e. a conic with base field
          `K=F(t)` and coefficients `a,b,c` such that `a,b,c \in F[t]`,
          `\gcd(a,b)=\gcd(b,c)=\gcd(c,a)=1` and `abc` is square-free.
          
          Assumes `self` is in diagonal form.
          
          OUTPUT:
          
          A tuple (coefficients, multipliers), the coefficients of the conic
          in reduced form and multipliers `\lambda, \mu, \nu \in F(t)^*` such
          that `(x,y,z) \in F(t)` is a solution of the reduced conic if and only
          if `(\lambda x, \mu y, \nu z)` is a solution of `self`.
          
          ALGORITHM:
          
          The algorithm used is the algorithm ReduceConic in [HC2006]_.
          
          EXAMPLES::
          
              sage: K.<t> = FractionField(PolynomialRing(QQ, 't'))
              sage: C = Conic(K, [t^2-2, 2*t^3, -2*t^3-13*t^2-2*t+18])
              sage: C._reduce_conic()
              ([t^2 - 2, 2*t, -2*t^3 - 13*t^2 - 2*t + 18], [t, 1, t])
          """
          
          # start with removing fractions
          coeff = [self.coefficients()[0], self.coefficients()[3],
                  self.coefficients()[5]]
          coeff = lcm(lcm(coeff[0].denominator(), coeff[1].denominator()),
                  coeff[2].denominator()) * vector(coeff)
          # go to base ring of fraction field
          coeff = [self.base().base()(x) for x in coeff]
          coeff = vector(coeff) / gcd(coeff)
          # remove common divisors
          labda = mu = nu = 1
          g1 = g2 = g3 = 0
          ca, cb, cc = coeff
          while g1 != 1 or g2 != 1 or g3 != 1:
              g1 = gcd(ca,cb); ca = ca/g1; cb = cb/g1; cc = cc*g1; nu = g1*nu
              g2 = gcd(ca,cc); ca = ca/g2; cc = cc/g2; cb = cb*g2; mu = g2*mu
              g3 = gcd(cb,cc); cb = cb/g3; cc = cc/g3; ca = ca*g3;
              labda = g3*labda
          coeff = [ca, cb, cc]
          multipliers = [labda, mu, nu]
          
          # remove squares
          for i, x in enumerate(coeff):
              if is_FractionField(x.parent()):
                  # go to base ring of fraction field
                  x = self.base().base()(x)
              
              try:
                  decom = x.squarefree_decomposition()
              except (NotImplementedError, AttributeError):
                  decom = x.factor()
              x = decom.unit()
              x2 = 1
              for factor in decom:
                  if factor[1] > 1:
                      if factor[1] % 2 == 0:
                          x2 *= factor[0] ** (factor[1] // 2)
                      else:
                          x *= factor[0]
                          x2 *= factor[0] ** ((factor[1] - 1) // 2)
                  else:
                      x *= factor[0]
              for j, y in enumerate(multipliers):
                  if j != i:
                      multipliers[j] = y * x2
              coeff[i] = self.base_ring().base().coerce(x);
          
          return (coeff, multipliers)
      
      def find_point(self, supports, roots, case, solution = 0):
          r"""
          Given a solubility certificate like in [HC2006]_, find a point on
          ``self``. Assumes ``self`` is in reduced form (see [HC2006]_ for a
          definition).
  
          If you don't have a solubility certificate and just want to find a
          point, use the function :meth:`has_rational_point` instead.
          
          INPUT:
          
          - ``self`` -- conic in reduced form.
          - ``supports`` -- 3-tuple where ``supports[i]`` is a list of all monic
            irreducible `p \in F[t]` that divide the `i`'th of the 3 coefficients.
          - ``roots`` -- 3-tuple containing lists of roots of all elements of
            ``supports[i]``, in the same order.
          - ``case`` -- 1 or 0, as in [HC2006]_.
          - ``solution`` -- (default: 0) a solution of (5) in [HC2006]_, if
            case = 0, 0 otherwise.
          
          OUTPUT:
          
          A point `(x,y,z) \in F(t)` of ``self``. Output is undefined when the
          input solubility certificate is incorrect.
  
          ALGORITHM:
          
          The algorithm used is the algorithm FindPoint in [HC2006]_, with
          a simplification from [ACKERMANS2016]_.
          
          EXAMPLES::
              
              sage: K.<t> = FractionField(QQ['t'])
              sage: C = Conic(K, [t^2-2, 2*t^3, -2*t^3-13*t^2-2*t+18])
              sage: C.has_rational_point(point=True) # indirect test
              (True, (-3 : (t + 1)/t : 1))
  
          Different solubility certificates give different points::
  
              sage: K.<t> = PolynomialRing(QQ, 't')
              sage: C = Conic(K, [t^2-2, 2*t, -2*t^3-13*t^2-2*t+18])
              sage: supp = [[t^2 - 2], [t], [t^3 + 13/2*t^2 + t - 9]]
              sage: tbar1 = QQ.extension(supp[0][0], 'tbar').gens()[0]
              sage: tbar2 = QQ.extension(supp[1][0], 'tbar').gens()[0]
              sage: tbar3 = QQ.extension(supp[2][0], 'tbar').gens()[0]
              sage: roots = [[tbar1 + 1], [1/3*tbar2^0], [2/3*tbar3^2 + 11/3*tbar3 - 3]]
              sage: C.find_point(supp, roots, 1)
              (3 : t + 1 : 1)
              sage: roots = [[-tbar1 - 1], [-1/3*tbar2^0], [-2/3*tbar3^2 - 11/3*tbar3 + 3]]
              sage: C.find_point(supp, roots, 1)
              (3 : -t - 1 : 1)
          """
          Ft = self.base().base()
          F = Ft.base()
          t, = Ft.gens()
          coefficients = [Ft(self.coefficients()[0]), Ft(self.coefficients()[3]),
              Ft(self.coefficients()[5])]
          deg = [coefficients[0].degree(), coefficients[1].degree(),
                  coefficients[2].degree()]
          # definitions as in [HC2006] and [ACKERMANS2016]
          A = ((deg[1] + deg[2]) / 2).ceil() - case
          B = ((deg[2] + deg[0]) / 2).ceil() - case
          C = ((deg[0] + deg[1]) / 2).ceil() - case
          
          # For all roots as calculated by has_rational_point(), we create
          # a system of linear equations. As in [ACKERMANS2016], we do this
          # by calculating the matrices for all phi_p, with basis consisting
          # of monomials of x, y and z in the space V of potential solutions:
          # t^0, ..., t^A, t^0, ..., t^B and t^0, ..., t^C.
          phi = []
          for (i, p) in enumerate(supports[0]):
              # lift to F[t] and map to R, with R as defined above
              if roots[0][i].parent().is_finite():
                  root = roots[0][i].polynomial()
              else:
                  root = roots[0][i].lift()
              alpha = root.parent().hom([t])(root)
              d = p.degree()
              # Calculate y - alpha*z mod p for all basis vectors
              phi_p = [[] for i in range(A+B+C+4)]
              phi_p[0:A+1] = [vector(F, d)] * (A+1)
              phi_p[A+1] = vector(F, d, {0: F(1)})
              lastpoly = F(1)
              for n in range(B):
                  lastpoly = (lastpoly * t) % p
                  phi_p[A+2+n] = vector(F, d, lastpoly.dict())
              lastpoly = -alpha % p
              phi_p[A+B+2] = vector(F, d, lastpoly.dict())
              for n in range(C):
                  lastpoly = (lastpoly * t) % p
                  phi_p[A+B+3+n] = vector(F, d, lastpoly.dict())
              phi_p[A+B+C+3] = vector(F, d)
              phi.append(matrix(phi_p).transpose())
          for (i, p) in enumerate(supports[1]):
              if roots[1][i].parent().is_finite():
                  root = roots[1][i].polynomial()
              else:
                  root = roots[1][i].lift()
              alpha = root.parent().hom([t])(root)
              d = p.degree()
              # Calculate z - alpha*x mod p for all basis vectors
              phi_p = [[] for i in range(A+B+C+4)]
              phi_p[A+1:A+B+2] = [vector(F, d)] * (B+1)
              phi_p[A+B+2] = vector(F, d, {0: F(1)})
              lastpoly = F(1)
              for n in range(C):
                  lastpoly = (lastpoly * t) % p
                  phi_p[A+B+3+n] = vector(F, d, lastpoly.dict())
              lastpoly = -alpha % p
              phi_p[0] = vector(F, d, lastpoly.dict())
              for n in range(A):
                  lastpoly = (lastpoly * t) % p
                  phi_p[1+n] = vector(F, d, lastpoly.dict())
              phi_p[A+B+C+3] = vector(F, d)
              phi.append(matrix(phi_p).transpose())
          for (i, p) in enumerate(supports[2]):
              if roots[2][i].parent().is_finite():
                  root = roots[2][i].polynomial()
              else:
                  root = roots[2][i].lift()
              alpha = root.parent().hom([t])(root)
              d = p.degree()
              # Calculate x - alpha*y mod p for all basis vectors
              phi_p = [[] for i in range(A+B+C+4)]
              phi_p[A+B+2:A+B+C+3] = [vector(F, d)] * (C+1)
              phi_p[0] = vector(F, d, {0: F(1)})
              lastpoly = F(1)
              for n in range(A):
                  lastpoly = (lastpoly * t) % p
                  phi_p[1+n] = vector(F, d, lastpoly.dict())
              lastpoly = -alpha % p
              phi_p[A+1] = vector(F, d, lastpoly.dict())
              for n in range(B):
                  lastpoly = (lastpoly * t) % p
                  phi_p[A+2+n] = vector(F, d, lastpoly.dict())
              phi_p[A+B+C+3] = vector(F, d)
              phi.append(matrix(phi_p).transpose())
          if case == 0:
              # We need three more equations
              lx = Ft(solution[0]).leading_coefficient()
              ly = Ft(solution[1]).leading_coefficient()
              lz = Ft(solution[2]).leading_coefficient()
              phi.append(matrix([vector(F, A+B+C+4, {A:1, A+B+C+3:-lx}),
                  vector(F, A+B+C+4, {A+B+1:1, A+B+C+3:-ly}),
                  vector(F, A+B+C+4, {A+B+C+2: 1, A+B+C+3:-lz})]))
          # Create the final matrix which we will solve
          M = block_matrix(phi, ncols = 1, subdivide = False)
          solution_space = M.right_kernel()
          for v in solution_space.basis():
              if v[:A+B+C+3] != 0: # we don't want to return a trivial solution
                  X = Ft(list(v[:A+1]))
                  Y = Ft(list(v[A+1:A+B+2]))
                  Z = Ft(list(v[A+B+2:A+B+C+3]))
                  return self.point([X,Y,Z])
  
          raise RuntimeError("No solution has been found: possibly incorrect\
   solubility certificate.")
  