Ticket #727: trac_727_more_conic_files3.patch

File trac_727_more_conic_files3.patch, 7.2 KB (added by Marco Streng, 12 years ago)

needs two previous tickets trac_727_more... and trac_9334-hilbert.patch

  • sage/schemes/plane_conics/conic_solve.py

    # HG changeset patch
    # User Marco Streng <marco.streng@gmail.com>
    # Date 1278581937 -7200
    # Node ID 510cd246e899886d8808693f7c8f257c71e13451
    # Parent  bd313996cb3dd157598db0057df8fccffe9f2aa1
    Trac 727: replace unfinished hilbert symbols by 9334
    
    diff -r bd313996cb3d -r 510cd246e899 sage/schemes/plane_conics/conic_solve.py
    a b  
    385385            return False
    386386        return True
    387387    if is_NumberField(B):
    388         from hilbert import hilbert_symbol_nf
    389         if hilbert_symbol_nf(a, b, p) == -1:
     388        from sage.rings.arith import generalized_hilbert_symbol
     389        if generalized_hilbert_symbol(a, b, p) == -1:
    390390            if self._local_obstruction == None:
    391391                self._local_obstruction = p
    392392            return False
  • deleted file sage/schemes/plane_conics/hilbert.py

    diff -r bd313996cb3d -r 510cd246e899 sage/schemes/plane_conics/hilbert.py
    + -  
    1 r"""
    2 Compute Hilbert symbols over number fields.
    3 
    4 AUTHORS:
    5 
    6 - Marco Streng (2009-08-07)
    7 
    8 """
    9 #*****************************************************************************
    10 #       Copyright (C) 2009 Marco Streng <marco.streng@gmail.com>
    11 #
    12 #  Distributed under the terms of the GNU General Public License (GPL)
    13 #
    14 #    This code is distributed in the hope that it will be useful,
    15 #    but WITHOUT ANY WARRANTY; without even the implied warranty of
    16 #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    17 #    General Public License for more details.
    18 #
    19 #  The full text of the GPL is available at:
    20 #
    21 #                  http://www.gnu.org/licenses/
    22 #*****************************************************************************
    23 
    24 from sage.rings.all import (ZZ, is_NumberFieldElement, is_RingHomomorphism,
    25                             is_RealField, is_ComplexField, is_NumberFieldIdeal)
    26 
    27 
    28 def hilbert_symbol_nf(a, b, p, algorithm="default"):
    29     r"""
    30     Returns the Hilbert symbol `(a, b / p)`, that is,
    31     returns 1 if `ax^2 + by^2` represents a nonzero
    32     `P`-adic square, otherwise returns `-1`. If either a or b
    33     is 0, returns 0.
    34    
    35     INPUT:
    36    
    37     -  ``p``          - prime element, prime ideal, or real or complex
    38                         embedding of a number field K
    39    
    40     -  ``a, b``       - elements of K
    41    
    42     -  ``algorithm``  - string
    43    
    44        -  ``'direct'``  - (default) use a Python implementation
    45 
    46        -  ``'magma'``   - use MAGMA (optional)
    47 
    48        -  ``'all'``     - use both algorithms and check that the
    49                           results agree, then return the common answer
    50 
    51        -  ``'default'`` - same as 'direct'
    52    
    53 
    54     OUTPUT: integer (0, -1, or 1)
    55    
    56     EXAMPLES::
    57    
    58         sage: from sage.schemes.plane_conics.hilbert import hilbert_symbol_nf
    59         sage: _.<x> = QQ[]
    60         sage: K.<a> = NumberField(x^3-2)
    61         sage: O = K.maximal_order()
    62         sage: p = (2*O).factor()[0][0]
    63         sage: q = (7*O)
    64         sage: (r, c) = K.places()
    65         sage: [[hilbert_symbol_nf(t, -2, s) \
    66                 for s in [p,q,r,c]] for t in [7, -7, 0]]
    67         [[-1, -1, 1, 1], [1, -1, -1, 1], [0, 0, 0, 0]]
    68         sage: [[hilbert_symbol_nf(t, -2, s, 'all') for s in [p,q,r,c]] for t in [7, -7, 0]] # optional: uses Magma, and long time (3 seconds)
    69         [[-1, -1, 1, 1], [1, -1, -1, 1], [0, 0, 0, 0]]   
    70     """
    71     one = ZZ(1)
    72     zero = ZZ(0)
    73     if is_NumberFieldElement(p):
    74         K = p.parent()
    75         O = K.maximal_order()
    76         p = p*O
    77     elif is_RingHomomorphism(p):
    78         if is_RealField(p.codomain()):
    79             if a == 0 or b == 0:
    80                 return zero
    81             if p(a) > 0 or p(b) > 0:
    82                 return one
    83             return -one
    84         if is_ComplexField(p.codomain()):
    85             if all([im.is_real() for im in p.im_gens()]):
    86                 raise TypeError, "Could not determine whether the place p " \
    87                                  "(= %s) is real or complex in " \
    88                                  "hilbert_symbol_nf" % p
    89             if a == 0 or b == 0:
    90                 return zero
    91             return one
    92     if (not is_NumberFieldIdeal(p)) or (not p.is_prime()):
    93         raise TypeError, "p (= %s) must be a prime of a number field in " \
    94                          "hilbert_symbol_nf" % p
    95     K = p.number_field()
    96     if not a in K and b in K:
    97         raise TypeError, "a (=%s) and b (=%s) must be elements of the number " \
    98                          "field of p (=%s)" % (a, b, p)
    99     a = K(a)
    100     b = K(b)
    101     if a == 0 or b == 0:
    102         return zero
    103     pi = [g for g in p.gens() if g.valuation(p) == 1][0]
    104     if algorithm == 'magma':
    105         from sage.interfaces.magma import magma
    106         if magma(a).HilbertSymbol(b, pi*magma(K).MaximalOrder()) == 1:
    107             return one
    108         return -one
    109     elif algorithm == 'all':
    110         magmaout = hilbert_symbol_nf(a, b, p, algorithm = 'magma')
    111         sageout = hilbert_symbol_nf(a, b, p, algorithm = 'direct')
    112         if magmaout != sageout:
    113             raise RuntimeError, "There is a bug in hilbert_symbol_nf: " \
    114                                 "different methods for computing the Hilbert " \
    115                                 "symbol (%s, %s / %s) disagree" % (a, b, p)
    116         return sageout
    117     elif algorithm != 'default' and algorithm != 'direct':
    118         raise ValueError, "Unknown algorithm %s in hilbert_symbol_nf" % \
    119                           algorithm
    120     av = a.valuation(p)
    121     bv = b.valuation(p)
    122     a *= pi**(-2*(av/2).floor())
    123     b *= pi**(-2*(bv/2).floor())
    124     coefficients_are_units = av%2==0 and bv%2==0
    125     if av%2==1 and bv%2==1:
    126         a = -a*b/pi**2
    127     elif av%2==1:
    128         t = b; b = a; a = t
    129     if p.smallest_integer() != 2:
    130         if coefficients_are_units:
    131             return one
    132         k = p.residue_field()
    133         if k(a).is_square():
    134             return one
    135         return -one
    136 
    137 
    138     res = list(p.residues())
    139    
    140     e = p.ramification_index()
    141 
    142     potential_solutions = [(a, b, -1, 0, 0, 0), \
    143                            (-1, a, b, 0, 0, 0), \
    144                            (b, -1, a, 0, 0, 0)]
    145     # If there is a solution, then it can be scaled
    146     # to (x:y:1), (y:1:x) or (1:x:y) with integral x, y.
    147     # This is thus a solution to fx^2+gy^2+h=0 with
    148     # (f, g, h) in {(a, b, -1), (-1, a, b), (b, -1, a)}.
    149     # For every solution to this equation, there is an
    150     # element (f, g, h, X, Y, k) in the set potential_solutions
    151     # such that X = x mod P^k and Y = y mod P^k hold and (f,g,h,X,Y)
    152     # is a solution modulo P^(k+e).
    153 
    154     while True:
    155         if len(potential_solutions) == 0:
    156             return -one
    157         (f, g, h, X, Y, k) = potential_solutions.pop()
    158         for u in res:
    159             Xn = X + u * pi^k
    160             for v in res:
    161                 Yn = Y + v * pi^k
    162                 if (f*Xn**2 + g*Yn**2 + h).valuation(p) < k+e:
    163                     if k == e:
    164                         # Found a solution mod P^(3+2*e), hence
    165                         # a solution exists by Hensel's lemma.
    166                         return one
    167                     potential_solutions.push((f, g, h, Xn, Yn, k+1))
    168 
    169    
    170