class TestHomogeneousWieghts(unittest.TestCase):

  def setUp(self):
    self.poly_ring = PolynomialRing(QQ,"x",3);
    self.x = self.poly_ring.gens()[0];
    self.y = self.poly_ring.gens()[1];
    self.z = self.poly_ring.gens()[2];
  
  def test_homogenous(self):
    x = self.x
    y = self.y
    z = self.z
    divisor = x*y*z + x**3 + z**3
    wieghts = homogenous_wieghts(divisor)
    #Test this works with a homogenous divisor
    self.assertEquals(wieghts,[3,1,1,1]);
    
  def test_weighted_homogenous(self):
    x = self.x
    y = self.y
    z = self.z
    divisor = x**2*y-z**2
    wieghts = homogenous_wieghts(divisor)
    #Test this works with a weighted homogenous divisor
    self.assertEquals(wieghts,[4,1,2,2])
  
  def test_not_homogenous(self):
    x = self.x
    y = self.y
    z = self.z
    divisor = x**2*y-x**3 + z +y**2;
    self.assertRaises(NotWieghtHomogeneousException,homogenous_wieghts,divisor)

def whitney_divisor(var=None):
  if var:
    poly_ring = PolynomialRing(QQ,3,var)
  else:
    poly_ring = PolynomialRing(QQ,3,"xyz")
  gens = poly_ring.gens()
  return gens[0]**2*gens[1]-gens[2]**2

    def genus(self):
        """
        Return the genus of this function field
        For now, the genus is computed using singular

        EXAMPLES::

            sage: K.<x> = FunctionField(QQ); R.<y> = K[]
            sage: L.<y> = K.extension(y^3 - (x^3 + 2*x*y + 1/x))
            sage: L.genus()
            3
        """
        # unfortunately singular can not compute the genus with the polynomial_ring()._singular_
        # object because genus method only accepts a ring of transdental degree 2 over a prime field
        # not a ring of transdental degree 1 over a rational function field of one variable

        if is_RationalFunctionField(self._base_field) and self._base_field.constant_field().is_prime_field():

            #Making the auxiliary ring which only has polynomials with integral coefficients.
            tmpAuxRing = PolynomialRing(self._base_field.constant_field(), str(self._base_field.gen())+','+str(self._ring.gen()))
            intMinPoly, d = self._make_monic_integral(self._polynomial)
            curveIdeal = tmpAuxRing.ideal(intMinPoly)

            singular.lib('normal.lib') #loading genus method in singular
            return int(curveIdeal._singular_().genus())

        else:
            raise NotImplementedError("Computation of genus over this rational function field not implemented yet")

    def quantum_group(self, q=None, c=None):
        r"""
        Return the quantum group of ``self``.

        The corresponding quantum group is the
        :class:`~sage.algebras.lie_algebras.onsager.QuantumOnsagerAlgebra`.
        The parameter `c` must be such that `c(1) = 1`

        INPUT:

        - ``q`` -- (optional) the quantum parameter; the default
          is `q \in R(q)`, where `R` is the base ring of ``self``
        - ``c`` -- (optional) the parameter `c`; the default is ``q``

        EXAMPLES::

            sage: O = lie_algebras.OnsagerAlgebra(QQ)
            sage: Q = O.quantum_group()
            sage: Q
            q-Onsager algebra with c=q over Fraction Field of
             Univariate Polynomial Ring in q over Rational Field
        """
        if q is None:
            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
            q = PolynomialRing(self.base_ring(), 'q').fraction_field().gen()
        if c is None:
            c = q
        else:
            c = q.parent()(c)
        return QuantumOnsagerAlgebra(self, q, c)

 def _interpolate(evaluation, num_of_var, order):
     if num_of_var == 0 or order == 0:
         return evaluation[0]
     base_field = polynomial_ring.base_ring()
     q = base_field.cardinality()
     n_by_q = q**(num_of_var - 1)
     d = min(order + 1, q)
     multipoint_evaluation_list = []
     uni_poly_ring = PolynomialRing(base_field, 'x')
     base_field_zero = base_field.zero()
     for k in range(n_by_q):
         iterator = iter(base_field)
         points = []
         for i in range(d):
             xcoordinate = next(iterator)
             points.append((xcoordinate, evaluation[k + i * n_by_q]))
         polyVector = uni_poly_ring.lagrange_polynomial(
             points).coefficients(sparse=False)
         if len(polyVector) < d:
             # adding zeros to represet a (d-1) degree polynomial
             polyVector += [base_field_zero] * (d - len(polyVector))
         multipoint_evaluation_list.append(polyVector)
     poly = polynomial_ring.zero()
     z = 1
     x = polynomial_ring.gen(num_of_var - 1)
     for k in range(d):  # computing the polynomial
         poly = poly + z * _interpolate([multipoint_evaluation_list[i][k]
                                         for i in range(n_by_q)], num_of_var - 1, order - k)
         z *= x
     return poly

