Ticket #9334: trac_9334-hilbert.patch

sage/rings/arith.py

@@ -3202 +3202 @@
     r"""
     The Legendre symbol `(x|p)`, for `p` prime.
     
-    .. note::
+    .. NOTE::
 
        The :func:`kronecker_symbol` command extends the Legendre
        symbol to composite moduli and `p=2`.
+       The :func:`generalized_legendre_symbol` command extends the Legendre
+       symbol to number fields.
     
     INPUT:
     
@@ -3241 +3243 @@
     if p == 2:
         raise ValueError, "p must be odd"
     return x.kronecker(p)
+
+def generalized_legendre_symbol(self,P):
+    r"""
+    Returns the generalized Legendre symbol `(a|\mathfrak{p})` for odd primes `\mathfrak{p}.`
+    
+    INPUT:
+
+    -``P`` --odd prime ideal of a number field
+    
+    OUTPUT:
+
+    The generalized legendre symbol.
+
+    ALGORITHM: Comes from an exercise in Cohen: A Course in Computational Number Theory.
+    
+    .. NOTE::
+
+        This is a naive implementation and is possibly VERY SLOW.  
+        The algorithm comes from the definition given in Cohen's A Course in Computational Number Theory:
+        `(a|\mathfrak{p}) = \epsilon`
+        where `\epsilon \in \{-1,+1,0\}`, `N(\mathfrak{p})` is the norm, and 
+        `a^{(N(\mathfrak{p})-1)/2} \equiv \epsilon \mod{\mathfrak{p}}.` 
+        
+    EXAMPLES::
+
+        sage: K.<a> = NumberField(x^2-23)
+        sage: P = K.ideal(5)
+        sage: b = a+5
+        sage: generalized_legendre_symbol(a,P)
+        -1
+        sage: generalized_legendre_symbol(K(7),P)
+        1
+    """
+    if not P.is_prime():
+        raise ValueError, 'P must be a prime ideal'
+    if self in P:
+        return ZZ(0)
+    n = P.norm()
+    K = P.ring()
+    k = P.residue_field() 
+    if n%2 == 0:
+        raise ValueError, 'P must be odd'
+    l = k(K(self)**((n-1)/2))
+    if k(-1) == l:
+        return ZZ(-1)
+    return ZZ(l)
+   
+
     
 def primitive_root(n):
     """
@@ -3872 +3922 @@
     
     
     OUTPUT: integer (0, -1, or 1)
+
+    .. NOTE::
+
+       The :func:`generalized_hilbert_symbol` command extends the Hilbert
+       symbol to number fields.
     
     EXAMPLES::
     
@@ -3948 +4003 @@
     else:
         raise ValueError, "Algorithm %s not defined"%algorithm
 
