Ticket #7585: 7585_1_FpT-orig.patch

File 7585_1_FpT-orig.patch, 25.0 KB (added by David Roe, 13 years ago)

Rebased against 4.3.rc0

  • module_list.py

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1260035621 28800
    # Node ID a1af5bd038a73daff7bf00be620c9ed56615ac25
    # Parent  380d1bf0c2d5e07b6bab0e3d0df3f6512248af8d
    Fast arithmatic in F_p(t).
    
    diff -r 380d1bf0c2d5 -r a1af5bd038a7 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 a1af5bd038a7 sage/rings/fraction_field_FpT.pyx
    - +  
     1import sys
     2
     3from sage.rings.all import GF
     4from sage.rings.ring cimport Field
     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, mod_pow_int
     10
     11class Fpt(Field):
     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 characteristic(self):
     33        return self.p
     34   
     35    def _coerce_map_from_(self, R):
     36        return self.poly_ring.has_coerce_map_from(R)
     37   
     38    def __iter__(self):
     39        return self.iter()
     40   
     41    def iter(self, bound=None, start=None):
     42        """
     43            sage: from sage.rings.fraction_field_FpT import *
     44            sage: R.<t> = Fpt(5)
     45            sage: list(R.iter(2))
     46        """
     47        return Fpt_iter(self, bound, start)
     48
     49cdef class FptElement(RingElement):
     50
     51    cdef long p
     52    cdef zmod_poly_t _numer, _denom
     53    cdef bint initalized
     54
     55    def __init__(self, parent, numer, denom=1):
     56        """
     57        EXAMPLES::
     58       
     59            sage: from sage.rings.fraction_field_FpT import *
     60            sage: R.<t> = Fpt(5)
     61            sage: R(7)
     62            2
     63           
     64        """
     65        RingElement.__init__(self, parent)
     66        numer = parent.poly_ring(numer)
     67        denom = parent.poly_ring(denom)
     68        self.p = parent.p
     69        zmod_poly_init(self._numer, self.p)
     70        zmod_poly_init(self._denom, self.p)
     71        self.initalized = True
     72        cdef long n
     73        for n, a in enumerate(numer):
     74            zmod_poly_set_coeff_ui(self._numer, n, a)
     75        for n, a in enumerate(denom):
     76            zmod_poly_set_coeff_ui(self._denom, n, a)
     77        normalize(self._numer, self._denom, self.p)
     78
     79    def __dealloc__(self):
     80        if self.initalized:
     81            zmod_poly_clear(self._numer)
     82            zmod_poly_clear(self._denom)
     83   
     84    cdef FptElement _new_c(self):
     85        cdef FptElement x = <FptElement>PY_NEW(FptElement)
     86        x._parent = self._parent
     87        x.p = self.p
     88        zmod_poly_init(x._numer, x.p)
     89        zmod_poly_init(x._denom, x.p)
     90        self.initalized = True
     91        return x
     92   
     93    def numer(self):
     94        cdef long n
     95        return self._parent.poly_ring(
     96            [zmod_poly_get_coeff_ui (self._numer, n) for n in range(0, zmod_poly_degree(self._numer)+1)])
     97
     98    def denom(self):
     99        cdef long n
     100        return self._parent.poly_ring(
     101            [zmod_poly_get_coeff_ui (self._denom, n) for n in range(0, zmod_poly_degree(self._denom)+1)])
     102
     103    def _repr_(self):
     104        """
     105            sage: from sage.rings.fraction_field_FpT import *
     106            sage: R.<t> = Fpt(17)
     107            sage: -t
     108            16*t
     109            sage: 1/t
     110            1/t
     111            sage: 1/(t+1)
     112            1/(t + 1)
     113            sage: 1-t/t
     114            0
     115            sage: (1-t)/t
     116            (16*t + 1)/t
     117        """
     118        if zmod_poly_degree(self._denom) == 0 and zmod_poly_get_coeff_ui(self._denom, 0) == 1:
     119            return repr(self.numer())
     120        else:
     121            numer_s = repr(self.numer())
     122            denom_s = repr(self.denom())
     123            if '-' in numer_s or '+' in numer_s:
     124                numer_s = "(%s)" % numer_s
     125            if '-' in denom_s or '+' in denom_s:
     126                denom_s = "(%s)" % denom_s
     127            return "%s/%s" % (numer_s, denom_s)
     128   
     129    def _latex_(self):
     130        if zmod_poly_degree(self._denom) == 0 and zmod_poly_get_coeff_ui(self._denom, 0) == 1:
     131            return self.numer()._latex_()
     132        else:
     133            return "{%s}/{%s}" % (self.numer()._latex_(), self.denom()._latex_())
     134        return
     135   
     136    def __cmp__(self, other):
     137        """
     138        TESTS::
     139       
     140            sage: from sage.rings.fraction_field_FpT import *
     141            sage: R.<t> = Fpt(7)
     142            sage: t == t
     143            True
     144            sage: t == -t
     145            False
     146            sage: -t == 6*t
     147            True
     148            sage: 1/t == 1/t
     149            True
     150            sage: 1/t == 1/(t+1)
     151            False
     152        """
     153        if isinstance(other, FptElement):
     154            # They are normalized.
     155            if (zmod_poly_equal(self._numer, (<FptElement>other)._numer) and
     156                    zmod_poly_equal(self._denom, (<FptElement>other)._denom)):
     157                return 0
     158            else:
     159                return -1
     160        else:
     161            return cmp(type(self), type(other))
     162   
     163    def __hash__(self):
     164        """
     165        TESTS::
     166           
     167            sage: from sage.rings.fraction_field_FpT import *
     168            sage: K.<t> = Fpt(7)
     169            sage: hash(K(0))
     170            0
     171            sage: hash(K(5))
     172            5
     173            sage: set([1, t, 1/t, t, t, 1/t, 1+1/t, t/t])
     174            set([1, (t + 1)/t, t, 1/t])
     175        """
     176        if self.denom() == 1:
     177            return hash(self.numer())
     178        return hash(str(self))
     179       
     180    def __neg__(self):
     181        cdef FptElement x = self._new_c()
     182        zmod_poly_neg(x._numer, self._numer)
     183        zmod_poly_set(x._denom, self._denom)
     184        return x
     185       
     186    def __invert__(self):
     187        if zmod_poly_degree(self._numer) == -1:
     188            raise ZeroDivisionError
     189        cdef FptElement x = self._new_c()
     190        zmod_poly_set(x._denom, self._numer)
     191        zmod_poly_set(x._numer, self._denom)
     192        return x
     193
     194    cpdef ModuleElement _add_(self, ModuleElement _other):
     195        """
     196        EXAMPLES::
     197       
     198            sage: from sage.rings.fraction_field_FpT import *
     199            sage: R.<t> = Fpt(7)
     200            sage: t + t
     201            2*t
     202            sage: (t + 3) + (t + 10)
     203            2*t + 6
     204            sage: sum([t] * 7)
     205            0
     206            sage: 1/t + t
     207            (t^2 + 1)/(t)
     208            sage: 1/t + 1/t^2
     209            (t + 1)/(t^2)
     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_add(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 ModuleElement _sub_(self, ModuleElement _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            0
     228            sage: (t + 3) - (t + 11)
     229            sage: 6
     230        """
     231        cdef FptElement other = <FptElement>_other
     232        cdef FptElement x = self._new_c()
     233        zmod_poly_mul(x._numer, self._numer, other._denom)
     234        zmod_poly_mul(x._denom, self._denom, other._numer) # use x._denom as a temp
     235        zmod_poly_sub(x._numer, x._numer, x._denom)
     236        zmod_poly_mul(x._denom, self._denom, other._denom)
     237        normalize(x._numer, x._denom, self.p)
     238        return x
     239
     240    cpdef RingElement _mul_(self, RingElement _other):
     241        """
     242        EXAMPLES::
     243       
     244            sage: from sage.rings.fraction_field_FpT import *
     245            sage: R.<t> = Fpt(7)
     246            sage: t * t
     247            t^2
     248            sage: (t + 3) * (t + 10)
     249            t^2 + 6*t + 2
     250        """
     251        cdef FptElement other = <FptElement>_other
     252        cdef FptElement x = self._new_c()
     253        zmod_poly_mul(x._numer, self._numer, other._numer)
     254        zmod_poly_mul(x._denom, self._denom, other._denom)
     255        normalize(x._numer, x._denom, self.p)
     256        return x
     257
     258    cpdef RingElement _div_(self, RingElement _other):
     259        """
     260        EXAMPLES::
     261       
     262            sage: from sage.rings.fraction_field_FpT import *
     263            sage: R.<t> = Fpt(5)
     264            sage: t / t
     265            1
     266            sage: (t + 3) / (t + 11)
     267            (t + 3)/(t + 1)
     268            sage: (t^2 + 2*t + 1) / (t + 1)
     269            t + 1
     270        """
     271        cdef FptElement other = <FptElement>_other
     272        if zmod_poly_degree(other._numer) == -1:
     273            raise ZeroDivisionError
     274        cdef FptElement x = self._new_c()
     275        zmod_poly_mul(x._numer, self._numer, other._denom)
     276        zmod_poly_mul(x._denom, self._denom, other._numer)
     277        normalize(x._numer, x._denom, self.p)
     278        return x
     279   
     280    cpdef FptElement next(self):
     281        # TODO: This misses denom > numer...
     282        """
     283        EXAMPLES::
     284       
     285            sage: from sage.rings.fraction_field_FpT import *
     286            sage: R.<t> = Fpt(3)
     287            sage: a = R(0)
     288            sage: for _ in range(15):
     289            ...       a = a.next()
     290            ...       print a
     291            ...   
     292            1
     293            2
     294            t
     295            t + 1
     296            (t + 1)/(t)
     297            t + 2
     298            (t + 2)/(t)
     299            (t + 2)/(t + 1)
     300            2*t
     301            (2*t)/(t + 1)
     302            (2*t)/(t + 2)
     303            2*t + 1
     304            (2*t + 1)/(t)
     305            (2*t + 1)/(t + 1)
     306            2*t + 2
     307        """
     308        cdef FptElement x = self._new_c()
     309        if zmod_poly_degree(self._numer) == -1:
     310            zmod_poly_set_coeff_ui(x._numer, 0, 1)
     311            zmod_poly_set_coeff_ui(x._denom, 0, 1)
     312            return x
     313        elif zmod_poly_degree(self._numer) == 0:
     314            zmod_poly_set_coeff_ui(x._numer, 0, zmod_poly_get_coeff_ui(self._numer, 0) + 1)
     315            zmod_poly_set_coeff_ui(x._denom, 0, 1)
     316            if zmod_poly_get_coeff_ui(x._numer, 0) == 0:
     317                zmod_poly_set_coeff_ui(x._numer, 1, 1)
     318            return x
     319        cdef zmod_poly_t g
     320        zmod_poly_init(g, self.p)
     321        zmod_poly_set(x._numer, self._numer)
     322        zmod_poly_set(x._denom, self._denom)
     323        while True:
     324            zmod_poly_inc(x._denom, True)
     325            if zmod_poly_degree(x._numer) < zmod_poly_degree(x._denom):
     326                zmod_poly_inc(x._numer, False)
     327                zmod_poly_zero(x._denom)
     328                zmod_poly_set_coeff_ui(x._denom, 0, 1)
     329            zmod_poly_gcd(g, x._numer, x._denom)
     330            if zmod_poly_is_one(g):
     331                zmod_poly_clear(g)
     332                return x
     333   
     334    cpdef _sqrt_or_None(self):
     335        """
     336        Returns the squre root of self, or None. Differs from sqrt() by not raising an exception.
     337       
     338        TESTS::
     339       
     340            sage: from sage.rings.fraction_field_FpT import *
     341            sage: R.<t> = Fpt(17)
     342            sage: sqrt(t^2)
     343            t
     344            sage: sqrt(1/t^2)
     345            1/t
     346            sage: sqrt((1+t)^2)
     347            t + 1
     348            sage: sqrt((1+t)^2 / t^2)
     349            (t + 1)/t
     350
     351            sage: sqrt(4 * (1+t)^2 / t^2)
     352            (2*t + 2)/t
     353           
     354            sage: sqrt(R(0))
     355            0
     356            sage: sqrt(R(-1))
     357            4
     358           
     359            sage: sqrt(t^4)
     360            t^2
     361            sage: sqrt(4*t^4/(1+t)^8)
     362            4*t^2/(t^4 + 4*t^3 + 6*t^2 + 4*t + 1)
     363           
     364            sage: R.<t> = Fpt(5)
     365            sage: [a for a in R.iter(2) if (a^2).sqrt() not in (a,-a)]
     366            []
     367            sage: [a for a in R.iter(2) if a.is_square() and a.sqrt()^2 != a]
     368            []
     369
     370        """
     371        if zmod_poly_degree(self._numer) == -1:
     372            return self
     373        if not zmod_poly_sqrt_check(self._numer) or not zmod_poly_sqrt_check(self._denom):
     374            return None
     375        cdef zmod_poly_t numer
     376        cdef zmod_poly_t denom
     377        cdef FptElement res
     378        try:
     379            zmod_poly_init(denom, self.p)
     380            zmod_poly_init(numer, self.p)
     381            if not zmod_poly_sqrt0(numer, self._numer):
     382                return None
     383            if not zmod_poly_sqrt0(denom, self._denom):
     384                return None
     385            res = self._new_c()
     386            zmod_poly_swap(numer, res._numer)
     387            zmod_poly_swap(denom, res._denom)
     388            return res
     389        finally:
     390            zmod_poly_clear(numer)
     391            zmod_poly_clear(denom)
     392   
     393    cpdef bint is_square(self):
     394        return self._sqrt_or_None() is not None
     395   
     396    def sqrt(self):
     397        """
     398        EXAMPLES::
     399       
     400            sage: from sage.rings.fraction_field_FpT import *
     401            sage: R.<t> = Fpt(7)
     402        """
     403        s = self._sqrt_or_None()
     404        if s is None:
     405            raise ValueError, "not a perfect square"
     406        else:
     407            return s
     408   
     409    def __pow__(FptElement self, Py_ssize_t e, dummy):
     410        """
     411        EXAMPLES::
     412       
     413            sage: from sage.rings.fraction_field_FpT import *
     414            sage: R.<t> = Fpt(7)
     415            sage: t^5
     416            t^5
     417            sage: t^-5
     418            1/t^5
     419           
     420            sage: a = (t+1)/(t-1); a
     421            (t + 1)/(t + 6)
     422            sage: a^5
     423            (t^5 + 5*t^4 + 3*t^3 + 3*t^2 + 5*t + 1)/(t^5 + 2*t^4 + 3*t^3 + 4*t^2 + 5*t + 6)
     424            sage: a^7
     425            (t^7 + 1)/(t^7 + 6)
     426            sage: a^14
     427            (t^14 + 2*t^7 + 1)/(t^14 + 5*t^7 + 1)
     428           
     429            sage: (a^2)^2 == a^4
     430            True
     431            sage: a^3 * a^2 == a^5
     432            True
     433            sage: a^47 * a^92 == a^(47+92)
     434            True
     435        """
     436        cdef long a
     437        assert dummy is None
     438        cdef FptElement x = self._new_c()
     439        if e >= 0:
     440            zmod_poly_pow(x._numer, self._numer, e)
     441            zmod_poly_pow(x._denom, self._denom, e)
     442        else:
     443            zmod_poly_pow(x._denom, self._numer, -e)
     444            zmod_poly_pow(x._numer, self._denom, -e)
     445            if zmod_poly_leading(x._denom) != 1:
     446                a = mod_inverse_int(zmod_poly_leading(x._denom), self.p)
     447                zmod_poly_scalar_mul(x._numer, x._numer, a)
     448                zmod_poly_scalar_mul(x._denom, x._denom, a)
     449        return x
     450
     451       
     452cdef class Fpt_iter:
     453   
     454    cdef parent
     455    cdef long degree
     456    cdef FptElement cur
     457   
     458    def __init__(self, parent, degree=None, FptElement start=None):
     459        if degree is None:
     460            raise NotImplementedError
     461        self.parent = parent
     462        self.cur = start
     463        self.degree = sys.maxint if degree is None else degree
     464   
     465    def __iter__(self):
     466        return self
     467   
     468    def __next__(self):
     469        """
     470        EXAMPLES::
     471       
     472            sage: from sage.rings.fraction_field_FpT import *
     473            sage: K.<t> = Fpt(3)
     474            sage: list(K.iter(1))
     475            [0,
     476             1,
     477             2,
     478             t,
     479             t + 1,
     480             t + 2,
     481             2*t,
     482             2*t + 1,
     483             2*t + 2,
     484             1/t,
     485             2/t,
     486             (t + 1)/t,
     487             (t + 2)/t,
     488             (2*t + 1)/t,
     489             (2*t + 2)/t,
     490             1/(t + 1),
     491             2/(t + 1),
     492             t/(t + 1),
     493             (t + 2)/(t + 1),
     494             2*t/(t + 1),
     495             (2*t + 1)/(t + 1),
     496             1/(t + 2),
     497             2/(t + 2),
     498             t/(t + 2),
     499             (t + 1)/(t + 2),
     500             2*t/(t + 2),
     501             (2*t + 2)/(t + 2)]
     502             
     503            sage: len(list(K.iter(3)))
     504            2187
     505
     506            sage: K.<t> = Fpt(5)
     507            sage: L = list(K.iter(3))
     508            78125
     509            sage: L[:10]
     510            [0, 1, 2, 3, 4, t, t + 1, t + 2, t + 3, t + 4]
     511            sage: L[2000]
     512            (3*t^3 + 3*t^2 + 3*t + 4)/(t + 2)
     513            sage: L[-1]
     514            (4*t^3 + 4*t^2 + 4*t + 4)/(t^3 + 4*t^2 + 4*t + 4)
     515        """
     516        cdef zmod_poly_t g
     517        cdef FptElement next
     518        if self.cur is None:
     519            self.cur = self.parent(0)
     520        else:
     521            next = self.cur._new_c()
     522            zmod_poly_set(next._numer, self.cur._numer)
     523            zmod_poly_set(next._denom, self.cur._denom)
     524            zmod_poly_init(g, next._numer.p)
     525            while True:
     526                zmod_poly_inc(next._numer, False)
     527                if zmod_poly_degree(next._numer) > self.degree:
     528                    zmod_poly_inc(next._denom, True)
     529                    if zmod_poly_degree(next._denom) > self.degree:
     530                        raise StopIteration
     531                        zmod_poly_clear(g)
     532                    zmod_poly_zero(next._numer)
     533                    zmod_poly_set_coeff_ui(next._numer, 0, 1)
     534                zmod_poly_gcd(g, next._numer, next._denom)
     535                if zmod_poly_is_one(g):
     536                    break
     537            zmod_poly_clear(g)
     538            self.cur = next
     539#            self.cur = self.cur.next()
     540#            if zmod_poly_degree(self.cur._numer) > self.degree:
     541#                raise StopIteration
     542        return self.cur
     543
     544
     545
     546cdef inline bint normalize(zmod_poly_t numer, zmod_poly_t denom, long p):
     547    cdef long a
     548    if zmod_poly_degree(numer) == -1:
     549        zmod_poly_truncate(denom, 0)
     550        zmod_poly_set_coeff_ui(denom, 0, 1)
     551    elif zmod_poly_degree(numer) == 0 or zmod_poly_degree(denom) == 0:
     552        if zmod_poly_leading(denom) != 1:
     553            a = mod_inverse_int(zmod_poly_leading(denom), p)
     554            zmod_poly_scalar_mul(numer, numer, a)
     555            zmod_poly_scalar_mul(denom, denom, a)
     556            return True
     557        return False
     558    cdef zmod_poly_t g
     559    cdef bint changed = False
     560    try:
     561        zmod_poly_init(g, p)
     562        zmod_poly_gcd(g, numer, denom)
     563        if zmod_poly_degree(g) != 0:
     564            # Divide knowing divisible by? Can we get these quotients as a byproduct of the gcd?
     565            zmod_poly_div(numer, numer, g)
     566            zmod_poly_div(denom, denom, g)
     567            changed = True
     568        if zmod_poly_leading(denom) != 1:
     569            a = mod_inverse_int(zmod_poly_leading(denom), p)
     570            zmod_poly_scalar_mul(numer, numer, a)
     571            zmod_poly_scalar_mul(denom, denom, a)
     572            changed = True
     573        return changed
     574    finally:
     575        zmod_poly_clear(g)
     576
     577cdef inline unsigned long zmod_poly_leading(zmod_poly_t poly):
     578    return zmod_poly_get_coeff_ui(poly, zmod_poly_degree(poly))
     579
     580cdef inline void zmod_poly_inc(zmod_poly_t poly, bint monic):
     581    cdef long n
     582    cdef long a
     583    cdef long p = poly.p
     584    for n from 0 <= n <= zmod_poly_degree(poly) + 1:
     585        a = zmod_poly_get_coeff_ui(poly, n) + 1
     586        if a == poly.p:
     587            zmod_poly_set_coeff_ui(poly, n, 0)
     588        else:
     589            zmod_poly_set_coeff_ui(poly, n, a)
     590            break
     591    if monic and a == 2 and n == zmod_poly_degree(poly):
     592        zmod_poly_set_coeff_ui(poly, n, 0)
     593        zmod_poly_set_coeff_ui(poly, n+1, 1)
     594   
     595cdef void zmod_poly_derivative(zmod_poly_t diff, zmod_poly_t poly):
     596    cdef long n
     597    zmod_poly_zero(diff)
     598    for n from zmod_poly_degree(poly) >= n >= 1:
     599        zmod_poly_set_coeff_ui(diff, n-1, n * zmod_poly_get_coeff_ui(poly, n))
     600
     601cdef void zmod_poly_pow(zmod_poly_t res, zmod_poly_t poly, unsigned long e):
     602    if zmod_poly_degree(poly) < 2:
     603        if zmod_poly_degree(poly) == -1:
     604            zmod_poly_zero(res)
     605            return
     606        elif zmod_poly_is_one(poly):
     607            zmod_poly_set(res, poly)
     608            return
     609        elif e > 1 and zmod_poly_degree(poly) == 1 and zmod_poly_get_coeff_ui(poly, 0) == 0 and zmod_poly_get_coeff_ui(poly, 1) == 1:
     610            zmod_poly_left_shift(res, poly, e-1)
     611            return
     612   
     613    # TODO: Could use the fact that (a+b)^p = a^p + b^p
     614    # Only seems to be a big gain for large p, large exponents...
     615    cdef zmod_poly_t pow2
     616   
     617    if e < 5:
     618        if e == 0:
     619            zmod_poly_zero(res)
     620            zmod_poly_set_coeff_ui(res, 0, 1)
     621        elif e == 1:
     622            zmod_poly_set(res, poly)
     623        elif e == 2:
     624            zmod_poly_sqr(res, poly)
     625        elif e == 3:
     626            zmod_poly_sqr(res, poly)
     627            zmod_poly_mul(res, res, poly)
     628        elif e == 4:
     629            zmod_poly_sqr(res, poly)
     630            zmod_poly_sqr(res, res)
     631    else:
     632        zmod_poly_init(pow2, poly.p)
     633        zmod_poly_zero(res)
     634        zmod_poly_set_coeff_ui(res, 0, 1)
     635        zmod_poly_set(pow2, poly)
     636        while True:
     637            if e & 1:
     638                zmod_poly_mul(res, res, pow2)
     639            e >>= 1
     640            if e == 0:
     641                break
     642            zmod_poly_sqr(pow2, pow2)
     643        zmod_poly_clear(pow2)
     644
     645
     646cdef long sqrt_mod_int(long a, long p) except -1:
     647    return GF(p)(a).sqrt()
     648
     649cdef bint zmod_poly_sqrt_check(zmod_poly_t poly):
     650    """
     651    Quick check to see if poly could possibly be a square.
     652    """
     653    return (zmod_poly_degree(poly) % 2 == 0
     654        and jacobi_int(zmod_poly_leading(poly), poly.p) == 1
     655        and jacobi_int(zmod_poly_get_coeff_ui(poly, 0), poly.p) != -1)
     656 
     657cdef bint zmod_poly_sqrt(zmod_poly_t res, zmod_poly_t poly):
     658    """
     659    Compute the square root of poly as res if res is a perfect square.
     660   
     661    Returns True on success, False on failure.
     662    """
     663    if not zmod_poly_sqrt_check(poly):
     664        return False
     665    return zmod_poly_sqrt0(res, poly)
     666
     667cdef bint zmod_poly_sqrt0(zmod_poly_t res, zmod_poly_t poly):
     668    """
     669    Compute the square root of poly as res if res is a perfect square, assuming that poly passes zmod_poly_sqrt_check.
     670   
     671    Returns True on success, False on failure.
     672    """
     673    if zmod_poly_degree(poly) == -1:
     674        zmod_poly_zero(res)
     675        return True
     676    if zmod_poly_degree(poly) == 0:
     677        zmod_poly_zero(res)
     678        zmod_poly_set_coeff_ui(res, 0, sqrt_mod_int(zmod_poly_get_coeff_ui(poly, 0), poly.p))
     679        return True
     680    cdef zmod_poly_t g
     681    cdef long p = poly.p
     682    cdef long n, leading
     683    try:
     684        zmod_poly_init(g, p)
     685        zmod_poly_derivative(res, poly)
     686        zmod_poly_gcd(g, res, poly)
     687        if 2*zmod_poly_degree(g) < zmod_poly_degree(poly):
     688            return False
     689        elif 2*zmod_poly_degree(g) > zmod_poly_degree(poly):
     690            # This is the messy case, for now we'll just factor the thing.
     691            # print "calling factor..."
     692            sage_poly = GF(p)['t']([zmod_poly_get_coeff_ui(poly, n) for n from 0 <= n <= zmod_poly_degree(poly)])
     693            sage_g = 1
     694            for f, e in sage_poly.factor():
     695                if e % 2:
     696                    return False
     697                sage_g *= f**(e//2)
     698            zmod_poly_zero(res)
     699            for n, a in enumerate(sage_g):
     700                zmod_poly_set_coeff_ui(res, n, a)
     701            zmod_poly_scalar_mul(res, res, sqrt_mod_int(zmod_poly_leading(poly), poly.p))
     702            return True
     703        else: # deg(g) == deg(poly)/2
     704            zmod_poly_sqr(res, g)
     705            leading = zmod_poly_leading(poly) * mod_inverse_int(zmod_poly_leading(res), p)
     706            zmod_poly_scalar_mul(res, res, leading)
     707            if zmod_poly_equal(res, poly):
     708                zmod_poly_scalar_mul(res, g, sqrt_mod_int(leading, p))
     709                return True
     710            else:
     711                return False
     712    finally:
     713        zmod_poly_clear(g)
  • sage/rings/integer_mod.pxd

    diff -r 380d1bf0c2d5 -r a1af5bd038a7 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_fast32_t mod_inverse_int(int_fast32_t x, int_fast32_t n) except 0
     52cdef int_fast32_t mod_pow_int(int_fast32_t base, int_fast32_t exp, int_fast32_t n)