def conv(N):
    file = "data/%s"%N
    if not os.path.exists(file):
        raise RuntimeError, "Data for level %s does not exist."%N

    F = open(file).read()
    i = F.find(":=")
    if i == -1:
        raise RuntimeError, "Syntax error in file for level %s."%N
    F = F[i+2:]
    TRANS = [("[*", "["), ("*]", "]"), ("<","["), (">","]"), \
             (";",""), ("\n",""), ("^","**")]
    for z,w in TRANS:
        F = F.replace(z,w)
    X = []
    # Define x so the eval below works.
    R = PolynomialRing(RationalField())
    x = R.gen()
    print "starting eval."
    #print "F = ", F
    for f in eval(F):
        print "creating object from f=",f[:4]
        cp = {}
        disc = 0        
        for z in f[5]:
            g = R(z[1])
            disc = GCD(disc,g.discriminant())
            cp[z[0]] = g
        X.append(ModularForm(f[0],f[1],f[2],f[3],f[4],cp,disc))
    return X

class TestConvertSymToPoly(unittest.TestCase):

  def setUp(self):
    self.poly_ring = PolynomialRing(QQ,"x",3);
    self.x = self.poly_ring.gens()[0];
    self.y = self.poly_ring.gens()[1];
    self.z = self.poly_ring.gens()[2];
    self.vars = var('x,y,z')
    
  def test_zero(self):
    zero = 0*self.vars[0]
    poly = convert_symbolic_to_polynomial(zero,self.poly_ring,self.vars)
    self.assertEqual(poly,self.poly_ring.zero())
    
  def test_convert(self):
    x = self.x
    y = self.y
    z = self.z
    sym_poly = 4*self.vars[0]**4 + self.vars[1]*self.vars[0]**12*self.vars[1]-self.vars[2]
    poly = convert_symbolic_to_polynomial(sym_poly,self.poly_ring,self.vars)
    self.assertEqual(poly,4*x**4+y*x**12*y-z)

  def test_convert_univarient(self):
    y = self.y
    sym_poly = 4*self.vars[1]**4 + self.vars[1]*self.vars[1]**12*self.vars[1]-self.vars[1]
    poly = convert_symbolic_to_polynomial(sym_poly,self.poly_ring,self.vars)
    self.assertEqual(poly,4*y**4+y*y**12*y-y)

  def test_convert_partial(self):
    y = self.y
    z = self.z
    sym_poly = 4*self.vars[2]**4 + self.vars[1]*self.vars[1]**12*self.vars[1]-self.vars[1]
    poly = convert_symbolic_to_polynomial(sym_poly,self.poly_ring,self.vars)
    self.assertEqual(poly,4*z**4+y*y**12*y-y)