+def generalized_hilbert_symbol(a,b,P):
+    
+    r"""
+    Returns the Hilbert symbol `(a,b)_\mathfrak{p}` for primes `\mathfrak{p}` over a number field.
+    
+    INPUT:
+
+    - ``a,b`` -- elements of a number field
+    - ``P`` -- odd prime ideal of the same number field
+    
+    OUTPUT:
+    The Hilbert symbol `(a,b)_{\mathfrak{p}}`.
+    
+    .. NOTE::
+
+        The algorithm follows Proposition 5.5 from John Voight: Identifying the Matrix Ring.
+        
+    EXAMPLES::
+
+        sage: K.<a> = NumberField(x^2-23)
+        sage: P = K.ideal(5)
+        sage: b = a+5
+        sage: generalized_hilbert_symbol(a,a+5,K.ideal(5))
+        1
+        sage: generalized_hilbert_symbol(a+1,13,K.ideal(-a-6))
+        -1
+        
+    ::
+
+        sage: K.<a> = NumberField(x^5-23)
+        sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+        sage: b = -a+5
+        sage: generalized_hilbert_symbol(a,a+5,K.ideal(2*a^4 + 3*a^3 + 4*a^2 + 15*a + 11))
+        1
+        sage: generalized_hilbert_symbol(a,b,K.ideal(2*a^4 + 3*a^3 + 4*a^2 + 15*a + 11))
+        1
+        sage: generalized_hilbert_symbol(a,b,P)
+        -1
+        sage: generalized_hilbert_symbol(a,2,P)
+        -1
+
+    ::
+
+        sage: K.<a> = NumberField(x^3+x+1)
+        sage: P = K.ideal(-a^2-a-2)
+        sage: generalized_hilbert_symbol(-6912, 24, P)
+        1
+
+    ::
+
+        sage: K.<a> = NumberField(x^5-23)
+        sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+        sage: Q = K.ideal(-7*a^4 + 13*a^3 - 13*a^2 - 2*a + 50)
+        sage: b = -a+5
+        sage: generalized_hilbert_symbol(a,b,P)
+        -1
+        sage: generalized_hilbert_symbol(a,b,Q)
+        1
+
+    ::
+        sage: K.<i> = QuadraticField(-1)
+        sage: O = K.maximal_order()
+        sage: generalized_hilbert_symbol(K(1/2), K(1), 3*O)
+        1
+        sage: generalized_hilbert_symbol(K(1/3), K(1), (1+i)*O)
+        1
+        sage: p = 1+i
+        sage: generalized_hilbert_symbol(p,p,p*O)
+        -1
+
+    """
+    from sage.rings.number_field.number_field import generalized_even_hilbert_symbol
+
+    if not P.is_prime():
+        raise ValueError, "P must be a prime"
+    
+    pi = P.uniformizer()
+    K = P.ring()
+
+    a = K(a)*K(a.denominator())
+    b = K(b)*K(b.denominator())
+    c = K(-a*b)
+
+    a = a / pi**(2*(a.valuation(P)//2))
+    b = b / pi**(2*(b.valuation(P)//2))
+    c = c / pi**(2*(c.valuation(P)//2))
+
+
+    if a.valuation(P)%2==1 and b.valuation(P)%2==1:
+       return generalized_hilbert_symbol(a, c, P) 
+
+    if P.divides(2):
+            return generalized_even_hilbert_symbol(a,b,P)
+
+   
+    ssa = a.square_symbol(P)
+    if ssa == 1:
+        return 1
+    ssb = b.square_symbol(P)
+    if ssb == 1:
+        return 1
+    ssc = c.square_symbol(P)
+    if ssc == 1:
+        return 1
+    if ssa == -1 and ssb == -1:
+        return 1
+
+    return -1
+
+
 
 def hilbert_conductor(a, b):
     """

sage/rings/number_field/number_field.py

@@ -7823 +7823 @@
     # field with lower precision will equal the lower-precision root!
     diffs = [(RC(ee(K.gen()))-old_root).abs() for ee in elist]
     return elist[min(izip(diffs,count()))[1]]
+
+def generalized_even_hilbert_symbol(a,b,P):
+    r"""
+    Returns the Hilbert symbol `(a,b)_P` for `a,b` elements of a number field `K` and `P` a prime ideal of `K`.
+    
+    INPUT:
+
+    - ``a,b`` --number field elements
+    - ``P`` --prime ideal
+    
+    OUTPUT:
+
+    - 1, -1, or 0
+    
+    ALGORITHM: This is Algorithm 6.6 from John Voight: Identifying the Matrix Ring. 
+        
+    EXAMPLES::
+
+        sage: from sage.rings.number_field.number_field import generalized_even_hilbert_symbol
+        sage: K.<a> = NumberField(x^5-23)
+        sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+        sage: b = a+5
+        sage: generalized_even_hilbert_symbol(a,b,P)
+        -1
+        sage: generalized_even_hilbert_symbol(a,2,P)
+        -1
+        sage: generalized_even_hilbert_symbol(a+5,2,P)
+        -1
+    """
+    
+    #step 0: checks
+    F = P.ring()
+    a,b = F(a),F(b)
+    if not P.divides(F.ideal(2)):
+        raise ValueError, 'P must be even'
+    
+    #step 1: assume a and b have been scaled
+    #step 2:
+    
+    x,y,z,w = _voight_alg_6_5(a,b,P)
+
+    Rp = F.residue_field(P)
+    pmap = Rp.lift_map()
+    R = F['T']
+    T = R.gen()
+    #the norm of i prime:
+    Niprime = (1-a*y**2-b*z**2+a*b*w**2)/4
+    #here f is the min poly of iprime
+    f = T**2-T+Niprime
+    for r in Rp:
+        if f(pmap(r))==0:
+            return 1
+    
+    bprime = (z**2)*(b**2)*a+(y**2)*(a**2)*b
+    if bprime.valuation(P) % 2 == 0:
+        return 1
+        
+    return -1   
+ 
+def _voight_alg_6_2(a,b,P):
+    r"""
+    This function returns a solution to the congruence 
+    `1-a*y^2-b*z^2 \equiv 0 (\mod P^{2*e})` where `e` is the ramification index of `P` over 2,
+    i.e. it returns a solution to an integral norm form (via a Hensel-like lifting).
+    This is a special case of the more general `x^2-a*y^2-b*z^2+a*b*w^2 \equiv 0 (\mod P^{2e}).`
+
+    INPUT:
+
+    - ``a,b`` --integers in the number field `F`
+    - ``P``  -- a prime ideal dividing 2 
+
+    OUTPUT:
+
+    a tuple `1,y,z,0` --integer elements of `F`
+
+    EXAMPLES::
+
+        sage: from sage.rings.number_field.number_field import _voight_alg_6_2
+        sage: K.<a> = NumberField(x^5-23)
+        sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+        sage: Q = K.ideal(-7*a^4 + 13*a^3 - 13*a^2 - 2*a + 50)
+        sage: b = -a+5
+        sage: _voight_alg_6_2(a,2,P)
+        (1, a^2, 1, 0)
+        sage: _voight_alg_6_2(b,2,P)
+        (1, a^2 + a + 1, a^3 + a^2, 0)
+        sage: _voight_alg_6_2(2,b,P)
+        (1, a^3 + a^2, a^2 + a + 1, 0)
+
+    """
+
+    #some initializations
+    F = P.number_field()
+    a = F(a)
+    b = F(b)
+
+    if not P.is_prime():
+        raise ValueError, "P must be prime"
+    if not P.divides(F.ideal(2)):
+        raise ValueError, "P must be even"
+    e = P.ramification_index()
+    k = P.residue_field()
+    E = 2*e
+    va = a.valuation(P)
+    vb = b.valuation(P)
+    switch = False
+    
+    if va != 0 and vb != 0:
+        raise ValueError, 'ord_P(a) or ord_P(b) must be zero'
+    if va == 0 and vb == 0:
+        raise ValueError, 'wrong case, use _voight_alg_6_5'
+    if vb == 0:
+        switch = True
+        a,b = b,a
+        va = a.valuation(P)
+        vb = b.valuation(P)
+
+    #step 1
+    f = P.residue_class_degree()
+    q = 2**f
+    pi = P.gens_reduced()[-1]
+ 
+
+    #step 2
+    y = a.solver_mod_p(P,1)
+    z = 0
+
+    #step 3
+
+    N = F(1-a*y**2-b*z**2).mod(F.ideal(4))
+
+    if N == 0:
+        return 1,y,z,0
+    t = N.valuation(P)
+
+    while t < E:
+        if E <= t:
+            break
+        if t%2==0:
+            ins1 = (N/(a*pi**t)).sqrt_mod_p(P)
+            y = y+ins1*pi**(t/2)
+            N = (1-a*y**2-b*z**2).mod(F.ideal(4))
+            t = N.valuation(P)
+
+        if E <= t:
+            break
+
+        if t%2==1:
+            ins2 = (N/(b*pi**(t-1))).sqrt_mod_p(P)
+            z = z+ins2*pi**((t/2).floor())
+            N = (1-a*y**2-b*z**2)
+            t = N.valuation(P)
+    
+        if E <= t:
+            break
+    #end of while
+    
+    if switch:
+        return 1,z,y,0
+    return 1,y,z,0    
+
+
+def _voight_alg_6_5(a,b,P):
+    r"""
+    Returns a solution to the congruence `1-a*y^2-b*z^2 +a*b*w^2 \equiv 0 (\mod P^{2e})`
+    where `e` is the ramification index of `P` over 2,
+    i.e. it returns a solution to an integral norm form (via a Hensel-like lifting).
+
+    INPUT:
+
+    - ``a,b`` --integers in the number field `F`
+    - ``P``  -- a prime ideal dividing 2 
+
+    OUTPUT:
+
+    a tuple `1,y,z,w` --integer elements of `F`
+
+    EXAMPLES::
+    
+        sage: from sage.rings.number_field.number_field import _voight_alg_6_5
+        sage: K.<a> = NumberField(x^5-23)
+        sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+        sage: b = -a+5
+        sage: _voight_alg_6_5(a,b,P)
+        (1, a^2, a^2 + a + 1, a^4 + a^3 + a^2)
+        sage: _voight_alg_6_5(a*b,b,P)
+        (1, a + 1, a^2 + a + 1, a^3 + 2*a^2 + 2*a + 1)
+
+
+
+    """
+
+    F = P.number_field()  
+    a,b = F(a),F(b)  
+    pi = P.uniformizer()
+    if not P.is_prime():
+        raise ValueError, 'P must be a prime'
+    if not P.divides(F.ideal(2)):
+        raise ValueError, 'P must lie over 2'
+        
+    va = a.valuation(P)
+    vb = b.valuation(P)
+    if va != 0 and vb != 0:
+        raise ValueError, 'ord_P(a) or ord_P(b) must be zero'
+    #fix these to reducing b and c:
+    if va == 0 and vb > 1:
+        b = b / pi**(2*(b.valuation(P)//2))
+        vb = b.valuation(P)
+    if vb == 0 and va > 1:
+        a = a / pi**(2*(a.valuation(P)//2))
+        va = a.valuation(P)
+    if va >1 or vb >1:
+        raise ValueError, 'the reduction did not work'
+        
+    Zp = P.residue_field()
+    e = P.ramification_index()
+    pi = P.gens_reduced()[0]    
+    E = 2*e
+
+    
+    #step 1
+    #first possibility
+    if va == 1 and vb == 0:
+        return _voight_alg_6_2(a,b,P)
+    #second possibility
+    if va == 0 and vb == 1:
+        return _voight_alg_6_2(a,b,P)
+    
+    #step 2
+    #check if we can solve x^2*c = 1 mod P^2 for c = a and c = b
+    #if not, returns solution = 0 
+    while True:
+        try:
+            a0 = a.solver_mod_p(P,e)
+            break
+        except ZeroDivisionError:
+            a0 = 0
+        break
+    while True:
+        try:
+            b0 = b.solver_mod_p(P,e)
+            break
+        except ZeroDivisionError:
+            b0 = 0
+        break
+    
+    if a0 and b0:
+        a0 = F(a).solver_mod_p(P,e)
+        b0 = F(b).solver_mod_p(P,e)
+        return 1,a0,b0,a0*b0
+        
+    #end step 2
+    
+    #step 3
+    #swap so that a0 is not invertible mod P^e
+    switch = False
+    if not b0:
+        a0,b0 = b0,a0
+        switch = True
+        
+    #now assuming that a0 is not invertible mod P^e
+    t = 1
+    aIn = a.solver_mod_p(P,1)
+    while aIn:
+        try:
+            t = t+1
+            aIn = a.solver_mod_p(P,t)
+            break
+        except ZeroDivisionError:
+            aIn = 0
+            break
+    
+
+    top = F(a)-F(a0**2)
+    at = top/(pi**t)
+    newb = -pi*at/b
+    x,y,z,w = _voight_alg_6_2(a,newb,P)
+    xnew = 1
+    
+    ynew = a.solver_mod_p(P,1)
+    zinv = z.inverse_mod(P)
+    znew = pi**((t/2).floor())*zinv*ynew
+    wnew = y*pi**((t/2).floor())*zinv*ynew
+    if switch:
+        return xnew,znew,ynew,wnew
+    return xnew,ynew,znew,wnew
+

sage/rings/number_field/number_field_element.pyx

@@ -2752 +2752 @@
             raise NotImplementedError, "inverse_mod is not implemented for non-integral elements"                
 
 
+
+    def sqrt_mod_p(self,P):
+        r"""
+       
+        Returns the square root of self mod the prime ideal `\mathfrak{p}`.
+
+        INPUT:
+
+        -  ``P`` -- is a prime ideal be an ideal of self.parent(). A ValueError
+           is raised if `\mathfrak{p}` is not prime or self is not a square mod `\mathfrak{p}`.
+
+        OUTPUT:
+
+        The square root of self mod `\mathfrak{p}`
+
+        .. NOTE:: 
+
+            This is not yet implemented for non-integral elements.
+
+        EXAMPLES::
+
+            sage: K.<a> = NumberField(x^5-23)
+            sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+            sage: Q = K.ideal(-7*a^4 + 13*a^3 - 13*a^2 - 2*a + 50)
+            sage: b = -a+5
+            sage: a.sqrt_mod_p(P)
+            a^3
+            sage: b.sqrt_mod_p(P)
+            a^3 + 1
+
+        ::
+
+            sage: k.<a> = NumberField(x^3 + 11)
+            sage: d = a + 3
+            sage: d.sqrt_mod_p(k.ideal(-a-2))
+            1
+
+        """
+        k = P.residue_field()
+        f = k.lift_map()
+        if not P.is_prime():
+            raise ValueError, 'P must be a prime ideal.'
+        if not k(self).is_square():
+            raise ValueError, 'must be a square mod P'
+        return f(k(self).sqrt())
+
+    def sqrt_mod(self,P,n, verb = False):
+        r"""
+       
+        Returns the square root of self mod the prime ideal `\mathfrak{p}^n`.
+        
+        INPUT:
+
+        -  ``P`` -- is a prime ideal be an ideal of self.parent(). A ValueError
+           is raised if `\mathfrak{p}` is not prime or self is not a square mod `\mathfrak{p}`.
+
+        OUTPUT:
+
+        The square root of self mod `\mathfrak{p}^n`
+
+        .. NOTE:: 
+
+            This is not yet implemented for non-integral elements.
+
+        EXAMPLES::
+
+            sage: K.<a> = NumberField(x^5-23)
+            sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+            sage: Q = K.ideal(-7*a^4 + 13*a^3 - 13*a^2 - 2*a + 50)
+            sage: b = -a+5
+            sage: a.sqrt_mod(P,1)
+            a^3
+            sage: b.sqrt_mod(P,1)
+            a^3 + 1
+            sage: ((a^2).sqrt_mod(P,2)^2 - a^2) in P^2
+            True
+
+        ::
+
+            sage: k.<a> = NumberField(x^3 + 11)
+            sage: d = a + 3
+            sage: d.sqrt_mod(k.ideal(-a-2),1)
+            1
+            sage: d.sqrt_mod(k.ideal(-a-2),2)
+            -a^2
+
+        """
+        if not P.is_prime():
+            raise ValueError, 'P must be a prime ideal.'
+        if n == 1:
+            return self.sqrt_mod_p(P)
+        Pn = P**n
+        I = Pn.residues()
+        for i in I:
+            if (i**2 - self) in Pn:
+                return i
+        raise ValueError, 'self is not a square root mod P^n'
+        
+    def solver_mod_p(self,P,n):
+        r"""
+        Returns `x` such that `x^2 a \equiv 1 \text{ } (\text{mod }{ \mathfrak{p}^n})` where `a =` self.
+        
+        INPUT:
+
+        -  ``P`` - is a prime ideal be an ideal of self.parent(). A ValueError
+           is raised if `\mathfrak{p}` is not prime or self has no solution to the congruence 
+           mod `\mathfrak{p}^n.`
+
+        - ``n`` - a positive integer
+
+        OUTPUT:
+
+        The solution to the above congruence.
+        
+        EXAMPLES::
+
+            sage: K.<a> = NumberField(x^5-23)
+            sage: P = K.ideal(-1105*a^4 + 1541*a^3 - 795*a^2 - 2993*a + 11853)
+            sage: Q = K.ideal(-7*a^4 + 13*a^3 - 13*a^2 - 2*a + 50)
+            sage: b = -a+5
+            sage: a.solver_mod_p(P,1)
+            a^2
+            sage: (b^2).solver_mod_p(P,2)
+            -a^3 - a
+        
+        """
+        
+        Ppow = P**n
+        a = self.inverse_mod(Ppow)
+        x = a.sqrt_mod(P,n)
+        return x
+
+    def square_symbol(self,P):
+        r"""
+        Returns the square symbol from John Voight's: Identifying the Matrix Ring.  
+        Let `K` be a number field, `a \in K` and `\mathfrak{p}` a prime ideal of `K` then:
+        `\{ \frac{a}{\mathfrak{p}}\} =  1` if `a \in K_{\mathfrak{p}}^{* 2}`, 
+        `\{ \frac{a}{\mathfrak{p}}\} = -1` if `a \not \in K_{\mathfrak{p}}^{* 2}` and 
+        `\text{ord}_{\mathfrak{p}}(a)` is even, 
+        and 
+        `\{ \frac{a}{\mathfrak{p}}\} =  0` if `a \not \in K_{\mathfrak{p}}^{* 2}` and 
+        `\text{ord}_{\mathfrak{p}}(a)` is odd.
+        
+        INPUT:
+
+        - ``P`` --prime ideal
+        
+        OUTPUT:
+
+        - 1,-1, or 0 
+        
+        ALGORITHM: The algorithm follows the section on Computing Local Hilbert Symbols from 
+        John Voight: Identifying the Matrix Ring.
+
+        EXAMPLE::
+
+            sage: K.<a> = NumberField(x^5-23)
+            sage: Q = K.ideal(-7*a^4 + 13*a^3 - 13*a^2 - 2*a + 50)
+            sage: b = -a+5
+            sage: a.square_symbol(Q)
+            1
+            sage: b.square_symbol(Q)
+            -1
+
+        """
+        from sage.rings.arith import generalized_legendre_symbol
+        if not P.is_prime():
+            raise ValueError, "P must be a prime"
+        K = P.ring()
+        if P.divides(2):
+            raise ValueError, "P must be odd"
+        pi = P.uniformizer()
+        a = K(self)
+        e = a.valuation(P)
+        if e%2 == 1:
+            square_symbol_a = 0
+        else:
+            pi = P.uniformizer()
+            a0 = a*pi**(-e/2)
+            square_symbol_a = generalized_legendre_symbol(a0,P)
+            
+        return square_symbol_a
+
+
+
 cdef class NumberFieldElement_absolute(NumberFieldElement):
 
     def _pari_(self, var='x'):

sage/rings/number_field/number_field_ideal.py

@@ -1167 +1167 @@
         """
         return self.number_field().galois_group().artin_symbol(self)
 
+    def uniformizer(self):
+        r"""
+        Returns an element `\pi` of `\mathfrak{p} =` self such that `\text{ord}(\pi)_{\mathfrak{p}} = 1`. 
+        This wraps K.uniformizer(self).
+        
+        
+        OUTPUT:  
+        `\pi` -- a uniformizer of `\mathfrak{p}`
+        
+        
+            
+        EXAMPLES::
+            
+            sage: K.<a> = NumberField(x^2-10)
+            sage: P = K.ideal(3).factor()[0][0]
+            sage: P.uniformizer()
+            -2*a + 1
+            
+        ::
+            
+            sage: K.<a> = NumberField(x^5-23)
+            sage: P = K.ideal(3).factor()[0][0]
+            sage: P.uniformizer()
+            a^4 - a^3 + a^2 - a + 1
+            
+        ::
+            
+            sage: K.<a> = NumberField(x^2+1489)
+            sage: P = K.ideal(3).factor()[0][0]
+            sage: P.uniformizer()
+            3
+                    
+        """
+        #this just raps K.uniformizer(self)
+        K = self.ring()
+        return K.uniformizer(self)
+
+
+
 def basis_to_module(B, K):
     r"""
     Given a basis `B` of elements for a `\ZZ`-submodule of a number