Ticket #727: trac_727-conic.patch

File trac_727-conic.patch, 17.5 KB (added by Marco Streng, 13 years ago)

Define a Conic class with an interface to Denis Simon's qfsolve

  • sage/schemes/plane_curves/constructor.py

    # HG changeset patch
    # User Marco Streng <marco.streng@gmail.com>
    # Date 1245316852 -3600
    # Node ID 0c3d0a49c689e35972e9ab37f1a3b354713236ce
    # Parent  5f46e0147001357af9779b11e9658b6c53fd6caa
    Integrate code of Denis Simon for finding points on rational conics.
    
    diff -r 5f46e0147001 -r 0c3d0a49c689 sage/schemes/plane_curves/constructor.py
    a b  
    2323#                  http://www.gnu.org/licenses/
    2424#*****************************************************************************
    2525
    26 from sage.rings.all import is_MPolynomial, is_MPolynomialRing, is_FiniteField
    27 
     26from sage.rings.all import PolynomialRing, is_MPolynomial, is_MPolynomialRing, is_FiniteField
    2827from sage.structure.all import Sequence
    29 
     28from sage.structure.element import is_Matrix
    3029from sage.schemes.generic.all import (is_AmbientSpace, is_AlgebraicScheme,
    3130                                      AffineSpace, ProjectiveSpace)
     31from sage.modules.free_module_element import vector
    3232
    3333import affine_curve
    3434
     
    3636                              ProjectiveSpaceCurve_generic,
    3737                              ProjectiveCurve_finite_field,
    3838                              ProjectiveCurve_prime_finite_field)
     39from projective_conic import (ProjectiveConic_generic)
    3940
    4041from affine_curve import (AffineCurve_generic,
    4142                          AffineSpaceCurve_generic,
     
    4445
    4546def Curve(F):
    4647    """
    47     Return the plane or space curve defined by `F`, where
    48     `F` can be either a multivariate polynomial, a list or
    49     tuple of polynomials, or an algebraic scheme.
     48    Return the plane or space curve defined by `F`, where `F` can be
     49    either a multivariate polynomial, a list or tuple of polynomials,
     50    an algebraic scheme, or a square matrix.
    5051   
    51     If `F` is in two variables the curve is affine, and if it
    52     is homogenous in `3` variables, then the curve is
    53     projective.
     52    If `F` is in two variables the curve is affine, and if it is
     53    homogenous in `3` variables, then the curve is projective.
    5454   
    5555    EXAMPLE: A projective plane curve
    5656   
     
    101101   
    102102        sage: x,y,z = QQ['x,y,z'].gens()
    103103        sage: Curve((x-y)*(x+y))
    104         Projective Curve over Rational Field defined by x^2 - y^2
     104        Projective Conic Curve over Rational Field defined by x^2 - y^2
    105105        sage: Curve((x-y)^2*(x+y)^2)
    106106        Projective Curve over Rational Field defined by x^4 - 2*x^2*y^2 + y^4
    107107   
     
    144144   
    145145        sage: x,y,z = QQ['x,y,z'].gens()
    146146        sage: Curve(x^2+y^2)
    147         Projective Curve over Rational Field defined by x^2 + y^2
     147        Projective Conic Curve over Rational Field defined by x^2 + y^2
    148148        sage: Curve(x^2+y^2+z)
    149149        Traceback (most recent call last):
    150150        ...
     
    157157        Traceback (most recent call last):
    158158        ...
    159159        ValueError: defining polynomial of curve must be nonzero
     160
     161    Plane curves are often represented by matrices::
     162
     163        sage: Curve(matrix(QQ, [[2, 3], [4, 5]]))
     164        Affine Curve over Rational Field defined by 2*x^2 + 7*x*y + 5*y^2
     165
     166        sage: x,y,z = GF(11)['x,y,z'].gens()
     167        sage: C = Curve(x^2+y^2-2*z^2); C
     168        Projective Conic Curve over Finite Field of size 11 defined by x^2 + y^2 - 2*z^2
     169        sage: Curve(C.symmetric_matrix())
     170        Projective Conic Curve over Finite Field of size 11 defined by x^2 + y^2 - 2*z^2
    160171    """
    161172    if is_AlgebraicScheme(F):
    162173        return Curve(F.defining_polynomials())
    163 
    164174    if isinstance(F, (list, tuple)):
    165175        if len(F) == 1:
    166176            return Curve(F[0])
     
    178188        A = ProjectiveSpace(P.ngens()-1, P.base_ring())
    179189        A._coordinate_ring = P
    180190        return ProjectiveSpaceCurve_generic(A, F)
    181            
     191
     192    if is_Matrix(F) and F.is_square():
     193        if F.ncols() == 2:
     194            temp_ring = PolynomialRing(F.parent().base_ring(), F.ncols(), 'x, y')
     195        elif F.ncols() == 3:
     196            temp_ring = PolynomialRing(F.parent().base_ring(), F.ncols(), 'x, y, z')
     197        else:
     198            temp_ring = PolynomialRing(F.parent().base_ring(), F.ncols(), 'x')
     199        F = vector(temp_ring.gens()) * F * vector(temp_ring.gens())
     200
    182201    if not is_MPolynomial(F):
    183202        raise TypeError, "F (=%s) must be a multivariate polynomial"%F
     203    if F == 0:
     204        raise ValueError, "defining polynomial of curve must be nonzero"
    184205
    185206    P = F.parent()
    186207    k = F.base_ring()
    187208    if F.parent().ngens() == 2:
    188         if F == 0:
    189             raise ValueError, "defining polynomial of curve must be nonzero"
    190209        A2 = AffineSpace(2, P.base_ring())
    191210        A2._coordinate_ring = P
    192211
     
    195214                return AffineCurve_prime_finite_field(A2, F)
    196215            else:
    197216                return AffineCurve_finite_field(A2, F)
    198         else:
    199             return AffineCurve_generic(A2, F)
    200217
    201     elif F.parent().ngens() == 3:
    202         if F == 0:
    203             raise ValueError, "defining polynomial of curve must be nonzero"
     218        return AffineCurve_generic(A2, F)
     219
     220    if F.parent().ngens() == 3:
    204221        P2 = ProjectiveSpace(2, P.base_ring())
    205222        P2._coordinate_ring = P
     223
     224        if F.total_degree() == 2:
     225            return ProjectiveConic_generic(P2, F)
     226       
    206227        if is_FiniteField(k):
    207228            if k.is_prime_field():
    208229                return ProjectiveCurve_prime_finite_field(P2, F)
    209230            else:
    210231                return ProjectiveCurve_finite_field(P2, F)
    211         else:
    212             return ProjectiveCurve_generic(P2, F)
     232        return ProjectiveCurve_generic(P2, F)
    213233
    214 
    215     else:
    216 
    217         raise TypeError, "Number of variables of F (=%s) must be 2 or 3"%F
    218 
    219 
     234    raise TypeError, "Number of variables of F (=%s) must be 2 or 3" % F
  • new file sage/schemes/plane_curves/projective_conic.py

    diff -r 5f46e0147001 -r 0c3d0a49c689 sage/schemes/plane_curves/projective_conic.py
    - +  
     1r"""
     2Projective plane conics over a general ring.
     3
     4Let M be a 3-by-3 matrix.  We denote by X_M the projective variety defined by
     5(x, y, z) M (x, y, z)^t = 0.
     6
     7AUTHORS:
     8    -- 2008-01-08, Nick Alexander <ncalexander@gmail.com>
     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#  The full text of the GPL is available at:
     17#
     18#                  http://www.gnu.org/licenses/
     19#*****************************************************************************
     20
     21from sage.libs.pari.gen import pari
     22from sage.misc.all import add, sage_eval
     23from sage.rings.all import (PolynomialRing, ZZ, QQ, MPolynomialRing,
     24                            degree_lowest_rational_function, is_PrimeField)
     25from sage.modules.free_module_element import vector
     26from sage.structure.sequence import Sequence
     27from sage.schemes.generic.projective_space import is_ProjectiveSpace
     28
     29from curve import Curve_generic_projective
     30from projective_curve import ProjectiveCurve_generic
     31from qfsolve import Qfsolve, Qfparam
     32from wamelen import diagsymm
     33
     34class ProjectiveConic_generic(ProjectiveCurve_generic):
     35    def __init__(self, A, f):
     36        ProjectiveCurve_generic.__init__(self, A, f)
     37        # We cache some things
     38        self._diagonal_matrix = None
     39        self._diagonal_transform = None
     40        self._symmetric_matrix = None
     41        self._rational_point = None
     42        self._parametrization = None
     43
     44    def _repr_type(self):
     45        return "Projective Conic"
     46
     47    def symmetric_matrix(self):
     48        r"""A symmetric matrix M such that (x y z) M (x y z)^t = self.
     49
     50        That is, X_M has defining equation the same as self.
     51
     52        EXAMPLES:
     53            sage: R.<x, y, z> = QQ[]
     54            sage: C = Curve(x^2 + y^2 + z^2)
     55            sage: C.symmetric_matrix()
     56            [1 0 0]
     57            [0 1 0]
     58            [0 0 1]
     59
     60            sage: C = Curve(x^2 + 2*x*y + y^2 + 3*x*z + z^2)
     61            sage: v = vector([x, y, z])
     62            sage: v * C.symmetric_matrix() * v
     63            x^2 + 2*x*y + y^2 + 3*x*z + z^2
     64
     65        TESTS:
     66            sage: C = Curve(2*x^2 + 3*y^2 + 4*z^2)
     67            sage: C.symmetric_matrix()
     68            [2 0 0]
     69            [0 3 0]
     70            [0 0 4]
     71        """
     72        if self._symmetric_matrix is not None:
     73            return self._symmetric_matrix
     74        from sage.matrix.constructor import matrix
     75        b = self.defining_polynomial()
     76        M = matrix([[ b[(2,0,0)]  , b[(1,1,0)]/2, b[(1,0,1)]/2 ],
     77                    [ b[(1,1,0)]/2, b[(0,2,0)]  , b[(0,1,1)]/2 ],
     78                    [ b[(1,0,1)]/2, b[(0,1,1)]/2, b[(0,0,2)]   ]])
     79        self._symmetric_matrix = M
     80        return self._symmetric_matrix
     81
     82    def normalized_symmetric_matrix(self):
     83        # Qfsolve barfs on matrices that have non-trivial denominator content
     84        from sage.rings.arith import lcm
     85        from constructor import Curve
     86        M = self.symmetric_matrix()
     87        M *= lcm([ t.denominator() for t in M.list() ])
     88        return M
     89
     90        # F = self.defining_polynomial()
     91        # return Curve(F * lcm([ t.denominator() for t in F.coefficients() ])).symmetric_matrix()
     92
     93    def determinant(self):
     94        r""" XXX
     95        """
     96
     97        return self.symmetric_matrix().determinant()
     98
     99    def discriminant(self):
     100        r""" XXX
     101        """
     102        return self.defining_polynomial().discriminant()
     103
     104    def diagonal_matrix(self):
     105        r"""A diagonal matrix D and a matrix Phi such that self and X_D are isomorphic
     106        via Phi.
     107
     108        EXAMPLES:
     109            sage: R.<x,y,z> = QQ[]
     110            sage: C = Curve(12*x^2 - 2*y^2 + 5*x*z + 5*z^2)
     111            sage: D, Phi = C.diagonal_matrix(); D
     112            [    12      0      0]
     113            [     0     -2      0]
     114            [     0      0 215/48]
     115           
     116            Let's test the isomorphism.  First, at a single point:
     117
     118            sage: P = C.rational_point(); P
     119            (1/3 : -2 : 1)
     120            sage: Q = vector(P) * Phi.inverse(); Q
     121            (13/24, -2, 1)
     122            sage: Q * D * Q
     123            0
     124
     125            Second, on a parametrization:
     126
     127            sage: X = Curve(D); X
     128            Projective Conic Curve over Rational Field defined by 12*x^2 - 2*y^2 + 215/48*z^2
     129            sage: param = C.parametrization(); param
     130            (3*t^2 - 14*t + 13, -7*t^2 + 31*t - 42, -2*t^2 + 12)
     131            sage: X.defining_polynomial()(*param * Phi.inverse())
     132            0
     133        """
     134        if self._diagonal_matrix is not None:
     135            return (self._diagonal_matrix, self._diagonal_transform)
     136        M = self.symmetric_matrix()
     137        P = diagsymm(M)
     138        self._diagonal_transform = P
     139        self._diagonal_matrix = P * M * P.transpose()
     140        return (self._diagonal_matrix, self._diagonal_transform)
     141
     142    standard_matrix = diagonal_matrix
     143
     144    def rational_point(self):
     145        r"""Return a rational point (x0, y0, z0) on self.
     146
     147        Raises ValueError if no rational point exists.  XXX Fixme
     148
     149        ALGORITHM:
     150            Uses Denis Simon's Qfsolve.
     151
     152        EXAMPLES:
     153            sage: R.<x,y,z> = QQ[]
     154            sage: C = Curve(7*x^2 + 2*y*z + z^2)
     155            sage: C.rational_point()
     156            (0 : 1 : 0)
     157
     158            sage: C = Curve(x^2 + 2*y^2 + z^2)
     159            sage: C.rational_point()
     160            Traceback (most recent call last):
     161            ...
     162            ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + 2*y^2 + z^2 has no rational points over Rational Field!
     163
     164            sage: C = Curve(x^2 + y^2 + 7*z^2)
     165            sage: C.rational_point()
     166            Traceback (most recent call last):
     167            ...
     168            ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + y^2 + 7*z^2 has no rational points over Rational Field!
     169        """
     170        if self._rational_point is not None:
     171            return self._rational_point
     172        pt = Qfsolve(self.normalized_symmetric_matrix())
     173        if pt in ZZ:
     174            raise ValueError, "Conic %s has no rational points over %s!" % (self, self.ambient_space().base_ring())
     175        self._rational_point = self.ambient_space()(pt[0], pt[1], pt[2])
     176        return self._rational_point
     177
     178    def parametrization(self):
     179        r"""Return a rational point (x0, y0, z0) on self.
     180
     181        Raises ValueError if no rational point exists.  XXX Fixme
     182
     183        ALGORITHM:
     184            Uses Denis Simon's Qfsolve and Qfparam.
     185
     186        EXAMPLES:
     187            sage: R.<x,y,z> = QQ[]
     188            sage: C = Curve(7*x^2 + 2*y*z + z^2)
     189            sage: P = C.parametrization(); P
     190            (-2*t, 7*t^2 + 1, -2)
     191            sage: C.defining_polynomial()(*P)
     192            0
     193
     194            sage: C = Curve(x^2 + 2*y^2 + z^2)
     195            sage: C.parametrization()
     196            Traceback (most recent call last):
     197            ...
     198            ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + 2*y^2 + z^2 has no rational points over Rational Field!
     199
     200            sage: C = Curve(x^2 + y^2 + 7*z^2)
     201            sage: C.parametrization()
     202            Traceback (most recent call last):
     203            ...
     204            ValueError: Conic Projective Conic Curve over Rational Field defined by x^2 + y^2 + 7*z^2 has no rational points over Rational Field!
     205        """
     206        if self._parametrization is not None:
     207            return self._parametrization
     208        return vector(list(Qfparam(self.symmetric_matrix(), list(self.rational_point()))))
     209    param = parametrization
  • new file sage/schemes/plane_curves/qfsolve.py

    diff -r 5f46e0147001 -r 0c3d0a49c689 sage/schemes/plane_curves/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