Ticket #727: trac_727_more_conic_files.patch

File trac_727_more_conic_files.patch, 86.4 KB (added by Marco Streng, 12 years ago)
  • new file sage/quadratic_forms/qfsolve.py

    # HG changeset patch
    # User Marco Streng <marco.streng@gmail.com>
    # Date 1278425230 -7200
    # Node ID 3fcab21a698fa0d667d07a8d380310c8e852c2cf
    # Parent  8dec8b43ccca5f104b1e280cb33c8f4c2c1b8f85
    Trac 727: Conic class using Simon's code for finding points over QQ and using
    Pari to find points over number fields. Bad version of Hilbert symbols over
    number fields.
    
    diff -r 8dec8b43ccca -r 3fcab21a698f sage/quadratic_forms/qfsolve.py
    - +  
     1"""
     2qfsolve: Programme de resolution des equations quadratiques.
     3
     4Interface to the GP quadratic forms code of Denis Simon.
     5
     6AUTHOR:
     7    Denis Simon (GP code)
     8    Nick Alexander (Sage interface)
     9"""
     10
     11#*****************************************************************************
     12#       Copyright (C) 2008 William Stein <wstein@gmail.com>
     13#
     14#  Distributed under the terms of the GNU General Public License (GPL)
     15#
     16#    This code is distributed in the hope that it will be useful,
     17#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     18#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     19#    General Public License for more details.
     20#
     21#  The full text of the GPL is available at:
     22#
     23#                  http://www.gnu.org/licenses/
     24#*****************************************************************************
     25
     26from sage.interfaces.gp import Gp
     27from sage.rings.all import ZZ, QQ
     28from sage.libs.pari.gen import pari
     29
     30_gp_for_simon_interpreter = None    # Global GP interpreter for Denis Simon's code
     31def _gp_for_simon():
     32    global _gp_for_simon_interpreter
     33    if _gp_for_simon_interpreter is None:
     34        _gp_for_simon_interpreter = Gp(script_subdirectory='simon')
     35        _gp_for_simon_interpreter.read("qfsolve.gp")
     36    return _gp_for_simon_interpreter
     37
     38# \\ - Qfsolve(G,factD): pour resoudre l'equation quadratique X^t*G*X = 0
     39# \\ G doit etre une matrice symetrique n*n, a coefficients dans Z.
     40# \\ S'il n'existe pas de solution, la reponse est un entier
     41# \\ indiquant un corps local dans lequel aucune solution n'existe
     42# \\ (-1 pour les reels, p pour Q_p).
     43# \\ Si on connait la factorisation de -abs(2*matdet(G)),
     44# \\ on peut la passer par le parametre factD pour gagner du temps.
     45# \\
     46# \\ - Qfparam(G,sol,fl): pour parametrer les solutions de la forme
     47# \\ quadratique ternaire G, en utilisant la solution particuliere sol.
     48# \\ si fl>0, la 'fl'eme forme quadratique est reduite.
     49
     50def Qfsolve(G, factD=None):
     51    r"""Find a solution to (x,y,z)^t*G*(x,y,z) = 0.
     52
     53    If successful, returns a triple of rational numbers (x0, y0, z0).
     54    If unsuccessful, returns -1 if no solutions exists over the reals and a
     55    prime p if no solution exists over the p-adic field Q_p.
     56
     57    EXAMPLES:
     58        sage: from sage.schemes.plane_curves.qfsolve import Qfsolve
     59        sage: R.<x,y,z> = QQ[]
     60        sage: M = Curve(-12*y^2 - 24*x*z - z^2).symmetric_matrix(); M
     61        [  0   0 -12]
     62        [  0 -12   0]
     63        [-12   0  -1]
     64        sage: sol = Qfsolve(M); sol
     65        (1, 0, 0)
     66        sage: sol[0].parent() is QQ
     67        True
     68       
     69        sage: C = Curve(x^2 + y^2 + z^2)
     70        sage: ret = Qfsolve(C.symmetric_matrix()); ret
     71        -1
     72        sage: ret.parent() is ZZ
     73        True
     74
     75        sage: C = Curve(x^2 + y^2 - 7*z^2)
     76        sage: Qfsolve(C.symmetric_matrix())
     77        7
     78    """
     79    gp = _gp_for_simon()
     80    if factD is not None:
     81        raise NotImplementedError
     82    ret = pari(gp('Qfsolve(%s)' % G._pari_()))
     83    if ret.type() == 't_COL':
     84        return (QQ(ret[0]), QQ(ret[1]), QQ(ret[2]))
     85    return ZZ(ret)
     86
     87def Qfparam(G, sol):
     88    r"""
     89
     90    EXAMPLES:
     91        sage: from sage.schemes.plane_curves.qfsolve import Qfsolve, Qfparam
     92        sage: R.<x,y,z> = QQ[]
     93        sage: M = Curve(-12*y^2 - 24*x*z - z^2).symmetric_matrix(); M
     94        [  0   0 -12]
     95        [  0 -12   0]
     96        [-12   0  -1]
     97        sage: sol = Qfsolve(M);
     98        sage: ret = Qfparam(M, sol); ret
     99        (-t^2 - 12, 24*t, 24*t^2)
     100        sage: ret[0].parent() is QQ['t']
     101        True
     102    """
     103    gp = _gp_for_simon()
     104    R = QQ['t']
     105    str = 'Qfparam(%s, (%s)~)*[t^2,t,1]~' % (G._pari_(), pari(gp(sol))._pari_())
     106    ret = pari(gp(str))
     107    return (R(ret[0]), R(ret[1]), R(ret[2]))
     108
     109qfsolve = Qfsolve
     110qfparam = Qfparam
  • sage/schemes/all.py

    diff -r 8dec8b43ccca -r 3fcab21a698f sage/schemes/all.py
    a b  
    2626
    2727from plane_curves.all import *
    2828
     29from plane_conics.all import *
     30
    2931from elliptic_curves.all import *
    3032
    3133from plane_quartics.all import *
  • new file sage/schemes/plane_conics/all.py

    diff -r 8dec8b43ccca -r 3fcab21a698f sage/schemes/plane_conics/all.py
    - +  
     1"""
     2Plane conics
     3"""
     4
     5#*****************************************************************************
     6#
     7#   SAGE: System for Algebra and Geometry Experimentation   
     8#
     9#       Copyright (C) 2005 William Stein <was@math.harvard.edu>
     10#
     11#  Distributed under the terms of the GNU General Public License (GPL)
     12#
     13#    This code is distributed in the hope that it will be useful,
     14#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     15#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     16#    General Public License for more details.
     17#
     18#  The full text of the GPL is available at:
     19#
     20#                  http://www.gnu.org/licenses/
     21#*****************************************************************************
     22
     23from constructor import Conic
     24
     25
  • new file sage/schemes/plane_conics/conic_solve.py

    diff -r 8dec8b43ccca -r 3fcab21a698f sage/schemes/plane_conics/conic_solve.py
    - +  
     1r"""
     2Find rational points on plane conics over fields.
     3
     4AUTHORS:
     5
     6- Marco Streng (2009-08-07)
     7
     8- Nick Alexander (2008-01-08)
     9
     10"""
     11#*****************************************************************************
     12#       Copyright (C) 2008 Nick Alexander <ncalexander@gmail.com>
     13#       Copyright (C) 2009 Marco Streng <marco.streng@gmail.com>
     14#
     15#  Distributed under the terms of the GNU General Public License (GPL)
     16#
     17#    This code is distributed in the hope that it will be useful,
     18#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     19#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     20#    General Public License for more details.
     21#
     22#  The full text of the GPL is available at:
     23#
     24#                  http://www.gnu.org/licenses/
     25#*****************************************************************************
     26
     27from sage.rings.all import (PolynomialRing, ZZ, is_FiniteField,
     28                            is_ComplexField, is_RealField,
     29                            is_pAdicField, is_RationalField,
     30                            is_NumberField, RR, is_RingHomomorphism)
     31from sage.modules.free_module_element import vector
     32from sage.structure.sequence import Sequence
     33from sage.schemes.generic.projective_space import ProjectiveSpace
     34from sage.matrix.constructor import Matrix
     35
     36
     37
     38
     39def has_rational_point(self, point = False, obstruction = False,
     40                       algorithm = 'default', read_cache = True):
     41    r"""
     42    Returns True if and only if self has a point defined over its base field B.
     43
     44    If point or obstruction is True, then returns a second output S:
     45    - if point is True and self has a rational point,
     46      then S is a rational point;
     47    - if obstruction is True, self has no rational point, and the base
     48      ring is a global field, then S is a prime or infinite place
     49      of the base ring such that no rational point exists
     50      over the localization at S.
     51   
     52    Points and obstructions are cached, even if they are not returned as
     53    output, but only if the algorithm finds them. If read_cache is True,
     54    then cached information is used for the output if available.
     55   
     56    ALGORITHM:
     57       
     58        If the base field is finite, try random y-coordinates.
     59       
     60        If the base field is a number field, the parameter ``algorithm``
     61        specifies the algorithm to be used:
     62       
     63            ``qfsolve``   - Use Denis Simon's Qfsolve (only over `\QQ)
     64            ``rnfisnorm`` - Use PARI's rnfisnorm (cannot be combined with
     65                            obstruction = True)
     66            ``local``     - Check if a local solution exists for all primes
     67                            and infinite places of K (cannot be combined with
     68                            point = True)
     69            ``default``   - Use ``qfsolve`` over `\QQ`. Use ``local`` over
     70                            other number fields. If the output is True and
     71                            point is True, then use ``rnfisnorm`` to find
     72                            the point.
     73            ``all``       - Use all applicable algorithms, check that the
     74                            outputs agree, and return the common answer.
     75       
     76    EXAMPLES:
     77
     78    Examples over QQ ::
     79
     80        sage: C = Conic(QQ, [1, 2, -3])
     81        sage: C.has_rational_point(point = True, obstruction = True)
     82        (True, (-1 : -1 : 1))
     83        sage: D = Conic(QQ, [1, 3, -5])
     84        sage: D.has_rational_point(point = True, obstruction = True, \
     85                                   algorithm = 'all')
     86        (False, 3)
     87
     88    The following would not terminate quickly with algorithm = 'rnfisnorm' ::
     89
     90        sage: C = Conic(QQ, [1, 113922743, -310146482690273725409])
     91        sage: C.has_rational_point(point = True)
     92        (True, (-76842858034579/5424 : -5316144401/5424 : 1))
     93        sage: C.has_rational_point(algorithm = 'local', read_cache = False)
     94        True
     95
     96    Examples over number fields ::
     97       
     98        sage: K.<i> = QuadraticField(-1)
     99        sage: C = Conic(K, [1, 3, -5])
     100        sage: C.has_rational_point(point = True, obstruction = True, \
     101                                   algorithm = 'all')
     102        (False, Fractional ideal (-i - 2))
     103
     104        sage: _.<x> = QQ[]
     105        sage: L.<b> = NumberField(x^3-5)
     106        sage: C = Conic(L, [1, 2, -3])
     107        sage: C.has_rational_point(point = True, obstruction = True, \
     108                                   algorithm = 'all') # long time (1 second)
     109        (True, (5/3 : -1/3 : 1))
     110   
     111    Examples over finite fields ::
     112   
     113        sage: Conic(FiniteField(37), [1, 2, 3, 4, 5, 6]).has_rational_point()
     114        True
     115       
     116        sage: C = Conic(FiniteField(2), [1, 1, 1, 1, 1, 0]); C
     117        Projective Conic Curve over Finite Field of size 2 defined by x^2 + x*y + y^2 + x*z + y*z
     118        sage: C.has_rational_point(point = True)  # output is random
     119        (True, (0 : 0 : 1))
     120       
     121        sage: p = next_prime(10^50)
     122        sage: F = FiniteField(p)
     123        sage: C = Conic(F, [1, 2, 3]); C
     124        Projective Conic Curve over Finite Field of size 100000000000000000000000000000000000000000000000151 defined by x^2 + 2*y^2 + 3*z^2
     125        sage: C.has_rational_point(point = True)  # output is random
     126        (True,
     127         (14971942941468509742682168602989039212496867586852 : 75235465708017792892762202088174741054630437326388 : 1)
     128 
     129        sage: F.<a> = FiniteField(7^20)
     130        sage: C = Conic([1, a, -5]); C
     131        Projective Conic Curve over Finite Field in a of size 7^20 defined by x^2 + (a)*y^2 + 2*z^2
     132        sage: C.has_rational_point(point = True)  # output is random
     133        (True,
     134         (a^18 + 2*a^17 + 4*a^16 + 6*a^13 + a^12 + 6*a^11 + 3*a^10 + 4*a^9 + 2*a^8 + 4*a^7 + a^6 + 4*a^4 + 6*a^2 + 3*a + 6 : 5*a^19 + 5*a^18 + 5*a^17 + a^16 + 2*a^15 + 3*a^14 + 4*a^13 + 5*a^12 + a^11 + 3*a^10 + 2*a^8 + 3*a^7 + 4*a^6 + 4*a^5 + 6*a^3 + 5*a^2 + 2*a + 4 : 1))
     135         
     136    Examples over real and complex fields ::
     137
     138        sage: Conic(RR, [1, 1, 1]).has_rational_point()
     139        False
     140        sage: Conic(CC, [1, 1, 1]).has_rational_point()
     141        True
     142       
     143        sage: Conic(RR, [1, 2, -3]).has_rational_point(point = True)
     144        (True, (1.73205080756888 : 0.000000000000000 : 1.00000000000000))
     145    """
     146    if read_cache:
     147        if self._rational_point is not None:
     148            if point or obstruction:
     149                return True, self._rational_point
     150            else:
     151                return True
     152        if self._local_obstruction is not None:
     153            if point or obstruction:
     154                return False, self._local_obstruction
     155            else:
     156                return False
     157        if (not point) and self._finite_obstructions == [] and \
     158           self._infinite_obstructions == []:
     159            return True
     160    B = self.base_ring()
     161    if algorithm == 'default':
     162        if is_RationalField(B):
     163            algorithm = 'qfsolve'
     164        elif is_NumberField(B):
     165            ret = self.has_rational_point(point = False, obstruction = True,
     166                                          algorithm = 'local')
     167            if ret[0]:
     168                if point:
     169                    ret = self.has_rational_point(point = True,
     170                                                  obstruction = False,
     171                                                  algorithm = 'rnfisnorm')
     172                    if not ret[0]:
     173                        raise RuntimeError, "Outputs of algorithms in " \
     174                                            "has_rational_point disagree for " \
     175                                            "conic %s" % self
     176                    return ret
     177                if obstruction:
     178                    return True, None
     179                return True
     180            if point or obstruction:
     181                return ret
     182            return False
     183    if algorithm == 'all':
     184        ret = []
     185        if is_RationalField(B):
     186            ret.append(self.has_rational_point(point = True,
     187                                               obstruction = True,
     188                                               algorithm = 'qfsolve',
     189                                               read_cache = False))
     190        if is_NumberField(B):
     191            ret.append(self.has_rational_point(obstruction = True,
     192                                               algorithm = 'local',
     193                                               read_cache = False))
     194            ret.append(self.has_rational_point(point = True,
     195                                               algorithm = 'rnfisnorm',
     196                                               read_cache = False))
     197        if all([r[0] for r in ret]):
     198            if point or obstruction:
     199                return [a for a in ret if a[1] != None][0]
     200            return True
     201        if all([not r[0] for r in ret]):
     202            if point or obstruction:
     203                return [a for a in ret if a[1] != None][0]
     204            return False
     205        raise RuntimeError, "Bug in has_rational_point: different " \
     206                            "algorithms yield different outputs for conic " \
     207                            "self (=%s)" % self
     208    if algorithm == 'qfsolve':
     209        if not is_RationalField(B):
     210            raise TypeError, "Algorithm qfsolve in has_rational_point only " \
     211                             "for conics over QQ, not over %s" % B
     212        from sage.quadratic_forms.qfsolve import Qfsolve
     213        from sage.rings.arith import lcm
     214        M = self.symmetric_matrix()
     215        M *= lcm([ t.denominator() for t in M.list() ])
     216        pt = Qfsolve(M)
     217        if pt in ZZ:
     218            if pt == -1:
     219                pt = B.embeddings(RR)[0]
     220            if self._local_obstruction == None:
     221                self._local_obstruction = pt
     222            if point or obstruction:
     223                return False, pt
     224            return False
     225        pt = self.point([pt[0], pt[1], pt[2]])
     226        if point or obstruction:
     227            return True, pt
     228        return True
     229    if algorithm == 'local':
     230        if point:
     231            raise ValueError, "Algorithm local cannot be combined with " \
     232                              "point = True in has_rational_point"
     233        obs = self.local_obstructions(infinite = True, finite = False,
     234                                      read_cache = read_cache)
     235        if obs != []:
     236            if obstruction:
     237                return False, obs[0]
     238            return False
     239        obs = self.local_obstructions(read_cache = read_cache)
     240        if obs == []:
     241            if obstruction:
     242                return True, None
     243            return True
     244        if obstruction:
     245            return False, obs[0]
     246        return False
     247    if algorithm == 'rnfisnorm':
     248        if obstruction:
     249            raise ValueError, "Algorithm rnfisnorm cannot be combined with " \
     250                              "obstruction = True in has_rational_point"
     251        D, T = self.diagonal_matrix()
     252        abc = [D[0,0], D[1,1], D[2,2]]
     253        for j in range(3):
     254            if abc[j] == 0:
     255                pt = self.point(T*vector({2:0,j:1}))
     256                if point or obstruction:
     257                    return True, pt
     258                return True
     259        if (-abc[1]/abc[0]).is_square():
     260            pt = self.point(T*vector([1, (-abc[1]/abc[0]).sqrt(), 0]))
     261            if point or obstruction:
     262                return True, pt
     263            return True
     264        if is_RationalField(B):
     265            K = B
     266            [KtoB, BtoK] = [K.hom(K) for _ in range(2)]
     267        elif is_NumberField(B):
     268            K = B.absolute_field('Y')
     269            [KtoB, BtoK] = K.structure()
     270        else:
     271            raise TypeError, "Algorithm rnfisnorm is only for QQ and number " \
     272                             "fields, not for %s" % B
     273        from rnfisnorm import rnfisnorm
     274        isnorm = rnfisnorm(BtoK(-abc[1]/abc[0]), BtoK(-abc[2]/abc[0]))
     275        if isnorm[0]:
     276            pt = self.point(T*vector([KtoB(isnorm[1][0]),
     277                                      KtoB(isnorm[1][1]), 1]))
     278            if point:
     279                return True, pt
     280            return True
     281        if point:
     282            return False, None
     283        return False       
     284    if is_FiniteField(B):
     285        if point:
     286            s, pt = self.has_singular_point(point = True)
     287            if s:
     288                return True, pt
     289            while True:
     290                x = B.random_element()
     291                Y = PolynomialRing(B,'Y').gen()
     292                r = self.defining_polynomial()([x,Y,1]).roots()
     293                if len(r) > 0:
     294                    return True, self.point([x,r[0][0],B(1)])
     295        if obstruction:
     296            return True, None
     297        return True
     298    if is_ComplexField(B):
     299        if point:
     300            [_,_,_,d,e,f] = self._coefficients
     301            if d == 0:
     302                return True, self.point([0,1,0])
     303            return True, self.point([0, ((e**2-4*d*f).sqrt()-e)/(2*d), 1],
     304                                    check = False)
     305        if obstruction:
     306            return True, None
     307        return True
     308    if is_RealField(B):
     309        D, T = self.diagonal_matrix()
     310        [a, b, c] = [D[0,0], D[1,1], D[2,2]]
     311        if a == 0:
     312            ret = True, self.point(T*vector([1,0,0]), check = False)
     313        elif a*c <= 0:
     314            ret = True, self.point(T*vector([(-c/a).sqrt(),0,1]),
     315                                   check = False)
     316        elif b == 0:
     317            ret = True, self.point(T*vector([0,1,0]), check = False)
     318        elif b*c <= 0:
     319            ret = True, self.point(T*vector([0,(-c/b).sqrt(),0,1]),
     320                                   check = False)
     321        else:
     322            ret = False, None
     323        if point or obstruction:
     324            return ret
     325        return ret[0]
     326    raise NotImplementedError, "has_rational_point not implemented for " \
     327                               "conics over base field %s" % B
     328
     329def is_locally_solvable(self, p):
     330    r"""
     331    Returns True if and only if self has a solution over the completion
     332    of the base field `B` of self at `p`. Here `p` is a finite prime
     333    or infinite place of `B`.
     334
     335    If `B` is `QQ`, then `p=-1` can be used to denote the infinite place.
     336
     337    EXAMPLES:
     338
     339    An example over QQ ::
     340
     341        sage: C = Conic(QQ, [1,2,3])
     342        sage: C.is_locally_solvable(-1)
     343        False
     344        sage: C.is_locally_solvable(2)
     345        False
     346        sage: C.is_locally_solvable(3)
     347        True
     348
     349    Example over a number field ::
     350
     351        sage: _.<x> = QQ[]
     352        sage: K.<a> = NumberField(x^3 + 5)
     353        sage: C = Conic(K, [1, 2, 3 - a])
     354        sage: [p1, p2] = K.places()
     355        sage: C.is_locally_solvable(p1)
     356        False
     357        sage: C.is_locally_solvable(p2)
     358        True
     359        sage: O = K.maximal_order()
     360        sage: f = (2*O).factor(); f
     361        (Fractional ideal (-a^2 - a + 1)) * (Fractional ideal (a^2 - 2*a + 3))
     362        sage: C.is_locally_solvable(f[0][0])
     363        True
     364        sage: C.is_locally_solvable(f[1][0])
     365        False
     366    """
     367    B = self.base_ring()
     368    D, T = self.diagonal_matrix()
     369    abc = [D[j, j] for j in range(3)]
     370    if abc[2] == 0:
     371        return True
     372    a = -abc[0]/abc[2]
     373    b = -abc[1]/abc[2]
     374    if is_RationalField(B):
     375        from sage.rings.arith import hilbert_symbol
     376        if is_RingHomomorphism(p):
     377            if p.domain == B and is_RealField(p.codomain):
     378                p = -1
     379            else:
     380                raise TypeError, "p (=%s) needs to be a prime of base field " \
     381                                 "B ( =`QQ`) in is_locally_solvable" % p
     382        if hilbert_symbol(a, b, p) == -1:
     383            if self._local_obstruction == None:
     384                self._local_obstruction = p
     385            return False
     386        return True
     387    if is_NumberField(B):
     388        from hilbert import hilbert_symbol_nf
     389        if hilbert_symbol_nf(a, b, p) == -1:
     390            if self._local_obstruction == None:
     391                self._local_obstruction = p
     392            return False
     393        return True
     394    raise TypeError, "Base field %s needs to be a number field or the " \
     395                     "rational numbers in is_locally_solvable" % B
     396
     397
     398def local_obstructions(self, finite = True, infinite = True, read_cache = True):
     399    r"""
     400    Returns the sequence of finite primes and/or infinite places
     401    such that self is locally solvable at those primes and places.
     402   
     403    If the base field is `QQ`, then the infinite place is denoted `-1`.
     404
     405    The parameters finite and infinite (both True by default) are
     406    used to specify whether to look at finite and/or infinite places.
     407    Note that finite = True involves factorization of the determinant
     408    of self, hence may be slow.
     409   
     410    Local obstructions are cached. The parameter read_cache can be used
     411    to specify whether to look at the cache before computing anything.
     412
     413    EXAMPLES:
     414
     415    Examples over QQ ::
     416
     417        sage: Conic(QQ, [1, 1, 1]).local_obstructions()
     418        [2, -1]
     419        sage: Conic(QQ, [1, 2, -3]).local_obstructions()
     420        []
     421        sage: Conic(QQ, [1, 2, 3, 4, 5, 6]).local_obstructions()
     422        [41, -1]
     423
     424
     425    Examples over number fields ::
     426
     427        sage: K.<i> = QuadraticField(-1)
     428        sage: Conic(K, [1, 1, 1]).local_obstructions()
     429        []
     430        sage: Conic(K, [1, 2, -3]).local_obstructions()
     431        []
     432        sage: Conic(K, [1, 2, 3, 4, 5, 6]).local_obstructions()
     433        [Fractional ideal (5*i - 4), Fractional ideal (-5*i - 4)]
     434
     435        sage: _.<x> = QQ[]
     436        sage: L.<b> = NumberField(x^4-2)
     437        sage: Conic(L, [1, 1, 1]).local_obstructions() # long time (2 seconds)
     438        [Ring morphism:
     439          From: Number Field in b with defining polynomial x^4 - 2
     440          To:   Real Field with 106 bits of precision
     441          Defn: b |--> -1.189207115002721066717492233629, Ring morphism:
     442          From: Number Field in b with defining polynomial x^4 - 2
     443          To:   Real Field with 106 bits of precision
     444          Defn: b |--> 1.189207115002721066717492233629]
     445        sage: Conic(L, [b, 1, 1]).local_obstructions() # long time (2 seconds)
     446        [Fractional ideal (b),
     447         Ring morphism:
     448         From: Number Field in b with defining polynomial x^4 - 2
     449         To:   Real Field with 106 bits of precision
     450         Defn: b |--> 1.189207115002721066717492233629]
     451    """
     452    obs0 = []
     453    obs1 = []
     454    B = self.base_ring()
     455    if infinite:
     456        if read_cache and self._infinite_obstructions != None:
     457            obs0 = self._infinite_obstructions
     458        else:
     459            if is_RationalField(B):
     460                if not self.is_locally_solvable(-1):
     461                    obs0 = [-1]
     462            elif is_NumberField(B):
     463                for b in B.places():
     464                    if not self.is_locally_solvable(b):
     465                        obs0.append(b)
     466            else:
     467                raise TypeError, "Base field (=%s) needs to be QQ or a" \
     468                                 "number field in local_obstructions" % B
     469            self._infinite_obstructions = obs0
     470    if finite:
     471        if read_cache and self._finite_obstructions != None:
     472            obs1 = self._finite_obstructions
     473        else:
     474            candidates = []
     475            if is_RationalField(B):
     476                if self.determinant() != 0:
     477                    for a in self.symmetric_matrix().list():
     478                        if a != 0:
     479                            for f in a.factor():
     480                                if f[1] < 0 and not f[0] in candidates:
     481                                    candidates.append(f[0])
     482                    for f in (2*self.determinant()).factor():
     483                        if f[1] > 0 and not f[0] in candidates:
     484                            candidates.append(f[0])
     485            elif is_NumberField(B):
     486                if self.determinant() != 0:
     487                    O = B.maximal_order()
     488                    for a in self.symmetric_matrix().list():
     489                        if a != 0:
     490                            for f in O.fractional_ideal(a).factor():
     491                                if f[1] < 0 and not f[0] in candidates:
     492                                    candidates.append(f[0])
     493                    for f in O.fractional_ideal(2*self.determinant()).factor():
     494                        if f[1] > 0 and not f[0] in candidates:
     495                            candidates.append(f[0])       
     496            else:
     497                raise TypeError, "Base field (=%s) needs to be QQ or a " \
     498                                 "number field in local_obstructions" % B
     499            for b in candidates:
     500                if not self.is_locally_solvable(b):
     501                   obs1.append(b)
     502            self._infinite_obstructions = obs1
     503    obs = obs1 + obs0
     504    if finite and infinite and len(obs)%2==1:
     505        raise RuntimeError, "Bug in local_obstructions or in " \
     506                            "is_locally_solvable: the number of bad places" \
     507                            "for conic self (=%s) is odd" % self
     508    return obs
     509
     510def parametrization(self, point=None, morphism=True):
     511    r"""
     512    Return a parametrization of self and the inverse
     513    of the parametrization.
     514
     515    If point is specified, then that point is used
     516    for the parametrization. Otherwise, use rational_point
     517    to find a point.
     518       
     519    If morphism is False, then returns the tuples of three
     520    polynomials that give the parametrization.
     521   
     522    Raises ValueError if no rational point exists or self is non-smooth.
     523
     524    ALGORITHM:
     525   
     526        Uses Denis Simon's Qfparam if the base field is `QQ`.
     527        Gives a straightforward non-optimized parametrization
     528        otherwise.
     529   
     530    EXAMPLES:
     531        sage: R.<x,y,z> = QQ[]
     532        sage: C = Curve(7*x^2 + 2*y*z + z^2)
     533        sage: (p, i) = C.parametrization(morphism = False); (p, i)
     534        ([-2*x*y, 7*x^2 + y^2, -2*y^2], [-1/2*x, -1/2*z])
     535        sage: C.defining_polynomial()(p)
     536        0
     537        sage: i[0](p) / i[1](p)
     538        x/y
     539
     540        sage: C = Curve(x^2 + 2*y^2 + z^2)
     541        sage: C.parametrization()
     542        Traceback (most recent call last):
     543        ...
     544        ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + 2*y^2 + z^2 has no rational points over Rational Field!
     545
     546        sage: C = Curve(x^2 + y^2 + 7*z^2)
     547        sage: C.parametrization()
     548        Traceback (most recent call last):
     549        ...
     550        ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + y^2 + 7*z^2 has no rational points over Rational Field!
     551
     552        sage: c = Conic([1,1,-1])
     553        sage: c.param()
     554        (Scheme morphism:
     555          From: Projective Space of dimension 1 over Rational Field
     556          To:   Projective Conic Curve over Rational Field defined by x^2 + y^2 - z^2
     557          Defn: Defined on coordinates by sending (x : y) to
     558                (2*x*y : x^2 - y^2 : x^2 + y^2),
     559         Scheme morphism:
     560           From: Projective Conic Curve over Rational Field defined by x^2 + y^2 - z^2
     561           To:   Projective Space of dimension 1 over Rational Field
     562           Defn: Defined on coordinates by sending (x : y : z) to
     563                 (1/2*x : -1/2*y + 1/2*z))
     564        sage: c = Conic(GF(2), [1,1,1,1,1,0])
     565        sage: c.param()
     566        (Scheme morphism:
     567          From: Projective Space of dimension 1 over Finite Field of size 2
     568          To:   Projective Conic Curve over Finite Field of size 2 defined by x^2 + x*y
     569        + y^2 + x*z + y*z
     570          Defn: Defined on coordinates by sending (x : y) to
     571                (x*y + y^2 : x^2 + x*y : x^2 + x*y + y^2),
     572         Scheme morphism:
     573          From: Projective Conic Curve over Finite Field of size 2 defined by x^2 + x*y
     574        + y^2 + x*z + y*z
     575          To:   Projective Space of dimension 1 over Finite Field of size 2
     576          Defn: Defined on coordinates by sending (x : y : z) to
     577                (y : x))
     578    """
     579    if (not self._parametrization is None) and not point:
     580        par = self._parametrization
     581    else:
     582        if not self.is_smooth():
     583            raise ValueError, "The conic self (=%s) is not smooth, hence does not have a parametrization.", self
     584        if point == None:
     585            point = self.rational_point()
     586        point = Sequence(point)
     587        B = self.base_ring()
     588        Q = PolynomialRing(B, 'x,y')
     589        [x, y] = Q.gens()
     590        gens = self.ambient_space().gens()
     591        if is_RationalField(B):
     592            from sage.quadratic_forms.qfsolve import Qfparam
     593            from sage.rings.arith import lcm
     594            M = self.symmetric_matrix()
     595            M *= lcm([ t.denominator() for t in M.list() ])
     596            par1 = Qfparam(M, point)
     597            B = Matrix([[par1[i][j] for j in range(3)] for i in range(3)])
     598            # self is in the image of B and does not lie on a line,
     599            # hence B is invertible
     600            A = B.inverse()
     601            par2 = [sum([A[i,j]*gens[j] for j in range(3)]) for i in [1,0]]
     602            par = ([Q(pol(x/y)*y**2) for pol in par1], par2)
     603        else:
     604            P = PolynomialRing(B, 4, ['X', 'Y', 'T0', 'T1'])
     605            [X, Y, T0, T1] = P.gens()
     606            c3 = [j for j in range(2,-1,-1) if point[j] != 0][0]
     607            c1 = [j for j in range(3) if j != c3][0]
     608            c2 = [j for j in range(3) if j != c3 and j != c1][0]
     609            L = [0,0,0]
     610            L[c1] = Y*T1*point[c1] + Y*T0
     611            L[c2] = Y*T1*point[c2] + X*T0
     612            L[c3] = Y*T1*point[c3]
     613            bezout = P(self.defining_polynomial()(L) / T0)
     614            t = [bezout([x,y,0,-1]),bezout([x,y,1,0])]
     615            par = (tuple([Q(p([x,y,t[0],t[1]])/y) for  p in L]),
     616                   tuple([gens[m]*point[c3]-gens[c3]*point[m]
     617                       for m in [c2,c1]]))
     618        if self._parametrization is None:
     619            self._parametrization = par
     620    if not morphism:
     621        return par
     622    P1 = ProjectiveSpace(self.base_ring(), 1, 'x,y')
     623    return P1.hom(par[0],self), self.hom(par[1],P1)
     624   
     625               
     626def rational_point(self, algorithm = 'default', read_cache = True):
     627    r"""Return a rational point (x0, y0, z0) on self.
     628
     629    Raises ValueError if no rational point exists.
     630
     631    See has_rational_point(self) for the algorithm used and for the use of
     632    the parameters algorithm and read_cache.
     633     
     634    EXAMPLES:
     635
     636    Examples over QQ ::
     637   
     638        sage: R.<x,y,z> = QQ[]
     639        sage: C = Curve(7*x^2 + 2*y*z + z^2)
     640        sage: C.rational_point()
     641        (0 : 1 : 0)
     642
     643        sage: C = Curve(x^2 + 2*y^2 + z^2)
     644        sage: C.rational_point()
     645        Traceback (most recent call last):
     646        ...
     647        ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + 2*y^2 + z^2 has no rational points over Rational Field!
     648
     649        sage: C = Curve(x^2 + y^2 + 7*z^2)
     650        sage: C.rational_point(algorithm = 'all')
     651        Traceback (most recent call last):
     652        ...
     653        ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + y^2 + 7*z^2 has no rational points over Rational Field!
     654
     655    Examples over number fields ::
     656       
     657        sage: _.<x> = QQ[]
     658        sage: L.<b> = NumberField(x^3-5)
     659        sage: C = Conic(L, [3, 2, -5])
     660        sage: C.rational_point(algorithm = 'all')  # long time (1 second)
     661        (-1 : -1 : 1)
     662   
     663    Examples over finite fields ::
     664   
     665        sage: F.<a> = FiniteField(7^20)
     666        sage: C = Conic([1, a, -5]); C
     667        Projective Conic Curve over Finite Field in a of size 7^20 defined by x^2 + (a)*y^2 + 2*z^2
     668        sage: C.rational_point()  # output is random
     669        (4*a^19 + 5*a^18 + 4*a^17 + a^16 + 6*a^15 + 3*a^13 + 6*a^11 + a^9 + 3*a^8 + 2*a^7 + 4*a^6 + 3*a^5 + 3*a^4 + a^3 + a + 6 : 5*a^18 + a^17 + a^16 + 6*a^15 + 4*a^14 + a^13 + 5*a^12 + 5*a^10 + 2*a^9 + 6*a^8 + 6*a^7 + 6*a^6 + 2*a^4 + 3 : 1)
     670         
     671    Examples over real and complex fields ::
     672
     673        sage: Conic(CC, [1, 2, 3]).rational_point()
     674        (0 : 1.22474487139159*I : 1)
     675
     676        sage: Conic(RR, [1, 1, 1]).rational_point()
     677        Traceback (most recent call last):
     678        ...
     679        ValueError: Conic Projective Conic Curve over Real Field with 53 bits of precision defined by x^2 + y^2 + z^2 has no rational points over Real Field with 53 bits of precision!
     680    """
     681    bl,pt = self.has_rational_point(point = True, algorithm = algorithm,
     682                                    read_cache = read_cache)
     683    if bl:
     684        return pt
     685    raise ValueError, "Conic %s has no rational points over %s!" % \
     686                      (self, self.ambient_space().base_ring())
     687
     688   
     689
  • new file sage/schemes/plane_conics/constructor.py

    diff -r 8dec8b43ccca -r 3fcab21a698f sage/schemes/plane_conics/constructor.py
    - +  
     1r"""
     2Plane conic constructors
     3
     4AUTHORS:
     5
     6- Marco Streng (2009-08-07)
     7
     8"""
     9
     10#*****************************************************************************
     11#       Copyright (C) 2009 Marco Streng <marco.streng@gmail.com>
     12#
     13#  Distributed under the terms of the GNU General Public License (GPL)
     14#
     15#    This code is distributed in the hope that it will be useful,
     16#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     17#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     18#    General Public License for more details.
     19#
     20#  The full text of the GPL is available at:
     21#
     22#                  http://www.gnu.org/licenses/
     23#*****************************************************************************
     24
     25from sage.matrix.constructor import Matrix
     26from sage.modules.free_module_element import vector
     27from sage.quadratic_forms.all import is_QuadraticForm
     28from sage.rings.all import (PolynomialRing, is_MPolynomial,
     29                            is_MPolynomialRing, is_Ring,
     30                            is_IntegralDomain)
     31from sage.schemes.generic.all import ProjectiveSpace
     32from sage.schemes.generic.morphism import (SchemeMorphism_affine_coordinates,
     33                            SchemeMorphism_projective_coordinates_field)
     34from sage.structure.all import Sequence
     35from sage.structure.element import is_Matrix
     36
     37from projective_conic import ProjectiveConic_generic
     38
     39def Conic(base_field, F=None, names=None, unique=True):
     40    r"""
     41    Return the plane projective conic curve defined by `F`
     42    over base_field, where `F` can be either a two- or
     43    three-variable polynomial, a list or matrix of coefficients,
     44    a ternary quadratic form, or a list or tuple of 5 points
     45    in the plane.
     46
     47    The input form Conic(F, names=None) is also accepted,
     48    in which case the fraction field of the base ring of `F`
     49    is used as base field.
     50
     51    The argument names is a list, tuple, or comma separated string
     52    of three variable names.
     53
     54    If `F` is a polynomial, then it needs to have degree 2
     55    and if it has three variables, then it needs to be
     56    homogeneous. If `F` is a polynomial or quadratic form,
     57    then the output is the curve in the projective plane
     58    defined by `F=0`.
     59
     60    If `F` is a matrix, then the output is the zero locus
     61    of `(x,y,z) F (x,y,z)^t`.
     62   
     63    If `F` is a list of coefficients, then it has
     64    length 3 or 6 and gives the coefficients of
     65    the monomials `x^2, y^2, z^2` or all 6 monomials
     66    in lexicographic order.
     67
     68    If `F` is a list of points in the plane, then the output
     69    is a conic through those points. If the argument unique
     70    is true and there is no unique conic through those points,
     71    then a ValueError is raised.
     72   
     73    EXAMPLE: Conic curves given by equations
     74   
     75    ::
     76
     77        sage: X,Y,Z = QQ['X,Y,Z'].gens()
     78        sage: Conic(X^2 - X*Y + Y^2 - Z^2)
     79        Projective Conic Curve over Rational Field defined by X^2 - X*Y + Y^2 - Z^2
     80        sage: x,y = GF(7)['x,y'].gens()
     81        sage: Conic(x^2 - x + 2*y^2 - 3, 'U,V,W')
     82        Projective Conic Curve over Finite Field of size 7 defined by U^2 + 2*V^2 - U*W - 3*W^2
     83
     84    EXAMPLE: Conic curves given by matrices
     85
     86    ::
     87
     88        sage: Conic(matrix(QQ, [[1, 2, 0], [4, 0, 0], [7, 0, 9]]), 'x,y,z')
     89        Projective Conic Curve over Rational Field defined by x^2 + 6*x*y + 7*x*z + 9*z^2
     90
     91        sage: x,y,z = GF(11)['x,y,z'].gens()
     92        sage: C = Conic(x^2+y^2-2*z^2); C
     93        Projective Conic Curve over Finite Field of size 11 defined by x^2 + y^2 - 2*z^2
     94        sage: Conic(C.symmetric_matrix(), 'x,y,z')
     95        Projective Conic Curve over Finite Field of size 11 defined by x^2 + y^2 - 2*z^2
     96
     97    EXAMPLE: Conics given by coefficients
     98
     99    ::
     100   
     101        sage: Conic(QQ, [1,2,3])
     102        Projective Conic Curve over Rational Field defined by x^2 + 2*y^2 + 3*z^2
     103        sage: Conic(GF(7), [1,2,3,4,5,6], 'X')
     104        Projective Conic Curve over Finite Field of size 7 defined by X0^2 + 2*X0*X1 - 3*X1^2 + 3*X0*X2 - 2*X1*X2 - X2^2
     105   
     106    EXAMPLE: The conic through a given set of points
     107
     108    ::
     109
     110        sage: C = Conic(QQ, [[10,2],[3,4],[-7,6],[7,8],[9,10]]); C
     111        Projective Conic Curve over Rational Field defined by x^2 + 13/4*x*y - 17/4*y^2 - 35/2*x*z + 91/4*y*z - 37/2*z^2
     112        sage: C.rational_point()
     113        (10 : 2 : 1)
     114        sage: C.point([3,4])
     115        (3 : 4 : 1)
     116
     117        sage: a=AffineSpace(GF(13),2)
     118        sage: Conic([a([x,x^2]) for x in range(5)])
     119        Projective Conic Curve over Finite Field of size 13 defined by x^2 - y*z
     120    """
     121    if not (is_IntegralDomain(base_field) or base_field == None):
     122        if names is None:
     123            names = F
     124        F = base_field
     125        base_field = None
     126    if isinstance(F, (list,tuple)):
     127        if len(F) == 1:
     128            return Conic(base_field, F[0], names)
     129        if names == None:
     130            names = 'x,y,z'
     131        if len(F) == 5:
     132            L=[]
     133            for f in F:
     134                if isinstance(f, SchemeMorphism_affine_coordinates):
     135                    C = Sequence(f, universe = base_field)
     136                    if len(C) != 2:
     137                        raise TypeError, "points in F (=%s) must be planar"%F
     138                    C.append(1)
     139                elif isinstance(f, SchemeMorphism_projective_coordinates_field):
     140                    C = Sequence(f, universe = base_field)
     141                elif isinstance(f, (list, tuple)):
     142                    C = Sequence(f, universe = base_field)
     143                    if len(C) == 2:
     144                        C.append(1)
     145                else:
     146                    raise TypeError, "F (=%s) must be a sequence of planar " \
     147                                      "points" % F
     148                if len(C) != 3:
     149                    raise TypeError, "points in F (=%s) must be planar" % F
     150                P = C.universe()
     151                if not is_IntegralDomain(P):
     152                    raise TypeError, "coordinates of points in F (=%s) must " \
     153                                     "be in an integral domain" % F
     154                L.append(Sequence([C[0]**2, C[0]*C[1], C[0]*C[2], C[1]**2,
     155                                   C[1]*C[2], C[2]**2], P.fraction_field()))
     156            M=Matrix(L)
     157            if unique and M.rank() != 5:
     158                raise ValueError, "points in F (=%s) do not define a unique " \
     159                                   "conic" % F
     160            con = Conic(base_field, Sequence(M.right_kernel().gen()), names)
     161            _ = con.point(F[0])
     162            return con
     163        F = Sequence(F, universe = base_field)
     164        base_field = F.universe().fraction_field()
     165        temp_ring = PolynomialRing(base_field, 3, names)
     166        (x,y,z) = temp_ring.gens()
     167        if len(F) == 3:
     168            return Conic(F[0]*x**2 + F[1]*y**2 + F[2]*z**2)
     169        if len(F) == 6:
     170            return Conic(F[0]*x**2 + F[1]*x*y + F[2]*x*z + F[3]*y**2 + \
     171                         F[4]*y*z + F[5]*z**2)
     172        raise TypeError, "F (=%s) must be a sequence of 3 or 6" \
     173                         "coefficients" % F
     174    if is_QuadraticForm(F):
     175        F = F.matrix()
     176    if is_Matrix(F) and F.is_square() and F.ncols() == 3:
     177        if names == None:
     178            names = 'x,y,z'
     179        temp_ring = PolynomialRing(F.base_ring(), 3, names)
     180        F = vector(temp_ring.gens()) * F * vector(temp_ring.gens())
     181
     182    if not is_MPolynomial(F):
     183        raise TypeError, "F (=%s) must be a three-variable polynomial or " \
     184                         "a sequence of points or coefficients" % F
     185
     186    if F.total_degree() != 2:
     187        return TypeError, "F (=%s) must have degree 2" % F
     188
     189    if base_field == None:
     190        base_field = F.base_ring()
     191    if not is_IntegralDomain(base_field):
     192        raise ValueError, "Base field (=%s) must be a field" % base_field
     193    base_field = base_field.fraction_field()
     194    if names == None:
     195        names = F.parent().variable_names()
     196    pol_ring = PolynomialRing(base_field, 3, names)
     197
     198    if F.parent().ngens() == 2:
     199        (x,y,z) = pol_ring.gens()
     200        F = pol_ring(F(x/z,y/z)*z**2)   
     201
     202    if F == 0:
     203        raise ValueError, "F must be nonzero over base field %s" % base_field
     204
     205    if F.total_degree() != 2:
     206        return TypeError, "F (=%s) must have degree 2 over base field %s" % \
     207                          (F, base_field)
     208
     209    if F.parent().ngens() == 3:
     210        P2 = ProjectiveSpace(2, base_field, names)
     211        return ProjectiveConic_generic(P2, F)
     212
     213    raise TypeError, "Number of variables of F (=%s) must be 2 or 3" % F
  • new file sage/schemes/plane_conics/hilbert.py

    diff -r 8dec8b43ccca -r 3fcab21a698f sage/schemes/plane_conics/hilbert.py
    - +  
     1r"""
     2Compute Hilbert symbols over number fields.
     3
     4AUTHORS:
     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
     24from sage.rings.all import (ZZ, is_NumberFieldElement, is_RingHomomorphism,
     25                            is_RealField, is_ComplexField, is_NumberFieldIdeal)
     26
     27
     28def 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
  • new file sage/schemes/plane_conics/projective_conic.py

    diff -r 8dec8b43ccca -r 3fcab21a698f sage/schemes/plane_conics/projective_conic.py
    - +  
     1r"""
     2Projective plane conics over a field.
     3
     4AUTHORS:
     5
     6- Marco Streng (2009-08-07)
     7
     8- Nick Alexander (2008-01-08)
     9
     10"""
     11#*****************************************************************************
     12#       Copyright (C) 2008 Nick Alexander <ncalexander@gmail.com>
     13#       Copyright (C) 2009 Marco Streng <marco.streng@gmail.com>
     14#
     15#  Distributed under the terms of the GNU General Public License (GPL)
     16#
     17#    This code is distributed in the hope that it will be useful,
     18#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     19#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     20#    General Public License for more details.
     21#
     22#  The full text of the GPL is available at:
     23#
     24#                  http://www.gnu.org/licenses/
     25#*****************************************************************************
     26
     27from sage.libs.pari.gen import pari
     28from sage.misc.all import add, sage_eval
     29from sage.rings.all import (PolynomialRing, ZZ, QQ, MPolynomialRing,
     30                            degree_lowest_rational_function,
     31                            is_PrimeField, is_FiniteField,
     32                            is_ComplexField, is_RealField,
     33                            is_pAdicField, is_Field,
     34                            is_RationalField)
     35from sage.modules.free_module_element import vector
     36from sage.structure.sequence import Sequence
     37from sage.structure.element import is_Vector
     38from sage.schemes.generic.projective_space import (ProjectiveSpace,
     39                                                   is_ProjectiveSpace)
     40from sage.matrix.constructor import Matrix
     41from sage.matrix.matrix import is_Matrix
     42
     43from sage.schemes.plane_curves.curve import Curve_generic_projective
     44from sage.schemes.plane_curves.projective_curve import ProjectiveCurve_generic
     45from sage.quadratic_forms.qfsolve import Qfsolve, Qfparam
     46
     47
     48class ProjectiveConic_generic(ProjectiveCurve_generic):
     49    def __init__(self, A, f):
     50        ProjectiveCurve_generic.__init__(self, A, f)
     51        self._coefficients = [f[(2,0,0)], f[(1,1,0)], f[(1,0,1)],
     52                                f[(0,2,0)], f[(0,1,1)], f[(0,0,2)]]
     53        self._rational_point = None
     54        self._parametrization = None
     55        self._diagonal_matrix = None
     56       
     57        # a single prime such that self has no point over the completion
     58        self._local_obstruction = None
     59        # all finite primes such that self has no point over the completion
     60        self._finite_obstructions = None
     61        # all infinite primes such that self has no point over the completion
     62        self._infinite_obstructions = None
     63
     64    from conic_solve import (has_rational_point, rational_point,
     65                             parametrization, local_obstructions,
     66                             is_locally_solvable)
     67
     68    param = parametrization
     69
     70
     71    def _repr_type(self):
     72        return "Projective Conic"
     73
     74    def _vector_(self):
     75        return vector(self.coefficients())
     76
     77    def base_extend(self, S):
     78        r"""
     79        Returns the conic over `S`, given by the same equation as self.
     80
     81        EXAMPLES:
     82
     83        ::
     84
     85            sage: c = Conic([1, 1, 1]); c
     86            Projective Conic Curve over Rational Field defined by x^2 + y^2 + z^2
     87            sage: c.has_rational_point()
     88            False
     89            sage: d = c.base_extend(QuadraticField(-1, 'i')); d
     90            Projective Conic Curve over Number Field in i with defining polynomial x^2 + 1 defined by x^2 + y^2 + z^2
     91            sage: d.rational_point()
     92            (-i : 1 : 0)
     93        """
     94        if is_Field(S):
     95            B = self.base_ring()
     96            if B == S:
     97                return self
     98            if not S.has_coerce_map_from(B):
     99                raise ValueError, "No natural map from the base ring of self " \
     100                                  "(= %s) to S (= %s)" % (self, S)
     101            from constructor import Conic
     102            con = Conic([S(c) for c in self.coefficients()], \
     103                        self.variable_names())
     104            if self._rational_point != None:
     105                pt = [S(c) for c in Sequence(self._rational_point)]
     106                if not pt == [0,0,0]:
     107                    _ = con.point(pt)
     108            return con
     109        return ProjectiveCurve_generic.base_extend(self, S)
     110   
     111    def cache_point(self, p):
     112        r"""
     113        Replace the point in the cache of self by `p` for use
     114        by self.rational_point and self.parametrization.
     115       
     116        EXAMPLES:
     117
     118        ::
     119
     120            sage: c = Conic([1, -1, 1])
     121            sage: c.point([15, 17, 8])
     122            (15/8 : 17/8 : 1)
     123            sage: c.rational_point()
     124            (15/8 : 17/8 : 1)
     125            sage: c.cache_point(c.rational_point(read_cache = False))
     126            sage: c.rational_point()
     127            (1 : 1 : 0)
     128        """
     129        if isinstance(p, (tuple, list)):
     130            p = self.point(p)
     131        self._rational_point = p
     132
     133    def coefficients(self):
     134        r"""
     135        Gives a the 6 coefficients of the conic
     136        in lexicographic order.
     137       
     138        EXAMPLES:
     139       
     140        ::
     141       
     142            sage: Conic(QQ, [1,2,3,4,5,6]).coefficients()
     143            [1, 2, 3, 4, 5, 6]
     144
     145            sage: P.<x,y,z> = GF(13)[]
     146            sage: a = Conic(x^2+5*x*y+y^2+z^2).coefficients(); a
     147            [1, 5, 0, 1, 0, 1]
     148            sage: Conic(a)
     149            Projective Conic Curve over Finite Field of size 13 defined by x^2 + 5*x*y + y^2 + z^2
     150        """
     151        return self._coefficients
     152       
     153    def count_points(self, n):
     154        r"""
     155        If the base field `B` of `self` is finite of order `q`,
     156        then returns the number of points over `GF{q}, ..., GF{q^n}`.
     157       
     158        EXAMPLES:
     159
     160            sage: P.<x,y,z> = GF(3)[]
     161            sage: c = Conic(x^2+y^2+z^2)
     162            sage: c.count_points(4)
     163            [4, 10, 28, 82]
     164        """       
     165        F = self.base_ring()
     166        if not F.is_finite():
     167            raise TypeError, "Point counting only defined for schemes over " \
     168                             "finite fields, not over the base ring of self " \
     169                             "(= %s)" % self
     170        q = F.cardinality()
     171        return [q**i+1 for i in range(1, n+1)]
     172       
     173    def derivative_matrix(self):
     174        r"""
     175        Gives the derivative of the defining polynomial,
     176        which is a linear map, as a `3 \times 3` matrix.
     177       
     178        Examples:
     179       
     180        In characteristic different from 2, the
     181        derivative matrix is twice the symmetric matrix:
     182
     183        ::
     184       
     185            sage: c = Conic(QQ, [1,1,1,1,1,0])
     186            sage: c.symmetric_matrix()
     187            [  1 1/2 1/2]
     188            [1/2   1 1/2]
     189            [1/2 1/2   0]
     190            sage: c.derivative_matrix()
     191            [2 1 1]
     192            [1 2 1]
     193            [1 1 0]
     194
     195        An example in characteristic 2:
     196
     197        ::
     198
     199            sage: P.<t> = GF(2)[]
     200            sage: c = Conic([t, 1, t^2, 1, 1, 0]); c
     201            Projective Conic Curve over Fraction Field of Univariate Polynomial Ring in t over Finite Field of size 2 (using NTL) defined by t*x^2 + x*y + y^2 + t^2*x*z + y*z
     202            sage: c.is_smooth()
     203            True
     204            sage: c.derivative_matrix()
     205            [  0   1 t^2]
     206            [  1   0   1]
     207            [t^2   1   0]
     208        """
     209        from sage.matrix.constructor import matrix
     210        [a,b,c,d,e,f] = self.coefficients()
     211        return matrix([[ 2*a ,   b ,   c ],
     212                       [   b , 2*d ,   e ],
     213                       [   c ,   e , 2*f ]])
     214       
     215    def determinant(self):
     216        r"""
     217        Gives the determinant of the symmetric matrix that defines the conic.
     218        This is defined only if the base field has characteristic different
     219        from 2.
     220
     221        EXAMPLES:
     222
     223        ::
     224
     225            sage: C = Conic([1,2,3,4,5,6])
     226            sage: C.determinant()
     227            41/4
     228            sage: C.symmetric_matrix().determinant()
     229            41/4
     230
     231        Determinants are only defined in characteristic different from 2::
     232
     233            sage: C = Conic(GF(2), [1, 1, 1, 1, 1, 0])
     234            sage: C.is_smooth()
     235            True
     236            sage: C.determinant()
     237            Traceback (most recent call last):
     238            ...
     239            ValueError: The conic self (= Projective Conic Curve over Finite Field of size 2 defined by x^2 + x*y + y^2 + x*z + y*z) has no symmetric matrix because the base field has characteristic 2
     240        """
     241        return self.symmetric_matrix().determinant()
     242
     243    det = determinant
     244
     245    def diagonal_matrix(self):
     246        r"""
     247        Returns a diagonal matrix `D` and a matrix `T` such that `T^t A T = D`
     248        holds, where `(x, y, z) A (x, y, z)^t` is the defining polynomial
     249        of the conic `self`.
     250
     251        EXAMPLE:
     252
     253        ::
     254
     255            sage: c = Conic(QQ, [1,2,3,4,5,6])
     256            sage: d, t = c.diagonal_matrix(); d, t
     257            ([    1     0     0]
     258            [    0     3     0]
     259            [    0     0 41/12],
     260             [   1   -1 -7/6]
     261             [   0    1 -1/3]
     262             [   0    0    1])
     263            sage: t.transpose()*c.symmetric_matrix()*t
     264            [    1     0     0]
     265            [    0     3     0]
     266            [    0     0 41/12]
     267         
     268        Diagonal matrices are only defined in characteristic different
     269        from 2:
     270       
     271        ::
     272
     273            sage: c = Conic(GF(4, 'a'), [0, 1, 1, 1, 1, 1])
     274            sage: c.is_smooth()
     275            True
     276            sage: c.diagonal_matrix()
     277            Traceback (most recent call last):
     278            ...
     279            ValueError: The conic self (= Projective Conic Curve over Finite Field in a of size 2^2 defined by x*y + y^2 + x*z + y*z + z^2) has no symmetric matrix because the base field has characteristic 2
     280        """
     281        A = self.symmetric_matrix()
     282        B = self.base_ring()
     283        basis = [vector(B,{2:0,i:1}) for i in range(3)]
     284        for i in range(3):
     285            zerovalue = (basis[i]*A*basis[i].transpose()== 0)
     286            if zerovalue:
     287                for j in range(i+1,3):
     288                    if basis[j]*A*basis[j].transpose() != 0:
     289                        b = basis[i]
     290                        basis[i] = basis[j]
     291                        basis[j] = b
     292                        zerovalue = False
     293            if zerovalue:
     294                for j in range(i+1,3):
     295                    if basis[i]*A*basis[j].transpose() != 0:
     296                        basis[i] = basis[i]+basis[j]
     297                        zerovalue = False
     298            if not zerovalue:
     299                l = (basis[i]*A*basis[i].transpose())
     300                for j in range(i+1,3):
     301                    basis[j] = basis[j] - \
     302                               (basis[i]*A*basis[j].transpose())/l * basis[i]
     303        T = Matrix(basis).transpose()
     304        return T.transpose()*A*T, T
     305       
     306    def diagonalization(self,names = None):
     307        r"""
     308        Returns a diagonal conic C, an isomorphism of schemes M: C -> self
     309        and the inverse N of M.
     310       
     311        EXAMPLES:
     312       
     313        ::
     314            sage: Conic(GF(5), [1,0,1,1,0,1]).diagonalization()           
     315            (Projective Conic Curve over Finite Field of size 5 defined by x^2 + y^2 + 2*z^2,
     316             Scheme morphism:
     317              From: Projective Conic Curve over Finite Field of size 5 defined by x^2 + y^2 + 2*z^2
     318              To:   Projective Conic Curve over Finite Field of size 5 defined by x^2 + y^2 + x*z + z^2
     319              Defn: Defined on coordinates by sending (x : y : z) to
     320                    (x + 2*z : y : z),
     321             Scheme morphism:
     322              From: Projective Conic Curve over Finite Field of size 5 defined by x^2 + y^2 + x*z + z^2
     323              To:   Projective Conic Curve over Finite Field of size 5 defined by x^2 + y^2 + 2*z^2
     324              Defn: Defined on coordinates by sending (x : y : z) to
     325                    (x - 2*z : y : z))
     326
     327        The diagonalization is only defined in characteristic different
     328        from 2:
     329       
     330        ::
     331
     332            sage: Conic(GF(2), [1,1,1,1,1,0]).diagonalization()
     333            Traceback (most recent call last):
     334            ...
     335            ValueError: The conic self (= Projective Conic Curve over Finite Field of size 2 defined by x^2 + x*y + y^2 + x*z + y*z) has no symmetric matrix because the base field has characteristic 2
     336        """
     337        if names == None:
     338            names = self.defining_polynomial().parent().variable_names()
     339        from constructor import Conic
     340        D, T = self.diagonal_matrix()
     341        con = Conic(D, names = names)
     342        return con, con.hom(T, self), self.hom(T.inverse(), con)
     343       
     344    def gens(self):
     345        r"""
     346        Returns the generators of the coordinate ring of self.
     347       
     348        EXAMPLES:
     349            sage: P.<x,y,z> = QQ[]
     350            sage: c = Conic(x^2+y^2+z^2)
     351            sage: c.gens()
     352            (xbar, ybar, zbar)
     353            sage: c.defining_polynomial()(c.gens())
     354            0
     355        """
     356        return self.coordinate_ring().gens()
     357       
     358    def has_singular_point(self, point = False):
     359        r"""
     360        Return True if and only if self has a rational
     361        singular point.
     362       
     363        If point is True, then also return a rational singular
     364        point (or None if no such point exists).
     365       
     366        EXAMPLES:
     367       
     368        ::
     369           
     370            sage: c = Conic(QQ, [1,0,1]); c
     371            Projective Conic Curve over Rational Field defined by x^2 + z^2
     372            sage: c.has_singular_point(point = True)
     373            (True, (0 : 1 : 0))
     374           
     375            sage: P.<x,y,z> = GF(7)[]
     376            sage: e = Conic((x+y+z)*(x-y+2*z)); e
     377            Projective Conic Curve over Finite Field of size 7 defined by x^2 - y^2 + 3*x*z + y*z + 2*z^2
     378            sage: e.has_singular_point(point = True)
     379            (True, (2 : 4 : 1))
     380
     381            sage: Conic([1, 1, -1]).has_singular_point()
     382            False
     383            sage: Conic([1, 1, -1]).has_singular_point(point = True)
     384            (False, None)
     385
     386        has_singular_point is not implemented over all fields
     387        of characteristic 2. It is implemented over finite fields.
     388
     389        ::
     390
     391            sage: F.<a> = FiniteField(8)
     392            sage: Conic([a, a+1, 1]).has_singular_point(point = True)
     393            (True, (a + 1 : 0 : 1))
     394
     395            sage: P.<t> = GF(2)[]
     396            sage: C = Conic(P, [t,t,1]); C
     397            Projective Conic Curve over Fraction Field of Univariate Polynomial Ring in t over Finite Field of size 2 (using NTL) defined by t*x^2 + t*y^2 + z^2
     398            sage: C.has_singular_point(point = False)
     399            Traceback (most recent call last):
     400            ...
     401            NotImplementedError: Sorry, find singular point on conics not implemented over all fields of characteristic 2.
     402        """
     403        if not point:
     404           ret = self.has_singular_point(point = True)
     405           return ret[0]
     406        B = self.base_ring()
     407        if B.characteristic() == 2:
     408            [a,b,c,d,e,f] = self.coefficients()
     409            if b == 0 and c == 0 and e == 0:
     410                for i in range(3):
     411                    if [a, d, f][i] == 0:
     412                        return True, self.point(vector(B, {2:0, i:1}))
     413                if hasattr(a/f, 'is_square') and hasattr(a/f, 'sqrt'):
     414                    if (a/f).is_square():
     415                        return True, self.point([1,0,(a/f).sqrt()])
     416                    if (d/f).is_square():
     417                        return True, self.point([0,1,(d/f).sqrt()])
     418                raise NotImplementedError, "Sorry, find singular point on conics not implemented over all fields of characteristic 2."
     419            pt = [e, c, b]
     420            if self.defining_polynomial()(pt) == 0:
     421                return True, self.point(pt)
     422            return False, None
     423        D = self.symmetric_matrix()
     424        if D.det() == 0:
     425            return True, self.point(Sequence(D.right_kernel().gen()))
     426        return False, None
     427
     428    def hom(self, x, Y=None):
     429        r"""
     430        Return the scheme morphism from self to Y defined by x.
     431        Here x can be a matrix or a sequence of polynomials.
     432        If Y is omitted, then a natural image is found if possible.
     433
     434        EXAMPLES:
     435       
     436        Morphisms given by matrices. In the first example, Y is omitted,
     437        in the second example, Y is specified.
     438
     439        ::
     440           
     441            sage: c = Conic([-1, 1, 1])
     442            sage: h = c.hom(Matrix([[1,1,0],[0,1,0],[0,0,1]])); h
     443            Scheme morphism:
     444              From: Projective Conic Curve over Rational Field defined by -x^2 + y^2 + z^2
     445              To:   Projective Conic Curve over Rational Field defined by -x^2 + 2*x*y + z^2
     446              Defn: Defined on coordinates by sending (x : y : z) to
     447                    (x + y : y : z)
     448            sage: h([-1, 1, 0])
     449            (0 : 1 : 0)
     450
     451            sage: c = Conic([-1, 1, 1])
     452            sage: d = Conic([4, 1, -1])
     453            sage: c.hom(Matrix([[0, 0, 1/2], [0, 1, 0], [1, 0, 0]]), d)
     454            Scheme morphism:
     455              From: Projective Conic Curve over Rational Field defined by -x^2 + y^2 + z^2
     456                To:   Projective Conic Curve over Rational Field defined by 4*x^2 + y^2 - z^2
     457                  Defn: Defined on coordinates by sending (x : y : z) to
     458                          (1/2*z : y : x)
     459
     460        A morphism given by polynomials:
     461
     462        ::
     463
     464            sage: c = Conic(GF(13), [-1, 1, 1])
     465            sage: P.<x, y, z> = GF(13)[]
     466            sage: h = c.hom([x + z, y - z], ProjectiveSpace(GF(13), 1)); h
     467            Scheme morphism:
     468              From: Projective Conic Curve over Finite Field of size 13 defined by -x^2 + y^2 + z^2
     469              To:   Projective Space of dimension 1 over Finite Field of size 13
     470              Defn: Defined on coordinates by sending (x : y : z) to
     471                    (x + z : y - z)
     472            sage: h([1, 1, 0])
     473            (1 : 1)
     474            sage: h([1, 0, 1])
     475            (11 : 1)
     476
     477        A ValueError is raised if the wrong codomain Y is specified:
     478
     479        ::
     480       
     481            sage: c = Conic([-1, 1, 1])
     482            sage: c.hom(Matrix([[0, 0, 1/2], [0, 1, 0], [1, 0, 0]]), c)
     483            Traceback (most recent call last):
     484            ...
     485            ValueError: The matrix x (= [  0   0 1/2]
     486            [  0   1   0]
     487            [  1   0   0]) does not define a map from self (= Projective Conic Curve over Rational Field defined by -x^2 + y^2 + z^2) to Y (= Projective Conic Curve over Rational Field defined by -x^2 + y^2 + z^2)
     488        """
     489        if is_Matrix(x):   
     490            from constructor import Conic
     491            y = x.inverse()
     492            A = y.transpose()*self.matrix()*y
     493            im = Conic(A)
     494            if Y != None:
     495                q = Y.defining_polynomial()/im.defining_polynomial()
     496                if not (q.numerator().is_constant()
     497                        and q.denominator().is_constant()):
     498                    raise ValueError, "The matrix x (= %s) does not define a " \
     499                                      "map from self (= %s) to Y (= %s)" % \
     500                                      (x, self, Y)
     501                im = Y
     502            return self.hom(Sequence(x*vector(self.ambient_space().gens())), im)
     503        return ProjectiveCurve_generic.hom(self, x, Y)           
     504           
     505    def irreducible_components(self):
     506        r"""
     507        Return the irreducible components of this conic, as
     508        subschemes of the same ambient space.
     509       
     510        EXAMPLES:
     511        ::
     512            sage: P.<x,y,z> = GF(7)[]
     513            sage: e = Conic((x+y+z)*(x-y+2*z)); e
     514            Projective Conic Curve over Finite Field of size 7 defined by x^2 - y^2 + 3*x*z
     515            + y*z + 2*z^2
     516           
     517            sage: e.irreducible_components()
     518            [
     519            Closed subscheme of Projective Space of dimension 2 over Finite Field of size 7
     520            defined by:
     521              x - y + 2*z,
     522            Closed subscheme of Projective Space of dimension 2 over Finite Field of size 7
     523            defined by:
     524              x + y + z
     525            ]
     526           
     527            sage: Conic(GF(7), [1,1,1]).irreducible_components()
     528            [Closed subscheme of Projective Space of dimension 2 over Finite Field of size 7
     529             defined by:
     530              x^2 + y^2 + z^2]
     531        """
     532        if self.is_smooth():
     533            return [self.ambient_space().subscheme(self.defining_polynomial())]
     534        return ProjectiveCurve_generic.irreducible_components(self)
     535       
     536    def is_diagonal(self):
     537        r"""
     538        Return True if and only if the conic has the form
     539        `a*x^2 + b*y^2 + c*z^2`.
     540
     541        EXAMPLES:
     542
     543        ::
     544
     545            sage: c=Conic([1,1,0,1,0,1]); c
     546            Projective Conic Curve over Rational Field defined by x^2 + x*y + y^2 + z^2
     547            sage: d,t = c.diagonal_matrix()
     548            sage: c.is_diagonal()
     549            False
     550            sage: c.diagonalization()[0].is_diagonal()
     551            True
     552        """
     553        return all([self.coefficients()[i] == 0 for i in [1,2,4]])
     554
     555    def is_isomorphic(self, c, isomorphism = False):
     556        r"""
     557        Return True if and only if self is isomorphic to
     558        the conic or projective line c.
     559       
     560        If the parameter isomorphic is True, then also return
     561        an isomorphism from self to c and its inverse.
     562        """
     563        B = self.base_field()
     564        Bc = c.base_field()
     565        if not B == Bc:
     566            raise TypeError, "c (= %s) must be a curve with the same base " \
     567                             "field as self in is_isomorphic" % c
     568        if is_ProjectiveSpace(c) and c.dimension == 1:
     569            if self.has_rational_point():
     570                if isomorphism:
     571                    par = self.parametrization()
     572                    return True, par[1], par[0]
     573                return True
     574            if isomorphism:
     575                return False, None
     576            return False
     577        if not isinstance(c, sage.schemes.plane_curves.projective_curve.
     578                                  ProjectiveConic_generic):
     579            raise TypeError, "c (= %s) must be a conic or projective line in " \
     580                             "is_isomorphic" % c
     581        obs = self.local_obstructions(finite = False)
     582        obsc = c.local_obstructions(finite = False)
     583        if not (all([o in obs for o in obsc]) and
     584                all([o in obsc for o in obs])):
     585            if isomorphism:
     586                return False, None
     587            return False
     588        obs = self.local_obstructions()
     589        obsc = c.local_obstructions()
     590        if not (all([o in obs for o in obsc]) and
     591                all([o in obsc for o in obs])):
     592            if isomorphism:
     593                return False, None
     594            return False
     595        if isomorphism:
     596            if obs != 0:
     597                raise NotImplementedError, "Find isomorphism for conics " \
     598                                           "without rational points not " \
     599                                           "yet implemented"
     600            par = self.parametrization()
     601            parc = c.parametrization()
     602            return (True, self.hom([s(par[1].defining_polynomials())
     603                         for s in parc[0]], c),
     604                         c.hom([s(parc[1].defining_polynomials())
     605                         for s in par[1]], self))
     606        return True
     607
     608    def is_smooth(self):
     609        r"""
     610        Returns True if and only if self is smooth.
     611
     612        EXAMPLES:
     613
     614        ::
     615
     616            sage: Conic([1,-1,0]).is_smooth()
     617            False
     618            sage: Conic(GF(2),[1,1,1,1,1,0]).is_smooth()
     619            True
     620        """
     621        if self.base_ring().characteristic() == 2:
     622            [a,b,c,d,e,f] = self.coefficients()
     623            if b == 0 and c == 0 and e == 0:
     624                return False
     625            return self.defining_polynomial()([e, c, b]) != 0
     626        return self.determinant() != 0
     627       
     628    def matrix(self):
     629        r"""
     630        Returns a matrix `M` such that `(x, y, z) M (x, y, z)^t`
     631        is the defining equation of self.
     632        The matrix `M` is upper triangular if the base field has
     633        characteristic 2 and symmetric otherwise.
     634
     635        EXAMPLES:
     636            sage: R.<x, y, z> = QQ[]
     637            sage: C = Conic(x^2 + x*y + y^2 + z^2)
     638            sage: C.matrix()
     639            [  1 1/2   0]
     640            [1/2   1   0]
     641            [  0   0   1]
     642
     643            sage: R.<x, y, z> = GF(2)[]
     644            sage: C = Conic(x^2 + x*y + y^2 + x*z + z^2)
     645            sage: C.matrix()
     646            [1 1 1]
     647            [0 1 0]
     648            [0 0 1]
     649        """
     650        if self.base_ring().characteristic() == 2:
     651            return self.upper_triangular_matrix()
     652        return self.symmetric_matrix()
     653       
     654    _matrix_ = matrix
     655
     656    def point(self, v, check=True):
     657        r"""
     658        Constructs a point on self corresponding to the input.
     659       
     660        If no rational point on self is known yet, then also caches the point
     661        for use by self.rational_point and self.parametrization
     662
     663        EXAMPLES:
     664
     665        ::
     666
     667            sage: c = Conic([1, -1, 1])
     668            sage: c.point([15, 17, 8])
     669            (15/8 : 17/8 : 1)
     670            sage: c.rational_point()
     671            (15/8 : 17/8 : 1)
     672            sage: d = Conic([1, -1, 1])
     673            sage: d.rational_point()
     674            (1 : 1 : 0)
     675        """
     676        if is_Vector(v):
     677            v = Sequence(v)
     678        p = ProjectiveCurve_generic.point(self, v, check=check)
     679        if self._rational_point is None:
     680            self._rational_point = p
     681        return p
     682
     683    def random_rational_point(self, bound=None):
     684        r"""
     685        Return a random rational point of this conic.
     686       
     687        The output is the image of a random point X on the
     688        projective line under a parametrization of self.
     689        If a bound is specified and the base field is `\QQ`,
     690        then that is a bound on the naive height of X.
     691        If the base field is a finite field, then the
     692        output is uniformly distributed over the rational
     693        points of self.
     694
     695        EXAMPLES:
     696           
     697        ::
     698            sage: c = Conic(GF(2), [1,1,1,1,1,0])
     699            sage: [c.random_rational_point() for _ in range(10)] # output is random
     700           
     701            [(0 : 0 : 1),
     702             (1 : 0 : 1),
     703             (1 : 0 : 1),
     704             (0 : 0 : 1),
     705             (1 : 0 : 1),
     706             (0 : 1 : 1),
     707             (0 : 1 : 1),
     708             (0 : 0 : 1),
     709             (0 : 0 : 1),
     710             (1 : 0 : 1)]
     711        """
     712        if not self.is_smooth():
     713            raise NotImplementedError, "Sorry, random points not implemented " \
     714                                       "for non-smooth conics"
     715        par = self.parametrization()
     716        x = 0
     717        y = 0
     718        B = self.base_ring()
     719        if is_RationalField(B) and bound != None:
     720            if bound < 1:
     721                raise ValueError, "The bound (= %s) needs to be at least 1" % \
     722                                  bound
     723            while x == 0 and y == 0:
     724                x = ZZ.random_element(-bound, bound+1)
     725                y = ZZ.random_element(-bound, bound+1)
     726        else:
     727            while x == 0 and y == 0:
     728                x = B.random_element()
     729                y = B.random_element()
     730        return par[0]([x,y])
     731
     732    def singular_point(self):
     733        r"""
     734        Returns a singular rational point of self
     735       
     736        EXAMPLES:
     737       
     738        ::
     739
     740            sage: Conic(GF(2), [1,1,1,1,1,1]).singular_point()
     741            (1 : 1 : 1)
     742
     743        A ValueError is raised if the conic has no rational singular point
     744           
     745        ::
     746
     747            sage: Conic(QQ, [1,1,1,1,1,1]).singular_point()
     748            Traceback (most recent call last):
     749            ...
     750            ValueError: The conic self (= Projective Conic Curve over Rational Field defined by x^2 + x*y + y^2 + x*z + y*z + z^2) has no rational singular point
     751        """
     752        b = self.has_singular_point(point = True)
     753        if not b[0]:
     754            raise ValueError, "The conic self (= %s) has no rational " \
     755                              "singular point" % self
     756        return b[1]
     757
     758    def symmetric_matrix(self):
     759        r"""
     760        The symmetric matrix M such that (x y z) M (x y z)^t
     761        is the defining equation of  self.
     762       
     763        EXAMPLES:
     764            sage: R.<x, y, z> = QQ[]
     765            sage: C = Curve(x^2 + x*y/2 + y^2 + z^2)
     766            sage: C.symmetric_matrix()
     767            [  1 1/4   0]
     768            [1/4   1   0]
     769            [  0   0   1]
     770
     771            sage: C = Curve(x^2 + 2*x*y + y^2 + 3*x*z + z^2)
     772            sage: v = vector([x, y, z])
     773            sage: v * C.symmetric_matrix() * v
     774            x^2 + 2*x*y + y^2 + 3*x*z + z^2
     775        """
     776        [a,b,c,d,e,f] = self.coefficients()
     777        if self.base_ring().characteristic() == 2:
     778            if b == 0 and c == 0 and e == 0:
     779                return matrix([[a,0,0],[0,d,0],[0,0,f]])
     780            raise ValueError, "The conic self (= %s) has no symmetric matrix " \
     781                              "because the base field has characteristic 2" % \
     782                              self
     783        from sage.matrix.constructor import matrix
     784        return matrix([[  a , b/2, c/2 ],
     785                       [ b/2,  d , e/2 ],
     786                       [ c/2, e/2,  f  ]])
     787 
     788    def upper_triangular_matrix(self):
     789        r"""
     790        The upper-triangular matrix M such that (x y z) M (x y z)^t
     791        is the defining equation of self.
     792
     793        EXAMPLES:
     794            sage: R.<x, y, z> = QQ[]
     795            sage: C = Curve(x^2 + x*y + y^2 + z^2)
     796            sage: C.upper_triangular_matrix()
     797            [1 1 0]
     798            [0 1 0]
     799            [0 0 1]
     800
     801            sage: C = Curve(x^2 + 2*x*y + y^2 + 3*x*z + z^2)
     802            sage: v = vector([x, y, z])
     803            sage: v * C.upper_triangular_matrix() * v
     804            x^2 + 2*x*y + y^2 + 3*x*z + z^2
     805        """
     806        from sage.matrix.constructor import matrix
     807        [a,b,c,d,e,f] = self.coefficients()
     808        return matrix([[ a, b, c ],
     809                       [ 0, d, e ],
     810                       [ 0, 0, f ]])
     811   
     812    def variable_names(self):
     813        r"""
     814        Returns the variable names of the defining polynomial
     815        of self
     816
     817        EXAMPLES:
     818
     819        ::
     820
     821            sage: c=Conic([1,1,0,1,0,1], 'x,y,z')
     822            sage: c.variable_names()
     823            ('x', 'y', 'z')
     824            sage: c.variable_name()
     825            'x'
     826        """
     827        return self.defining_polynomial().parent().variable_names()
     828   
     829
  • new file sage/schemes/plane_conics/rnfisnorm.py

    diff -r 8dec8b43ccca -r 3fcab21a698f sage/schemes/plane_conics/rnfisnorm.py
    - +  
     1r"""
     2Solve relative quadratic norm equations using PARI.
     3
     4The function rnfisnorm is a wrapper for the PARI function with the same name.
     5On input a relative quadratic number field `L` and an element `x` of its
     6base field `K`, the function rnfisnorm(L, x) determines whether `x` is the
     7relative norm `N_{L/K}(y)` of an element `y\in L`.
     8
     9AUTHORS:
     10
     11- Marco Streng (2009-08-07)
     12
     13"""
     14#*****************************************************************************
     15#       Copyright (C) 2009 Marco Streng <marco.streng@gmail.com>
     16#
     17#  Distributed under the terms of the GNU General Public License (GPL)
     18#
     19#    This code is distributed in the hope that it will be useful,
     20#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     21#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     22#    General Public License for more details.
     23#
     24#  The full text of the GPL is available at:
     25#
     26#                  http://www.gnu.org/licenses/
     27#*****************************************************************************
     28
     29from sage.libs.pari.gen import pari
     30from sage.rings.all import (NumberField, PolynomialRing, is_NumberField,
     31                            is_NumberFieldElement, is_AbsoluteNumberField,
     32                            is_RationalField, ZZ, QQ)
     33from sage.structure.sequence import Sequence
     34
     35
     36
     37def rnfisnorm(L, x, galois_check = 2, extra_primes = 0):
     38    r"""
     39    For a relative number field `L` and an element `x` of its base field `K`,
     40    determine whether `x` is the relative norm `N_{L/K}(y)` of an element
     41    `y\in L`.
     42   
     43    The output takes the form (True, y) if the answer is positive, and
     44    (False, None) otherwise.
     45   
     46    The output is guaranteed if `L/K` is Galois. Otherwise, the parameters
     47    galois_check and extra_primes can be used as follows.
     48   
     49    galois_check can be 0, 1, or 2.
     50        0 do not care if L/K is Galois
     51        1 assume L/K is Galois
     52        2 let PARI determine whether L/K is Galois
     53   
     54    extra_primes:
     55        The norm equation is solved in S-units, where `S` is a set determined
     56        by the PARI. If extra_primes is a positive integer, add all primes
     57        <= extra_primes to `S`. If extra_primes is a negative integer, add
     58        all primes dividing extra_primes to `S`. If `S` contains all primes
     59        less than `12`log(disc(`M`)), where `M` is the normal closure of `L/K`,
     60        then the output is guaranteed under GRH.
     61   
     62    ALGORITHM:
     63   
     64        Uses PARI's rnfisnorm or bnfisnorm to find a solution in the `S`-unit
     65        group for a suitable set of primes `S` of `L`.
     66   
     67    EXAMPLES:
     68   
     69    Absolute number fields are interpreted as relative fields over `QQ`.
     70   
     71    ::
     72        sage: from sage.schemes.plane_conics.rnfisnorm import rnfisnorm
     73        sage: rnfisnorm(QuadraticField(-1, 'i'), 2)
     74        (True, i + 1)
     75        sage: rnfisnorm(CyclotomicField(7), 7)
     76        (True,  zeta7 - 1)
     77        sage: rnfisnorm(CyclotomicField(7), 3)
     78        (False, None)
     79
     80    A higher degree example:
     81   
     82    ::
     83   
     84        sage: from sage.schemes.plane_conics.rnfisnorm import rnfisnorm
     85        sage: P.<x> = QQ[]
     86        sage: K.<a> = NumberField(x^3 + x + 1)
     87        sage: Q.<X> = K[]
     88        sage: L.<b> = NumberField(X^4 + a)
     89        sage: t = rnfisnorm(L, -a); t
     90        (True, b^3 + 1)
     91        sage: t[1].norm(K) == -a
     92        True
     93        sage: rnfisnorm(K, 2)
     94        (False, None)
     95        sage: t = rnfisnorm(L, 2); t
     96        (True, -a*b^3 + b^2 + (a^2 + 1)*b + a^2 + 1)
     97        sage: t[1].norm(K) == 2
     98        True
     99
     100    A relative quadratic number field `L` over `K` can be given by an
     101    element `z` of `K` such that `L = K(\sqrt{z})`. In that case, the
     102    output `y \in L` is given by the tuple `(a, b)` such that
     103    `y = a + b\sqrt{z}`.
     104   
     105    ::
     106   
     107        sage: from sage.schemes.plane_conics.rnfisnorm import rnfisnorm
     108        sage: rnfisnorm(-1, 2)
     109        (True, (1, 1))
     110        sage: K.<b> = NumberField(x^5+x+3)
     111        sage: rnfisnorm(b, 1 - b)
     112        (True, (-1, 1))
     113        sage: rnfisnorm(b/2, -5)
     114        (False, None)
     115        sage: rnfisnorm(b, (1 - b)/3)
     116        (False, None)
     117    """
     118    if is_NumberFieldElement(L) or L in QQ:
     119        s = Sequence([L, x])
     120        K = s.universe()
     121        if K == ZZ:
     122            K = QQ
     123        elif not (K == QQ or is_AbsoluteNumberField(K)):
     124            raise ValueError, "L (=%s) and x (=%s) must lie in the same " \
     125                              "number field" % (L, x)
     126        X = PolynomialRing(K, 'X').gen()
     127        x = K(x)
     128        den = L.denominator()
     129        L = NumberField(X**2-L*den**2, K.variable_name() + 'L')
     130        ret = rnfisnorm(L, x, galois_check = galois_check,
     131                      extra_primes = extra_primes)
     132        if ret[0]:
     133            ret1 = ret[1].vector()
     134            return (True, (ret1[0], ret1[1]*den))
     135        return ret
     136    elif is_NumberField(L):
     137        K = L.base_field()
     138        x = K(x)
     139    else:
     140        raise TypeError, "L (=%s) must be a number field or an element of a " \
     141                         "number field" % L   
     142    if is_RationalField(K) or (is_AbsoluteNumberField(K) and K.degree() == 1):
     143        pol1 = 'y'
     144        pol2 = L.pari_polynomial()
     145        xpari = QQ(x)._pari_()
     146    elif is_AbsoluteNumberField(K):
     147        pol1 = K.pari_polynomial('y')
     148        pol2 = L.pari_relative_polynomial().lift()
     149        xpari = x._pari_('y')
     150    else:
     151        raise ValueError, "Base field of L (=%s) needs to be an absolute" \
     152                          "number field in rnfisnorm" % L
     153    T = pari("rnfisnorm(rnfisnorminit(%s, %s, %s), %s, %s)" %
     154              (pol1, pol2, galois_check, xpari, extra_primes))
     155    if T[1] == 1:
     156        Q = PolynomialRing(K, 'X')
     157        return True, L(Q(T[0].lift()))
     158    return False, None
     159   
     160   
     161   
  • setup.py

    diff -r 8dec8b43ccca -r 3fcab21a698f setup.py
    a b  
    907907                     'sage.schemes.generic',
    908908                     'sage.schemes.jacobians',
    909909                     'sage.schemes.plane_curves',
     910                     'sage.schemes.plane_conics',
    910911                     'sage.schemes.plane_quartics',
    911912                     'sage.schemes.elliptic_curves',
    912913                     'sage.schemes.hyperelliptic_curves',