def automorphy_factor_vector(p, a, c, k, chi, p_prec, var_prec, R):
    """
    EXAMPLES::
        
        sage: from sage.modular.pollack_stevens.families_util import automorphy_factor_vector
        sage: automorphy_factor_vector(3, 1, 3, 0, None, 4, 3, PowerSeriesRing(ZpCA(3), 'w'))
        [1 + O(3^20), O(3^21) + (3 + 3^2 + 2*3^3 + O(3^21))*w + (3^2 + 2*3^3 + O(3^22))*w^2, O(3^22) + (3^2 + 2*3^3 + O(3^22))*w + (2*3^2 + O(3^22))*w^2, O(3^22) + (3^2 + 3^3 + O(3^22))*w + (2*3^3 + O(3^23))*w^2]
    """
    S = PolynomialRing(R, 'z')
    z = S.gens()[0]
    w = R.gen()
    aut = S(1)
    for n in range(1, var_prec):
        LB = logpp_binom(n, p, p_prec)
        ta = ZZ(Qp(p, 2 * max(p_prec, var_prec)).teichmuller(a))
        arg = (a / ta - 1) / p + c / (p * ta) * z
        aut += LB(arg) * (w ** n)
    aut *= (ta ** k)
    if not (chi is None):
        aut *= chi(a)
    aut = aut.list()
    len_aut = len(aut)
    if len_aut == p_prec:
        return aut
    elif len_aut > p_prec:
        return aut[:p_prec]
    return aut + [R.zero_element()] * (p_prec - len_aut)

    def generator_relations(self, K) :
        """
        An ideal `I` in the attach polynomial ring `R`, such that the ring of
        modular forms is a quotient of `R / I`. This ideal must be unique for `K`.
        
        INPUT:
            - `K`      -- A ring or module; The ring of Fourier coefficients.
        
        OUTPUT:
            An ideal in a polynomial ring.
            
        TESTS::
            sage: from hermitianmodularforms import *
            sage: HermitianModularFormD2_Gamma(-3).generator_relations(QQ)
            Ideal (0) of Multivariate Polynomial Ring in HE4, HE6, HE10, HE12, Hphi9 over Rational Field
            sage: HermitianModularFormD2_Gamma(-3).generator_relations(GF(2))
            Traceback (most recent call last):
            ...
            NotImplementedError: Only Fourier coefficients in a number fields are implemented.
        """
        if self.__D == -3 :
            if K is QQ or K in NumberFields() :
                R = PolynomialRing(K, self._generator_names(K))
                return R.ideal(0)

            raise NotImplementedError( "Only Fourier coefficients in a number fields are implemented." )
        
        raise NotImplementedError( "Discriminant %s is not implemented." % (self.__D,) )

class TestMonomialsOfOrder(unittest.TestCase):
 
  def setUp(self):
    self.poly_ring = PolynomialRing(QQ,"x",3);
    self.x = self.poly_ring.gens()[0];
    self.y = self.poly_ring.gens()[1];
    self.z = self.poly_ring.gens()[2];
  
  def test_zero(self):
    mons = [mon for mon in monomials_of_order(0,self.poly_ring,[1,1,1])]
    self.assertEqual(mons,[self.poly_ring.one()])

  def test_homogeneous_3(self):
    x = self.x;
    y = self.y;
    z = self.z;
    true_mons = Set([x**3,y**3,z**3,x**2*y,x**2*z,y**2*x,y**2*z,z**2*x,z**2*y,x*y*z])
    mons = [mon for mon in monomials_of_order(3,self.poly_ring,[1,1,1])]
    self.assertEqual(true_mons,Set(mons))

  def test_non_homogeneous_4(self):
    x = self.x;
    y = self.y;
    z = self.z;
    true_mons = Set([x**4,x**2*y,x*z,y**2])
    mons = [mon for mon in monomials_of_order(4,self.poly_ring,[1,2,3])]
    self.assertEqual(true_mons,Set(mons))

