Ticket #11685: 11685.patch

File 11685.patch, 17.5 KB (added by jdemeyer, 11 years ago)

New alternative patch, rewrites PARI->Integer coercion, fast and many features

  • sage/libs/pari/gen.pyx

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1313570400 -7200
    # Node ID 6167385ad10e180852f361a063812a172e588105
    # Parent  2a2abbcad325ccca9399981ceddf5897eb467e64
    Rewrite PARI -> Integer coercion
    Fix creation of zero elements for class FiniteField_ext_pari
    Re-enable use of PARI's exponentiation in FiniteField_ext_pari
    
    diff -r 2a2abbcad325 -r 6167385ad10e sage/libs/pari/gen.pyx
    a b  
    18321832            sage: a.gequal(c)
    18331833            False
    18341834        """
    1835         global t0
    18361835        t0GEN(b)
    18371836        sig_on()
    1838         ret = gequal(a.g, t0)
     1837        cdef int ret = gequal(a.g, t0)
     1838        sig_off()
     1839        return ret != 0
     1840
     1841    def gequal0(gen a):
     1842        r"""
     1843        Check whether `a` is equal to zero.
     1844       
     1845        EXAMPLES::
     1846           
     1847            sage: pari(0).gequal0()
     1848            True
     1849            sage: pari(1).gequal0()
     1850            False
     1851            sage: pari(1e-100).gequal0()
     1852            False
     1853            sage: pari("0.0 + 0.0*I").gequal0()
     1854            True
     1855            sage: pari(GF(3^20,'t')(0)).gequal0()
     1856            True
     1857        """
     1858        sig_on()
     1859        cdef int ret = gequal0(a.g)
    18391860        sig_off()
    18401861        return ret != 0
    18411862
     
    18621883            False
    18631884        """
    18641885        sig_on()
    1865         ret = gequalsg(b, a.g)
     1886        cdef int ret = gequalsg(b, a.g)
    18661887        sig_off()
    18671888        return ret != 0
    18681889   
  • sage/rings/finite_rings/element_ext_pari.py

    diff -r 2a2abbcad325 -r 6167385ad10e sage/rings/finite_rings/element_ext_pari.py
    a b  
    7979        sage: a - a == K(0)
    8080        True
    8181    """
    82     def __init__(self, parent, value, check=True):
     82    def __init__(self, parent, value, value_from_pari=False):
    8383        """
    8484        Create element of a finite field.
    8585       
     
    114114            2*t
    115115            sage: pari(x)
    116116            Mod(Mod(2, 5)*a, Mod(1, 5)*a^2 + Mod(3, 5))
     117
     118        When initializing from a list, the elements are first coerced
     119        to the prime field (#11685)::
     120
     121            sage: k = FiniteField_ext_pari(3^11, 't')
     122            sage: k([ 0, 1/2 ])
     123            2*t
     124            sage: k([ k(0), k(1) ])
     125            t
     126            sage: k([ GF(3)(2), GF(3^5,'u')(1) ])
     127            t + 2
     128            sage: R.<x> = PolynomialRing(k)
     129            sage: k([ R(-1), x/x ])
     130            t + 2
     131
     132        TESTS:
     133
     134        Check that zeros are created correctly (#11685)::
     135
     136            sage: from sage.rings.finite_rings.finite_field_ext_pari import FiniteField_ext_pari
     137            sage: K = FiniteField_ext_pari(3^11, 't'); a = K.0
     138            sage: v = 0; pari(K(v)).lift()
     139            Mod(0, 3)
     140            sage: v = Mod(0,3); pari(K(v)).lift()
     141            Mod(0, 3)
     142            sage: v = pari(0); pari(K(v)).lift()
     143            Mod(0, 3)
     144            sage: v = pari("Mod(0,3)"); pari(K(v)).lift()
     145            Mod(0, 3)
     146            sage: v = []; pari(K(v)).lift()
     147            Mod(0, 3)
     148            sage: v = [0]; pari(K(v)).lift()
     149            Mod(0, 3)
     150            sage: v = [0,0]; pari(K(v)).lift()
     151            Mod(0, 3)
     152            sage: v = pari("Pol(0)"); pari(K(v)).lift()
     153            Mod(0, 3)
     154            sage: v = pari("Mod(0, %s)"%K._pari_modulus()); pari(K(v)).lift()
     155            Mod(0, 3)
     156            sage: v = pari("Mod(Pol(0), %s)"%K._pari_modulus()); pari(K(v)).lift()
     157            Mod(0, 3)
     158            sage: v = K(1) - K(1); pari(K(v)).lift()
     159            Mod(0, 3)
     160            sage: v = K([1]) - K([1]); pari(K(v)).lift()
     161            Mod(0, 3)
     162            sage: v = a - a; pari(K(v)).lift()
     163            Mod(0, 3)
     164            sage: v = K(1)*0; pari(K(v)).lift()
     165            Mod(0, 3)
     166            sage: v = K([1])*K([0]); pari(K(v)).lift()
     167            Mod(0, 3)
     168            sage: v = a*0; pari(K(v)).lift()
     169            Mod(0, 3)
     170
    117171        """
    118172        field_element.FieldElement.__init__(self, parent)
    119173        self.__class__ = dynamic_FiniteField_ext_pariElement
    120         if isinstance(value, str):
    121             raise TypeError, "value must not be a string"
    122         if not check:
    123             if not value:  # see comment below about this workaround
    124                 self.__value = pari(0).Mod(parent._pari_modulus())*parent._pari_one()
    125             else:
    126                 self.__value = value
     174
     175        # If value_from_pari is True, directly set self.__value to value.
     176        # This assumes that value is a POLMOD with the correct modulus
     177        # whose lift is either an INTMOD or a POL with INTMOD
     178        # coefficients. In practice, this means that value comes from a
     179        # PARI calculation with other elements of this finite field.
     180        # This assumption is not checked.
     181        if value_from_pari:
     182            self.__value = value
    127183            return
     184
    128185        try:
    129186            if isinstance(value, pari_gen):
    130187                try:
    131188                    # In some cases we get two different versions of
    132189                    # the 0 element of a finite field.  This has to do
    133                     # with how PARI's Mod function works -- it treats
    134                     # 0 differently.  In particular, if value is a a
    135                     # pari t_POLMOD that is 0, modding it simply
    136                     # doesn't work correctly.  We fix this by changing
    137                     # the value in the 0 case to the standard pari 0,
    138                     # which works correctly.  Note -- probably all
    139                     # this code should be replaced by much faster
    140                     # finite field arithmetic programmed against NTL.
    141                     if not value:
    142                         value = pari(0)
    143                     if value.type() == "t_POLMOD":
    144                         self.__value = value * parent._pari_one()
    145                     else:
    146                         self.__value = value.Mod(parent._pari_modulus())*parent._pari_one()
     190                    # with how PARI's Mod() function works -- it treats
     191                    # polynomials different from integers.
     192                    # Also, we need to fix things like ``1/Mod(3,5)`` when
     193                    # we want ``Mod(2,5)``.
     194                    # These issues are solved by first simplifying the
     195                    # given value.
     196                    value = value.simplify()
     197                    if value.type() != "t_POLMOD":
     198                        value = value.Mod(parent._pari_modulus())
     199                    self.__value = value * parent._pari_one()
    147200                except RuntimeError:
    148201                    raise TypeError, "no possible coercion implemented"
    149202            elif isinstance(value, FiniteField_ext_pariElement):
     
    158211                for i in range(len(value)):
    159212                    self.__value = self.__value + pari(int(value[i])).Mod(parent._pari_modulus())*pari("a^%s"%i)
    160213            elif isinstance(value, list):
    161                 # Make the list into a polynomial with variable "a",
    162                 # multiply it with _pari_one (which equals Mod(1,p)) to
    163                 # fix the characteristic and then mod out by the modulus.
    164                 # See #10486.
    165                 self.__value = (pari(value).Polrev("a") * parent._pari_one()).Mod(parent._pari_modulus())
     214                # AUTHOR: Jeroen Demeyer
     215                # See Trac #10486 and #11685 for comments about this
     216                # implementation
     217
     218                # We need a special case for the empty list because
     219                # PARI's ``Pol([])`` returns an exact 0 while we want
     220                # ``Mod(0,3)``.
     221                if not value:
     222                    value = [0]
     223                try:
     224                    # First, try the conversion directly in PARI.  This
     225                    # should cover the most common cases, like converting
     226                    # from integers or intmods.
     227                    # Convert the list to PARI, then mod out the
     228                    # characteristic (PARI can do this directly for lists),
     229                    # convert to a polynomial with variable "a" and finally
     230                    # mod out the field modulus.
     231                    self.__value = pari(value).Mod(parent.characteristic()).Polrev("a").Mod(parent._pari_modulus())
     232                except RuntimeError:
     233                    # That didn't work, do it in a more general but also
     234                    # slower way: first convert all list elements to the
     235                    # prime field.
     236                    GFp = parent.prime_subfield()
     237                    self.__value = pari([GFp(c) for c in value]).Polrev("a").Mod(parent._pari_modulus())
     238            elif isinstance(value, str):
     239                raise TypeError, "value must not be a string"
    166240            else:
    167241                try:
    168242                    self.__value = pari(value).Mod(parent._pari_modulus())*parent._pari_one()
     
    378452            sage: a is a
    379453            True
    380454        """
    381         return FiniteField_ext_pariElement(self.parent(), self.__value, check=False)
     455        return FiniteField_ext_pariElement(self.parent(), self.__value, value_from_pari=True)
    382456
    383457    def _pari_(self, var=None):
    384458        """
     
    489563            raise TypeError, "Parents of finite field elements must be equal."
    490564       
    491565    def _add_(self, right):
    492         return FiniteField_ext_pariElement(self.parent(), self.__value + right.__value, check=False)
     566        return FiniteField_ext_pariElement(self.parent(), self.__value + right.__value, value_from_pari=True)
    493567   
    494568    def _sub_(self, right):
    495         return FiniteField_ext_pariElement(self.parent(), self.__value - right.__value, check=False)
     569        return FiniteField_ext_pariElement(self.parent(), self.__value - right.__value, value_from_pari=True)
    496570   
    497571    def _mul_(self, right):
    498         return FiniteField_ext_pariElement(self.parent(), self.__value * right.__value, check=False)
     572        return FiniteField_ext_pariElement(self.parent(), self.__value * right.__value, value_from_pari=True)
    499573
    500574    def _div_(self, right):
    501575        if right.__value == 0:
    502576            raise ZeroDivisionError
    503         return FiniteField_ext_pariElement(self.parent(), self.__value / right.__value, check=False)
     577        return FiniteField_ext_pariElement(self.parent(), self.__value / right.__value, value_from_pari=True)
    504578
    505579    def __int__(self):
    506580        try:
     
    523597        except ValueError:
    524598            raise TypeError, "cannot coerce to float"
    525599
    526     # commented out because PARI (used for .__value) prints
    527     # out crazy warnings when the exponent is LARGE -- this
    528     # is even a problem in gp!!!
    529     # (Commenting out causes this to use a generic algorithm)
    530     #def __pow__(self, _right):
    531     #    right = int(_right)
    532     #    if right != _right:
    533     #         raise ValueError
    534     #    return FiniteField_ext_pariElement(self.parent(), self.__value**right, check=False)
     600    def __pow__(self, _right):
     601        """
     602        TESTS::
     603
     604            sage: K.<a> = GF(5^10)
     605            sage: n = (2*a)/a
     606
     607        Naively compute n^-15 in PARI, note that the result is `1/3`.
     608        This is mathematically correct (modulo 5), but not what we want.
     609        In particular, comparisons will fail::
     610
     611            sage: pari(n)^-15
     612            Mod(1/Mod(3, 5), Mod(1, 5)*a^10 + Mod(3, 5)*a^5 + Mod(3, 5)*a^4 + Mod(2, 5)*a^3 + Mod(4, 5)*a^2 + Mod(1, 5)*a + Mod(2, 5))
     613
     614        We need to ``simplify()`` the result (which is done in the
     615        ``FiniteField_ext_pariElement()`` constructor::
     616
     617            sage: n^-15
     618            2
     619
     620        Large exponents are not a problem::
     621
     622            sage: e = 3^10000
     623            sage: a^e
     624            2*a^9 + a^5 + 4*a^4 + 4*a^3 + a^2 + 3*a
     625            sage: a^(e % (5^10 - 1))
     626            2*a^9 + a^5 + 4*a^4 + 4*a^3 + a^2 + 3*a
     627        """
     628        right = int(_right)
     629        if right != _right:
     630             raise ValueError
     631        # It is important to set value_from_pari=False, see doctest above!
     632        return FiniteField_ext_pariElement(self.parent(), self.__value**right, value_from_pari=False)
    535633
    536634    def __neg__(self):
    537         return FiniteField_ext_pariElement(self.parent(), -self.__value, check=False)
     635        return FiniteField_ext_pariElement(self.parent(), -self.__value, value_from_pari=True)
    538636   
    539637    def __pos__(self):
    540638        return self
     
    552650            a + 2
    553651            sage: (a+1)*a
    554652            2*a + 1
     653            sage: ~((2*a)/a)
     654            2
    555655        """
    556656
    557657        if self.__value == 0:
    558658            raise ZeroDivisionError, "Cannot invert 0"
    559         return FiniteField_ext_pariElement(self.parent(), ~self.__value, check=False)
     659        return FiniteField_ext_pariElement(self.parent(), ~self.__value, value_from_pari=True)
    560660   
    561661    def lift(self):
    562662        """
  • sage/rings/integer.pyx

    diff -r 2a2abbcad325 -r 6167385ad10e sage/rings/integer.pyx
    a b  
    528528            2
    529529       
    530530        ::
    531        
    532             sage: t = pari(0*ZZ[x].0 + 3)
    533             sage: t.type()
    534             't_POL'
    535             sage: ZZ(t)
    536             3
    537531
    538532            sage: ZZ(float(2.0))
    539533            2
     
    560554            4
    561555            sage: ZZ(MyFloat(5))
    562556            5
     557
     558        Test conversion from PARI (#11685)::
     559       
     560            sage: ZZ(pari(-3))
     561            -3
     562            sage: ZZ(pari("-3.0"))
     563            -3
     564            sage: ZZ(pari("-3.5"))
     565            Traceback (most recent call last):
     566            ...
     567            TypeError: Attempt to coerce non-integral real number to an Integer
     568            sage: ZZ(pari("1e100"))
     569            Traceback (most recent call last):
     570            ...
     571            PariError: precision too low (10)
     572            sage: ZZ(pari("10^50"))
     573            100000000000000000000000000000000000000000000000000
     574            sage: ZZ(pari("Pol(3)"))
     575            3
     576            sage: ZZ(GF(3^20,'t')(1))
     577            1
     578            sage: ZZ(pari(GF(3^20,'t')(1)))
     579            1
     580            sage: x = polygen(QQ)
     581            sage: K.<a> = NumberField(x^2+3)
     582            sage: ZZ(a^2)
     583            -3
     584            sage: ZZ(pari(a)^2)
     585            -3
     586            sage: ZZ(pari("Mod(x, x^3+x+1)"))   # Note error message refers to lifted element
     587            Traceback (most recent call last):
     588            ...
     589            TypeError: Unable to coerce PARI x to an Integer
    563590        """
    564591
    565592        # TODO: All the code below should somehow be in an external
     
    574601
    575602        cdef Integer tmp
    576603        cdef char* xs
     604        cdef int paritype
    577605       
    578606        cdef Element lift
    579607       
     
    604632                    raise TypeError, "Cannot convert non-integral float to integer"
    605633
    606634            elif PY_TYPE_CHECK(x, pari_gen):
    607                
    608                 if typ((<pari_gen>x).g) == t_INT:
    609                     t_INT_to_ZZ(self.value, (<pari_gen>x).g)
    610                    
    611                 else:
    612                     if typ((<pari_gen>x).g) == t_INTMOD:
     635                # Simplify and lift until we get an integer
     636                while typ((<pari_gen>x).g) != t_INT:
     637                    x = x.simplify()
     638                    paritype = typ((<pari_gen>x).g)
     639                    if paritype == t_INT:
     640                        break
     641                    elif paritype == t_REAL:
     642                        # Check that the fractional part is zero
     643                        if not x.frac().gequal0():
     644                            raise TypeError, "Attempt to coerce non-integral real number to an Integer"
     645                        # floor yields an integer
     646                        x = x.floor()
     647                        break
     648                    elif paritype == t_PADIC:
     649                        # Lifting a PADIC yields an integer
    613650                        x = x.lift()
    614                     # TODO: figure out how to convert to pari integer in base 16 ?
    615                    
    616                     # todo: having this "s" variable around here is causing
    617                     # Cython to play games with refcount for the None object, which
    618                     # seems really stupid.
    619                
    620                     try:
    621                         s = hex(x)
    622                         base = 16
    623                     except:
    624                         s = str(x)
    625                         base = 10
    626                        
    627                     if mpz_set_str(self.value, s, base) != 0:
    628                         raise TypeError, "Unable to coerce PARI %s to an Integer."%x
     651                        break
     652                    elif paritype == t_INTMOD:
     653                        # Lifting an INTMOD yields an integer
     654                        x = x.lift()
     655                        break
     656                    elif paritype == t_POLMOD:
     657                        x = x.lift()
     658                    else:
     659                        raise TypeError, "Unable to coerce PARI %s to an Integer"%x
     660
     661                # Now we have a true PARI integer, convert it to Sage
     662                t_INT_to_ZZ(self.value, (<pari_gen>x).g)
    629663
    630664            elif PyString_Check(x):
    631665                if base < 0 or base > 36:
  • sage/rings/rational.pyx

    diff -r 2a2abbcad325 -r 6167385ad10e sage/rings/rational.pyx
    a b  
    429429            sage: Rational(pari('x'))
    430430            Traceback (most recent call last):
    431431            ...
    432             TypeError: Unable to coerce PARI x to an Integer.
     432            TypeError: Unable to coerce PARI x to an Integer
    433433        """
    434434        mpq_set_str(self.value, s, 32)
    435435