Ticket #7585: 7585_2_FpT-more.patch

File 7585_2_FpT-more.patch, 8.1 KB (added by David Roe, 13 years ago)

Rebased against 4.3.rc0

  • sage/libs/flint/zmod_poly.pxd

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1260259767 28800
    # Node ID fdab549b2310da46ca701a8e32cd5ccb70617852
    # Parent  a1af5bd038a73daff7bf00be620c9ed56615ac25
    More work on F_p(t)
    
    Fix memory leak, misc speedups.
    
    diff -r 0a8ab0fcd2e8 sage/libs/flint/zmod_poly.pxd
    a b  
    242242    cdef int zmod_poly_gcd_invert(zmod_poly_t res, zmod_poly_t poly1, zmod_poly_t poly2)
    243243    cdef void zmod_poly_xgcd(zmod_poly_t res, zmod_poly_t s, zmod_poly_t t, zmod_poly_t poly1, zmod_poly_t poly2)
    244244
     245    #
     246    # Derivative
     247    #
     248   
     249    cdef void zmod_poly_derivative(zmod_poly_t diff, zmod_poly_t poly)
     250
    245251    # Composition / evaluation
    246252
    247253    cdef unsigned long zmod_poly_evaluate(zmod_poly_t, unsigned long)
     
    265271
    266272    # FLINT 1.1 will have:
    267273    # cdef void zmod_poly_powmod(zmod_poly_t res, zmod_poly_t pol, long exp, zmod_poly_t f)
     274
     275    #
     276    # Factorisation/Irreducibility
     277    #
     278   
     279    ctypedef struct zmod_poly_factor_struct:
     280        zmod_poly_t* factors
     281        unsigned long * exponents
     282        unsigned long alloc
     283        unsigned long num_factors
     284
     285    ctypedef zmod_poly_factor_struct zmod_poly_factor_t[1]
     286
     287    void zmod_poly_factor_init(zmod_poly_factor_t fac)
     288    void zmod_poly_factor_clear(zmod_poly_factor_t fac)
     289
     290    bint zmod_poly_isirreducible(zmod_poly_t f)
     291    void zmod_poly_factor_add(zmod_poly_factor_t fac, zmod_poly_t poly)
     292    void zmod_poly_factor_concat(zmod_poly_factor_t res, zmod_poly_factor_t fac)
     293    void zmod_poly_factor_print(zmod_poly_factor_t fac)
     294    void zmod_poly_factor_pow(zmod_poly_factor_t fac, unsigned long exp)
     295   
     296    void zmod_poly_factor_square_free(zmod_poly_factor_t res, zmod_poly_t f)
     297    void zmod_poly_factor_berlekamp(zmod_poly_factor_t factors, zmod_poly_t f)
     298    unsigned long zmod_poly_factor(zmod_poly_factor_t result, zmod_poly_t input)
  • sage/rings/fraction_field_FpT.pyx

    diff -r 0a8ab0fcd2e8 sage/rings/fraction_field_FpT.pyx
    a b  
    8585        cdef FptElement x = <FptElement>PY_NEW(FptElement)
    8686        x._parent = self._parent
    8787        x.p = self.p
    88         zmod_poly_init(x._numer, x.p)
    89         zmod_poly_init(x._denom, x.p)
    90         self.initalized = True
     88        zmod_poly_init_precomp(x._numer, x.p, self._numer.p_inv)
     89        zmod_poly_init_precomp(x._denom, x.p, self._numer.p_inv)
     90        x.initalized = True
     91        return x
     92   
     93    cdef FptElement _copy_c(self):
     94        cdef FptElement x = <FptElement>PY_NEW(FptElement)
     95        x._parent = self._parent
     96        x.p = self.p
     97        zmod_poly_init2_precomp(x._numer, x.p, self._numer.p_inv, self._numer.length)
     98        zmod_poly_init2_precomp(x._denom, x.p, self._denom.p_inv, self._denom.length)
     99        zmod_poly_set(x._numer, self._numer)
     100        zmod_poly_set(x._denom, self._denom)
     101        x.initalized = True
    91102        return x
    92103   
    93104    def numer(self):
     
    454465    cdef parent
    455466    cdef long degree
    456467    cdef FptElement cur
     468    cdef zmod_poly_t g
    457469   
    458470    def __init__(self, parent, degree=None, FptElement start=None):
    459471        if degree is None:
     
    462474        self.cur = start
    463475        self.degree = sys.maxint if degree is None else degree
    464476   
     477    def __cinit__(self, parent, *args, **kwds):
     478        zmod_poly_init(self.g, parent.characteristic())
     479   
     480    def __dealloc__(self):
     481        zmod_poly_clear(self.g)
     482   
    465483    def __iter__(self):
    466484        return self
    467485   
     
    513531            sage: L[-1]
    514532            (4*t^3 + 4*t^2 + 4*t + 4)/(t^3 + 4*t^2 + 4*t + 4)
    515533        """
    516         cdef zmod_poly_t g
    517534        cdef FptElement next
    518535        if self.cur is None:
    519536            self.cur = self.parent(0)
    520537        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)
     538            next = self.cur._copy_c()
    525539            while True:
    526540                zmod_poly_inc(next._numer, False)
    527541                if zmod_poly_degree(next._numer) > self.degree:
    528542                    zmod_poly_inc(next._denom, True)
    529543                    if zmod_poly_degree(next._denom) > self.degree:
    530544                        raise StopIteration
    531                         zmod_poly_clear(g)
    532545                    zmod_poly_zero(next._numer)
    533546                    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):
     547                zmod_poly_gcd(self.g, next._numer, next._denom)
     548                if zmod_poly_is_one(self.g):
    536549                    break
    537             zmod_poly_clear(g)
    538550            self.cur = next
    539551#            self.cur = self.cur.next()
    540552#            if zmod_poly_degree(self.cur._numer) > self.degree:
     
    558570    cdef zmod_poly_t g
    559571    cdef bint changed = False
    560572    try:
    561         zmod_poly_init(g, p)
     573        zmod_poly_init_precomp(g, p, numer.p_inv)
    562574        zmod_poly_gcd(g, numer, denom)
    563575        if zmod_poly_degree(g) != 0:
    564576            # Divide knowing divisible by? Can we get these quotients as a byproduct of the gcd?
     
    591603    if monic and a == 2 and n == zmod_poly_degree(poly):
    592604        zmod_poly_set_coeff_ui(poly, n, 0)
    593605        zmod_poly_set_coeff_ui(poly, n+1, 1)
    594    
    595 cdef 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))
    600606
    601607cdef void zmod_poly_pow(zmod_poly_t res, zmod_poly_t poly, unsigned long e):
    602608    if zmod_poly_degree(poly) < 2:
     
    680686    cdef zmod_poly_t g
    681687    cdef long p = poly.p
    682688    cdef long n, leading
     689    cdef zmod_poly_factor_t factors
    683690    try:
    684691        zmod_poly_init(g, p)
    685692        zmod_poly_derivative(res, poly)
    686693        zmod_poly_gcd(g, res, poly)
    687694        if 2*zmod_poly_degree(g) < zmod_poly_degree(poly):
    688695            return False
     696           
    689697        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
     698            try:
     699                zmod_poly_factor_init(factors)
     700                leading = zmod_poly_leading(poly)
     701                if leading == 1:
     702                    zmod_poly_factor_square_free(factors, poly)
     703                else:
     704                    zmod_poly_scalar_mul(res, poly, mod_inverse_int(leading, p))
     705                    zmod_poly_factor_square_free(factors, res)
     706                for n in range(factors.num_factors):
     707                    if factors.exponents[n] % 2 != 0:
     708                        return False
     709                zmod_poly_pow(res, factors.factors[0], factors.exponents[0]//2)
     710                for n in range(1, factors.num_factors):
     711                    zmod_poly_pow(factors.factors[n], factors.factors[n], factors.exponents[n]//2)
     712                    zmod_poly_mul(res, res, factors.factors[n])
     713                if leading != 1:
     714                    zmod_poly_scalar_mul(res, res, sqrt_mod_int(zmod_poly_leading(poly), p))
     715                return True
     716            finally:
     717                zmod_poly_factor_clear(factors)
     718           
    703719        else: # deg(g) == deg(poly)/2
    704720            zmod_poly_sqr(res, g)
    705721            leading = zmod_poly_leading(poly) * mod_inverse_int(zmod_poly_leading(res), p)