class TestWeightedMinDegree(unittest.TestCase):

  def setUp(self):
    self.poly_ring = PolynomialRing(QQ,"x",3);
    self.x = self.poly_ring.gens()[0];
    self.y = self.poly_ring.gens()[1];
    self.z = self.poly_ring.gens()[2];

  def test_zero(self):
    self.assertEqual(wieghted_min_degree(self.poly_ring.zero(),[1,1,1]),0);
  
  def test_homogeneous(self):
    x = self.x;
    y = self.y;
    z = self.z;
    f = x**4 + 4*y*z**3 - z**2;
    self.assertEqual(wieghted_min_degree(f,[1,1,1]),2);

  def test_non_homogeneous(self):
    x = self.x;
    y = self.y;
    z = self.z;
    f = x**4 + 4*y*z**3 - z**2;
    self.assertEqual(wieghted_min_degree(f,[2,1,2]),4);

  def test_negative_wieghts(self):
    x = self.x;
    y = self.y;
    z = self.z;
    f = x**4 + 4*y*z**3 - z**2;
    self.assertEqual(wieghted_min_degree(f,[-1,-1,1]),-4);

    def _roots_univariate_polynomial(self, p, ring=None, multiplicities=None, algorithm=None):
        r"""
        Return a list of pairs ``(root,multiplicity)`` of roots of the polynomial ``p``.

        If the argument ``multiplicities`` is set to ``False`` then return the
        list of roots.

        .. SEEALSO::

            :meth:`_factor_univariate_polynomial`

        EXAMPLES::

            sage: R.<x> = PolynomialRing(GF(5),'x')
            sage: K = GF(5).algebraic_closure('t')

            sage: sorted((x^6 - 1).roots(K,multiplicities=False))
            [1, 4, 2*t2 + 1, 2*t2 + 2, 3*t2 + 3, 3*t2 + 4]
            sage: ((K.gen(2)*x - K.gen(3))**2).roots(K)
            [(3*t6^5 + 2*t6^4 + 2*t6^2 + 3, 2)]

            sage: for _ in range(10):
            ....:     p = R.random_element(degree=randint(2,8))
            ....:     for r in p.roots(K, multiplicities=False):
            ....:         assert p(r).is_zero(), "r={} is not a root of p={}".format(r,p)

        """
        from sage.arith.all import lcm
        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

        # first build a polynomial over some finite field
        coeffs = [v.as_finite_field_element(minimal=True) for v in p.list()]
        l = lcm([c[0].degree() for c in coeffs])
        F, phi = self.subfield(l)
        P = p.parent().change_ring(F)

        new_coeffs = [self.inclusion(c[0].degree(), l)(c[1]) for c in coeffs]

        polys = [(g,m,l,phi) for g,m in P(new_coeffs).factor()]
        roots = []    # a list of pair (root,multiplicity)
        while polys:
            g,m,l,phi = polys.pop()
        
            if g.degree() == 1: # found a root
                r = phi(-g.constant_coefficient())
                roots.append((r,m))
            else: # look at the extension of degree g.degree() which contains at
                  # least one root of g
                ll = l * g.degree()
                psi = self.inclusion(l, ll)
                FF, pphi = self.subfield(ll)
                # note: there is no coercion from the l-th subfield to the ll-th
                # subfield. The line below does the conversion manually.
                g = PolynomialRing(FF, 'x')([psi(_) for _ in g])
                polys.extend((gg,m,ll,pphi) for gg,_ in g.factor())

        if multiplicities:
            return roots
        else:
            return [r[0] for r in roots]

def rand_w_hom_divisor(n,degs=None,mon_num=None,var="z"):
  if degs==None:
    degs = [randrange(2,6) for _ in range(n)]
  deg = sum(degs)
  if mon_num==None:
    mon_num = randrange(2,8)
  poly_ring = PolynomialRing(QQ,n,var)
  div = poly_ring.zero()
  min_w = min(degs)
  for i in range(mon_num):
    expo = [0]*n
    cur_deg = 0
    while cur_deg!=deg:
      if cur_deg>deg:
        expo = [0]*n
        cur_deg = 0
      if deg-cur_deg<min_w:
        expo = [0]*n
        cur_deg = 0
      next_g = randrange(0,n)
      expo[next_g] += 1
      cur_deg += degs[next_g]
    coeff = randrange(-n,n)/n
    mon = poly_ring.one()
    for i,e in enumerate(expo):
      mon *= poly_ring.gens()[i]**e
    div += coeff*mon
  return div

def breiskorn_pham_divisor(n,var="z"):
  poly_ring = PolynomialRing(QQ,n,var)
  for i,g in enumerate(poly_ring.gens()):
    if i==0:
      div = g**3
    else:
      div += g**2
  return div

def comp_prod(P, Q):
    n = Q.degree()
    K = P.base_ring()
    A = PolynomialRing(K, 'X')
    X = A.gen()
    AA = PolynomialRing(K, 'Y,Z')
    Y, Z = AA.gens()
    return P(Y).resultant(AA(Y**n * Q(Z/Y)), Y)(1,X)

def braid_divisor(n,var="z"):
  poly_ring = PolynomialRing(QQ,n,var)
  div = poly_ring.one()
  gens = poly_ring.gens()
  for i in range(n):
    for j in range(i+1,n):
      div *= (gens[i]-gens[j])
  return div

 def test_p_module_n_crossing(self):
   #Make sure this doesnt throw an error - fix bug
   for i in range(4,5):
     p_ring = PolynomialRing(QQ,i,"z")
     crossing = p_ring.one()
     for g in p_ring.gens():
       crossing *= g
     logdf = LogarithmicDifferentialForms(crossing)
     logdf.p_module(i-1)

 def eu(p):
     """
     local euler factor
     """
     f = rho.local_factor(p)
     co = [ZZ(round(x)) for x in f.coeffs()]
     R = PolynomialRing(QQ, "T")
     T = R.gens()[0]
     return sum( co[n] * T**n for n in range(len(co)))

