Ticket #7585: 7585_FpT.patch

File 7585_FpT.patch, 19.4 KB (added by Robert Bradshaw, 13 years ago)
  • module_list.py

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1259780572 28800
    # Node ID ed7c7bc46a5b14ae3cf44feddc5d18638b3ce095
    # Parent  380d1bf0c2d5e07b6bab0e3d0df3f6512248af8d
    Fast arithmatic in F_p(t).
    
    diff -r 380d1bf0c2d5 -r ed7c7bc46a5b module_list.py
    a b  
    10111011    Extension('sage.rings.fraction_field_element',
    10121012              sources = ['sage/rings/fraction_field_element.pyx']),
    10131013   
     1014    Extension('sage.rings.fraction_field_FpT',
     1015              sources = ['sage/rings/fraction_field_FpT.pyx'],
     1016              libraries = ["csage", "flint", "gmp", "gmpxx", "ntl", "zn_poly"],
     1017              extra_compile_args=["-std=c99", "-D_XPG6"],
     1018              include_dirs = [SAGE_ROOT+'/local/include/FLINT/'],
     1019              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
     1020   
    10141021    Extension('sage.rings.integer_mod',
    10151022              sources = ['sage/rings/integer_mod.pyx'],
    10161023              libraries = ['gmp']),
  • new file sage/rings/fraction_field_FpT.pyx

    diff -r 380d1bf0c2d5 -r ed7c7bc46a5b sage/rings/fraction_field_FpT.pyx
    - +  
     1import sys
     2
     3from sage.rings.all import GF
     4from sage.rings.ring cimport Ring
     5from sage.libs.flint.zmod_poly cimport *
     6from sage.structure.element cimport Element, ModuleElement, RingElement
     7import sage.algebras.algebra
     8
     9from sage.rings.integer_mod cimport jacobi_int, mod_inverse_int
     10
     11class Fpt(Ring):
     12    def __init__(self, long p, names='t'):
     13        assert 2 < p < 2**16
     14        if isinstance(names, tuple):
     15            names, = names
     16        self.p = p
     17        base_ring = GF(p)
     18        self.poly_ring = base_ring[names]
     19        sage.algebras.algebra.Algebra.__init__(self, base_ring, names=names, normalize=True)
     20        self._populate_coercion_lists_(element_constructor=FptElement)
     21   
     22#    def _element_constructor_(self, *args):
     23#        return FptElement(self, *args)
     24   
     25    def ngens(self):
     26        return 1
     27   
     28    def gen(self, ix):
     29        assert ix == 0
     30        return FptElement(self, self.poly_ring.gen())
     31   
     32    def _coerce_map_from_(self, R):
     33        return self.poly_ring.has_coerce_map_from(R)
     34   
     35    def __iter__(self):
     36        return self.iter()
     37   
     38    def iter(self, bound=None, start=None):
     39        """
     40            sage: from sage.rings.fraction_field_FpT import *
     41            sage: R.<t> = Fpt(5)
     42            sage: list(R.iter(2))
     43        """
     44        return Fpt_iter(self, bound, start)
     45
     46cdef class FptElement(RingElement):
     47
     48    cdef long p
     49    cdef zmod_poly_t _numer, _denom
     50    cdef bint initalized
     51
     52    def __init__(self, parent, numer, denom=1):
     53        """
     54        EXAMPLES::
     55       
     56            sage: from sage.rings.fraction_field_FpT import *
     57            sage: R.<t> = Fpt(5)
     58            sage: R(7)
     59            2
     60           
     61        """
     62        RingElement.__init__(self, parent)
     63        numer = parent.poly_ring(numer)
     64        denom = parent.poly_ring(denom)
     65        self.p = parent.p
     66        zmod_poly_init(self._numer, self.p)
     67        zmod_poly_init(self._denom, self.p)
     68        self.initalized = True
     69        cdef long n
     70        for n, a in enumerate(numer):
     71            zmod_poly_set_coeff_ui(self._numer, n, a)
     72        for n, a in enumerate(denom):
     73            zmod_poly_set_coeff_ui(self._denom, n, a)
     74        normalize(self._numer, self._denom, self.p)
     75
     76    def __dealloc__(self):
     77        if self.initalized:
     78            zmod_poly_clear(self._numer)
     79            zmod_poly_clear(self._denom)
     80   
     81    def _new_c(self):
     82        cdef FptElement x = <FptElement>PY_NEW(FptElement)
     83        x._parent = self._parent
     84        x.p = self.p
     85        zmod_poly_init(x._numer, x.p)
     86        zmod_poly_init(x._denom, x.p)
     87        self.initalized = True
     88        return x
     89   
     90    def numer(self):
     91        cdef long n
     92        return self._parent.poly_ring(
     93            [zmod_poly_get_coeff_ui (self._numer, n) for n in range(0, zmod_poly_degree(self._numer)+1)])
     94
     95    def denom(self):
     96        cdef long n
     97        return self._parent.poly_ring(
     98            [zmod_poly_get_coeff_ui (self._denom, n) for n in range(0, zmod_poly_degree(self._denom)+1)])
     99
     100    def _repr_(self):
     101        """
     102            sage: from sage.rings.fraction_field_FpT import *
     103            sage: R.<t> = Fpt(17)
     104            sage: -t
     105            16*t
     106            sage: 1/t
     107            1/t
     108            sage: 1/(t+1)
     109            1/(t + 1)
     110            sage: 1-t/t
     111            0
     112            sage: (1-t)/t
     113            (16*t + 1)/t
     114        """
     115        if zmod_poly_degree(self._denom) == 0 and zmod_poly_get_coeff_ui(self._denom, 0) == 1:
     116            return repr(self.numer())
     117        else:
     118            numer_s = repr(self.numer())
     119            denom_s = repr(self.denom())
     120            if '-' in numer_s or '+' in numer_s:
     121                numer_s = "(%s)" % numer_s
     122            if '-' in denom_s or '+' in denom_s:
     123                denom_s = "(%s)" % denom_s
     124            return "%s/%s" % (numer_s, denom_s)
     125   
     126    def _latex_(self):
     127        if zmod_poly_degree(self._denom) == 0 and zmod_poly_get_coeff_ui(self._denom, 0) == 1:
     128            return self.numer()._latex_()
     129        else:
     130            return "{%s}/{%s}" % (self.numer()._latex_(), self.denom()._latex_())
     131        return
     132   
     133    def __cmp__(self, other):
     134        """
     135        TESTS::
     136       
     137            sage: from sage.rings.fraction_field_FpT import *
     138            sage: R.<t> = Fpt(7)
     139            sage: t == t
     140            True
     141            sage: t == -t
     142            False
     143            sage: -t == 6*t
     144            True
     145            sage: 1/t == 1/t
     146            True
     147            sage: 1/t == 1/(t+1)
     148            False
     149        """
     150        if isinstance(other, FptElement):
     151            # They are normalized.
     152            if (zmod_poly_equal(self._numer, (<FptElement>other)._numer) and
     153                    zmod_poly_equal(self._denom, (<FptElement>other)._denom)):
     154                return 0
     155            else:
     156                return -1
     157        else:
     158            return cmp(type(self), type(other))
     159       
     160    def __neg__(self):
     161        cdef FptElement x = self._new_c()
     162        zmod_poly_neg(x._numer, self._numer)
     163        zmod_poly_set(x._denom, self._denom)
     164        return x
     165       
     166    def __invert__(self):
     167        if zmod_poly_degree(self._numer) == -1:
     168            raise ZeroDivisionError
     169        cdef FptElement x = self._new_c()
     170        zmod_poly_set(x._denom, self._numer)
     171        zmod_poly_set(x._numer, self._denom)
     172        return x
     173
     174    cpdef ModuleElement _add_(self, ModuleElement _other):
     175        """
     176        EXAMPLES::
     177       
     178            sage: from sage.rings.fraction_field_FpT import *
     179            sage: R.<t> = Fpt(7)
     180            sage: t + t
     181            2*t
     182            sage: (t + 3) + (t + 10)
     183            2*t + 6
     184            sage: sum([t] * 7)
     185            0
     186            sage: 1/t + t
     187            (t^2 + 1)/(t)
     188            sage: 1/t + 1/t^2
     189            (t + 1)/(t^2)
     190        """
     191        cdef FptElement other = <FptElement>_other
     192        cdef FptElement x = self._new_c()
     193        zmod_poly_mul(x._numer, self._numer, other._denom)
     194        zmod_poly_mul(x._denom, self._denom, other._numer) # use x._denom as a temp
     195        zmod_poly_add(x._numer, x._numer, x._denom)
     196        zmod_poly_mul(x._denom, self._denom, other._denom)
     197        normalize(x._numer, x._denom, self.p)
     198        return x
     199
     200    cpdef ModuleElement _sub_(self, ModuleElement _other):
     201        """
     202        EXAMPLES::
     203       
     204            sage: from sage.rings.fraction_field_FpT import *
     205            sage: R.<t> = Fpt(7)
     206            sage: t - t
     207            0
     208            sage: (t + 3) - (t + 11)
     209            sage: 6
     210        """
     211        cdef FptElement other = <FptElement>_other
     212        cdef FptElement x = self._new_c()
     213        zmod_poly_mul(x._numer, self._numer, other._denom)
     214        zmod_poly_mul(x._denom, self._denom, other._numer) # use x._denom as a temp
     215        zmod_poly_sub(x._numer, x._numer, x._denom)
     216        zmod_poly_mul(x._denom, self._denom, other._denom)
     217        normalize(x._numer, x._denom, self.p)
     218        return x
     219
     220    cpdef RingElement _mul_(self, RingElement _other):
     221        """
     222        EXAMPLES::
     223       
     224            sage: from sage.rings.fraction_field_FpT import *
     225            sage: R.<t> = Fpt(7)
     226            sage: t * t
     227            t^2
     228            sage: (t + 3) * (t + 10)
     229            t^2 + 6*t + 2
     230        """
     231        cdef FptElement other = <FptElement>_other
     232        cdef FptElement x = self._new_c()
     233        zmod_poly_mul(x._numer, self._numer, other._numer)
     234        zmod_poly_mul(x._denom, self._denom, other._denom)
     235        normalize(x._numer, x._denom, self.p)
     236        return x
     237
     238    cpdef RingElement _div_(self, RingElement _other):
     239        """
     240        EXAMPLES::
     241       
     242            sage: from sage.rings.fraction_field_FpT import *
     243            sage: R.<t> = Fpt(5)
     244            sage: t / t
     245            1
     246            sage: (t + 3) / (t + 11)
     247            (t + 3)/(t + 1)
     248            sage: (t^2 + 2*t + 1) / (t + 1)
     249            t + 1
     250        """
     251        cdef FptElement other = <FptElement>_other
     252        if zmod_poly_degree(other._numer) == -1:
     253            raise ZeroDivisionError
     254        cdef FptElement x = self._new_c()
     255        zmod_poly_mul(x._numer, self._numer, other._denom)
     256        zmod_poly_mul(x._denom, self._denom, other._numer)
     257        normalize(x._numer, x._denom, self.p)
     258        return x
     259   
     260    cpdef FptElement next(self):
     261        """
     262        EXAMPLES::
     263       
     264            sage: from sage.rings.fraction_field_FpT import *
     265            sage: R.<t> = Fpt(3)
     266            sage: a = R(0)
     267            sage: for _ in range(15):
     268            ...       a = a.next()
     269            ...       print a
     270            ...   
     271            1
     272            2
     273            t
     274            t + 1
     275            (t + 1)/(t)
     276            t + 2
     277            (t + 2)/(t)
     278            (t + 2)/(t + 1)
     279            2*t
     280            (2*t)/(t + 1)
     281            (2*t)/(t + 2)
     282            2*t + 1
     283            (2*t + 1)/(t)
     284            (2*t + 1)/(t + 1)
     285            2*t + 2
     286        """
     287        cdef FptElement x = self._new_c()
     288        if zmod_poly_degree(self._numer) == -1:
     289            zmod_poly_set_coeff_ui(x._numer, 0, 1)
     290            zmod_poly_set_coeff_ui(x._denom, 0, 1)
     291            return x
     292        elif zmod_poly_degree(self._numer) == 0:
     293            zmod_poly_set_coeff_ui(x._numer, 0, zmod_poly_get_coeff_ui(self._numer, 0) + 1)
     294            zmod_poly_set_coeff_ui(x._denom, 0, 1)
     295            if zmod_poly_get_coeff_ui(x._numer, 0) == 0:
     296                zmod_poly_set_coeff_ui(x._numer, 1, 1)
     297            return x
     298        cdef zmod_poly_t g
     299        zmod_poly_init(g, self.p)
     300        zmod_poly_set(x._numer, self._numer)
     301        zmod_poly_set(x._denom, self._denom)
     302        while True:
     303            zmod_poly_inc(x._denom, True)
     304            if zmod_poly_equal(x._numer, x._denom) or zmod_poly_degree(x._numer) < zmod_poly_degree(x._denom):
     305                zmod_poly_inc(x._numer, False)
     306                zmod_poly_zero(x._denom)
     307                zmod_poly_set_coeff_ui(x._denom, 0, 1)
     308            zmod_poly_gcd(g, x._numer, x._denom)
     309            if zmod_poly_is_one(g):
     310                zmod_poly_clear(g)
     311                return x
     312   
     313    cpdef _sqrt_or_None(self):
     314        """
     315        Returns the squre root of self, or None. Differs from sqrt() by not raising an exception.
     316       
     317        TESTS::
     318       
     319            sage: from sage.rings.fraction_field_FpT import *
     320            sage: R.<t> = Fpt(17)
     321            sage: sqrt(t^2)
     322            t
     323            sage: sqrt(1/t^2)
     324            1/t
     325            sage: sqrt((1+t)^2)
     326            t + 1
     327            sage: sqrt((1+t)^2 / t^2)
     328            (t + 1)/t
     329
     330            sage: sqrt(4 * (1+t)^2 / t^2)
     331            (2*t + 2)/t
     332           
     333            sage: sqrt(R(0))
     334            0
     335            sage: sqrt(R(-1))
     336            4
     337           
     338            sage: sqrt(t^4)
     339            t^2
     340            sage: sqrt(4*t^4/(1+t)^8)
     341            4*t^2/(t^4 + 4*t^3 + 6*t^2 + 4*t + 1)
     342           
     343            sage: R.<t> = Fpt(5)
     344            sage: [a for a in R.iter(2) if (a^2).sqrt() not in (a,-a)]
     345            []
     346            sage: [a for a in R.iter(2) if a.is_square() and a.sqrt()^2 != a]
     347            []
     348
     349        """
     350        if zmod_poly_degree(self._numer) == -1:
     351            return self
     352        if not zmod_poly_sqrt_check(self._numer) or not zmod_poly_sqrt_check(self._denom):
     353            return None
     354        cdef zmod_poly_t numer
     355        cdef zmod_poly_t denom
     356        cdef FptElement res
     357        try:
     358            zmod_poly_init(denom, self.p)
     359            zmod_poly_init(numer, self.p)
     360            if not zmod_poly_sqrt0(numer, self._numer):
     361                return None
     362            if not zmod_poly_sqrt0(denom, self._denom):
     363                return None
     364            res = self._new_c()
     365            zmod_poly_swap(numer, res._numer)
     366            zmod_poly_swap(denom, res._denom)
     367            return res
     368        finally:
     369            zmod_poly_clear(numer)
     370            zmod_poly_clear(denom)
     371   
     372    cpdef bint is_square(self):
     373        return self._sqrt_or_None() is not None
     374   
     375    def sqrt(self):
     376        """
     377        EXAMPLES::
     378       
     379            sage: from sage.rings.fraction_field_FpT import *
     380            sage: R.<t> = Fpt(7)
     381        """
     382        s = self._sqrt_or_None()
     383        if s is None:
     384            raise ValueError, "not a perfect square"
     385        else:
     386            return s
     387   
     388#    def __pow__(FptElement self, Py_ssize_t n):
     389#       
     390
     391
     392       
     393cdef class Fpt_iter:
     394   
     395    cdef parent
     396    cdef long degree
     397    cdef FptElement cur
     398   
     399    def __init__(self, parent, degree=None, FptElement start=None):
     400        self.parent = parent
     401        self.cur = start
     402        self.degree = sys.maxint if degree is None else degree
     403   
     404    def __iter__(self):
     405        return self
     406   
     407    def __next__(self):
     408        if self.cur is None:
     409            self.cur = self.parent(0)
     410        else:
     411            self.cur = self.cur.next()
     412            if zmod_poly_degree(self.cur._numer) > self.degree:
     413                raise StopIteration
     414        return self.cur
     415
     416
     417
     418cdef inline bint normalize(zmod_poly_t numer, zmod_poly_t denom, long p):
     419    cdef long a
     420    if zmod_poly_degree(numer) == -1:
     421        zmod_poly_truncate(denom, 0)
     422        zmod_poly_set_coeff_ui(denom, 0, 1)
     423    elif zmod_poly_degree(numer) == 0 or zmod_poly_degree(denom) == 0:
     424        if zmod_poly_leading(denom) != 1:
     425            a = mod_inverse_int(zmod_poly_leading(denom), p)
     426            zmod_poly_scalar_mul(numer, numer, a)
     427            zmod_poly_scalar_mul(denom, denom, a)
     428            return True
     429        return False
     430    cdef zmod_poly_t g
     431    cdef bint changed = False
     432    try:
     433        zmod_poly_init(g, p)
     434        zmod_poly_gcd(g, numer, denom)
     435        if zmod_poly_degree(g) != 0:
     436            # Divide knowing divisible by? Can we get these quotients as a byproduct of the gcd?
     437            zmod_poly_div(numer, numer, g)
     438            zmod_poly_div(denom, denom, g)
     439            changed = True
     440        if zmod_poly_leading(denom) != 1:
     441            a = mod_inverse_int(zmod_poly_leading(denom), p)
     442            zmod_poly_scalar_mul(numer, numer, a)
     443            zmod_poly_scalar_mul(denom, denom, a)
     444            changed = True
     445        return changed
     446    finally:
     447        zmod_poly_clear(g)
     448
     449cdef inline unsigned long zmod_poly_leading(zmod_poly_t poly):
     450    return zmod_poly_get_coeff_ui(poly, zmod_poly_degree(poly))
     451
     452cdef inline void zmod_poly_inc(zmod_poly_t poly, bint monic):
     453    cdef long n
     454    cdef long a
     455    cdef long p = poly.p
     456    for n from 0 <= n <= zmod_poly_degree(poly) + 1:
     457        a = zmod_poly_get_coeff_ui(poly, n) + 1
     458        if a == poly.p:
     459            zmod_poly_set_coeff_ui(poly, n, 0)
     460        else:
     461            zmod_poly_set_coeff_ui(poly, n, a)
     462            break
     463    if monic and a == 2 and n == zmod_poly_degree(poly):
     464        zmod_poly_set_coeff_ui(poly, n, 0)
     465        zmod_poly_set_coeff_ui(poly, n+1, 1)
     466   
     467cdef void zmod_poly_derivative(zmod_poly_t diff, zmod_poly_t poly):
     468    cdef long n
     469    zmod_poly_zero(diff)
     470    for n from zmod_poly_degree(poly) >= n >= 1:
     471        zmod_poly_set_coeff_ui(diff, n-1, n * zmod_poly_get_coeff_ui(poly, n))
     472
     473cdef long sqrt_mod_int(long a, long p) except -1:
     474    return GF(p)(a).sqrt()
     475
     476cdef bint zmod_poly_sqrt_check(zmod_poly_t poly):
     477    """
     478    Quick check to see if poly could possibly be a square.
     479    """
     480    return (zmod_poly_degree(poly) % 2 == 0
     481        and jacobi_int(zmod_poly_leading(poly), poly.p) == 1
     482        and jacobi_int(zmod_poly_get_coeff_ui(poly, 0), poly.p) != -1)
     483 
     484cdef bint zmod_poly_sqrt(zmod_poly_t res, zmod_poly_t poly):
     485    """
     486    Compute the square root of poly as res if res is a perfect square.
     487   
     488    Returns True on success, False on failure.
     489    """
     490    if not zmod_poly_sqrt_check(poly):
     491        return False
     492    return zmod_poly_sqrt0(res, poly)
     493
     494cdef bint zmod_poly_sqrt0(zmod_poly_t res, zmod_poly_t poly):
     495    """
     496    Compute the square root of poly as res if res is a perfect square, assuming that poly passes zmod_poly_sqrt_check.
     497   
     498    Returns True on success, False on failure.
     499    """
     500    if zmod_poly_degree(poly) == -1:
     501        zmod_poly_zero(res)
     502        return True
     503    if zmod_poly_degree(poly) == 0:
     504        zmod_poly_zero(res)
     505        zmod_poly_set_coeff_ui(res, 0, sqrt_mod_int(zmod_poly_get_coeff_ui(poly, 0), poly.p))
     506        return True
     507    cdef zmod_poly_t g
     508    cdef long p = poly.p
     509    cdef long n, leading
     510    try:
     511        zmod_poly_init(g, p)
     512        zmod_poly_derivative(res, poly)
     513        zmod_poly_gcd(g, res, poly)
     514        if 2*zmod_poly_degree(g) < zmod_poly_degree(poly):
     515            return False
     516        elif 2*zmod_poly_degree(g) > zmod_poly_degree(poly):
     517            # This is the messy case, for now we'll just factor the thing.
     518            # print "calling factor..."
     519            sage_poly = GF(p)['t']([zmod_poly_get_coeff_ui(poly, n) for n from 0 <= n <= zmod_poly_degree(poly)])
     520            sage_g = 1
     521            for f, e in sage_poly.factor():
     522                if e % 2:
     523                    return False
     524                sage_g *= f**(e//2)
     525            zmod_poly_zero(res)
     526            for n, a in enumerate(sage_g):
     527                zmod_poly_set_coeff_ui(res, n, a)
     528            zmod_poly_scalar_mul(res, res, sqrt_mod_int(zmod_poly_leading(poly), poly.p))
     529            return True
     530        else: # deg(g) == deg(poly)/2
     531            zmod_poly_sqr(res, g)
     532            leading = zmod_poly_leading(poly) * mod_inverse_int(zmod_poly_leading(res), p)
     533            zmod_poly_scalar_mul(res, res, leading)
     534            if zmod_poly_equal(res, poly):
     535                zmod_poly_scalar_mul(res, g, sqrt_mod_int(leading, p))
     536                return True
     537            else:
     538                return False
     539    finally:
     540        zmod_poly_clear(g)
  • sage/rings/integer_mod.pxd

    diff -r 380d1bf0c2d5 -r ed7c7bc46a5b sage/rings/integer_mod.pxd
    a b  
    4646    cdef void set_from_int(IntegerMod_int64 self, int_fast64_t value)
    4747    cdef int_fast64_t get_int_value(IntegerMod_int64 self)
    4848    cdef IntegerMod_int64 _new_c(self, int_fast64_t value)
     49
     50cdef int jacobi_int(int_fast32_t a, int_fast32_t m) except -2
     51cdef int mod_inverse_int(int_fast32_t x, int_fast32_t n) except 0