def BezoutianQuadraticForm(f, g):
    r"""
    Compute the Bezoutian of two polynomials defined over a common base ring.  This is defined by

    .. MATH::

        {\rm Bez}(f, g) := \frac{f(x) g(y) - f(y) g(x)}{y - x}

    and has size defined by the maximum of the degrees of `f` and `g`.

    INPUT:

    - `f`, `g` -- polynomials in `R[x]`, for some ring `R`

    OUTPUT:

    a quadratic form over `R`

    EXAMPLES::

        sage: R = PolynomialRing(ZZ, 'x')
        sage: f = R([1,2,3])
        sage: g = R([2,5])
        sage: Q = BezoutianQuadraticForm(f, g) ; Q
        Quadratic form in 2 variables over Integer Ring with coefficients:
        [ 1 -12 ]
        [ * -15 ]

    AUTHORS:

    - Fernando Rodriguez-Villegas, Jonathan Hanke -- added on 11/9/2008

    """
    ## Check that f and g are polynomials with a common base ring
    if not is_Polynomial(f) or not is_Polynomial(g):
        raise TypeError("Oops!  One of your inputs is not a polynomial. =(")
    if f.base_ring() != g.base_ring():                   ## TO DO:  Change this to allow coercion!
        raise TypeError("Oops!  These polynomials are not defined over the same coefficient ring.")

    ## Initialize the quadratic form
    R = f.base_ring()
    P = PolynomialRing(R, ['x','y'])
    a, b = P.gens()
    n = max(f.degree(), g.degree())
    Q = QuadraticForm(R, n)

    ## Set the coefficients of Bezoutian
    bez_poly = (f(a) * g(b) - f(b) * g(a)) // (b - a)    ## Truncated (exact) division here
    for i in range(n):
        for j in range(i, n):
            if i == j:
                Q[i,j] = bez_poly.coefficient({a:i,b:j})
            else:
                Q[i,j] = bez_poly.coefficient({a:i,b:j}) * 2

    return Q

def _eta_relations_helper(eta1, eta2, degree, qexp_terms, labels, verbose):
    r"""
    Helper function used by eta_poly_relations. Finds a basis for the
    space of linear relations between the first qexp_terms of the
    `q`-expansions of the monomials
    `\eta_1^i * \eta_2^j` for `0 \le i,j < degree`,
    and calculates a Groebner basis for the ideal generated by these
    relations.

    Liable to return meaningless results if qexp_terms isn't at least
    `1 + d*(m_1,m_2)` where

    .. math::

       m_i = min(0, {\text degree of the pole of $\eta_i$ at $\infty$})

    as then 1 will be in the ideal.

    EXAMPLE::

        sage: from sage.modular.etaproducts import _eta_relations_helper
        sage: r,s = EtaGroup(4).basis()
        sage: _eta_relations_helper(r,s,4,100,['a','b'],False)
        [a*b - a + 16]
        sage: _eta_relations_helper(EtaProduct(26, {2:2,13:2,26:-2,1:-2}),EtaProduct(26, {2:4,13:2,26:-4,1:-2}),3,12,['a','b'],False) # not enough terms, will return rubbish
        [1]
    """

    indices = [(i,j) for j in range(degree) for i in range(degree)]
    inf = CuspFamily(eta1.level(), 1)

    pole_at_infinity = -(min([0, eta1.order_at_cusp(inf)]) + min([0,eta2.order_at_cusp(inf)]))*degree
    if verbose: print "Trying all coefficients from q^%s to q^%s inclusive" % (-pole_at_infinity, -pole_at_infinity + qexp_terms - 1)

    rows = []
    for j in xrange(qexp_terms):
        rows.append([])
    for i in indices:
        func = (eta1**i[0]*eta2**i[1]).qexp(qexp_terms)
        for j in xrange(qexp_terms):
            rows[j].append(func[j - pole_at_infinity])
    M = matrix(rows)
    V = M.right_kernel()
    if V.dimension() == 0:
        if verbose: print "No polynomial relation of order %s valid for %s terms" % (degree, qexp_terms)
        return None
    if V.dimension() >= 1:
        #print "Found relation: "
        R = PolynomialRing(QQ, 2, labels)
        x,y = R.gens()
        relations = []
        for c in V.basis():
            relations.append(sum( [ c[v] * x**indices[v][0] * y**indices[v][1] for v in xrange(len(indices))]))
            #print relations[-1], " = 0"
        id = R.ideal(relations)
        return id.groebner_basis()