Ticket #6941: trac_6941_monicgcd.patch

File trac_6941_monicgcd.patch, 9.2 KB (added by spancratz, 9 years ago)
  • sage/libs/flint/zmod_poly_linkage.pxi

    # HG changeset patch
    # User Sebastian Pancratz <sage@pancratz@org>
    # Date 1253105442 -7200
    # Node ID 8525a44c11824a40b4e5e448b1ab8c5fd6eacc00
    # Parent  70bbd9ed6103083d718d454582d667ad12ab9f39
    GCD and XGCD for polynomial rings via templating return monic versions.
    
    diff -r 70bbd9ed6103 -r 8525a44c1182 sage/libs/flint/zmod_poly_linkage.pxi
    a b  
    547547        sage: (G//d)*d == G
    548548        True
    549549    """
    550     if celement_is_zero(b, n):
    551         zmod_poly_set(res, a)
    552         return 0
    553 
     550    cdef unsigned long lead, modulus
     551   
     552    if celement_is_zero(a, n):
     553        if celement_is_zero(b, n):
     554            zmod_poly_zero(res)
     555            return 0
     556        else:
     557            lead = zmod_poly_get_coeff_ui(b, zmod_poly_degree(b))
     558            modulus = zmod_poly_modulus(b)
     559            if z_gcd(modulus, lead) == 1:
     560                zmod_poly_make_monic(res, b)
     561            else:
     562                zmod_poly_set(res, b)
     563            return 0
     564    else:
     565        if celement_is_zero(b, n):
     566            lead = zmod_poly_get_coeff_ui(a, zmod_poly_degree(a))
     567            modulus = zmod_poly_modulus(a)
     568            if z_gcd(modulus, lead) == 1:
     569                zmod_poly_make_monic(res, a)
     570            else:
     571                zmod_poly_set(res, a)
     572            return 0
     573        else:
     574            pass
     575   
    554576    zmod_poly_gcd(res, a, b)
    555     cdef unsigned long leadcoeff = zmod_poly_get_coeff_ui(res, zmod_poly_degree(res))
    556     cdef unsigned long modulus = zmod_poly_modulus(res)
    557     if z_gcd(modulus,leadcoeff) == 1:
     577    lead = zmod_poly_get_coeff_ui(res, zmod_poly_degree(res))
     578    modulus = zmod_poly_modulus(res)
     579    if z_gcd(modulus, lead) == 1:
    558580        zmod_poly_make_monic(res, res)
    559581
    560582cdef inline int celement_xgcd(zmod_poly_t res, zmod_poly_t s, zmod_poly_t t, zmod_poly_t a, zmod_poly_t b, unsigned long n) except -2:
     
    596618        sage: (G//d)*d == G
    597619        True
    598620    """
     621    cdef unsigned long lead, leadinv, modulus
     622   
     623    # Deal with the special cases
     624    #   a == 0, b == 0
     625    #   a == 0, b != 0
     626    #   a != 0, b == 0
     627    #
     628    if celement_is_zero(a, n):
     629        if celement_is_zero(b, n):
     630            zmod_poly_zero(res)
     631            zmod_poly_zero(s)
     632            zmod_poly_zero(t)
     633            return 0
     634        else:
     635            lead = zmod_poly_get_coeff_ui(b, zmod_poly_degree(b))
     636            modulus = zmod_poly_modulus(b)
     637            if z_gcd(modulus, lead) == 1:
     638                leadinv = z_invert(lead , modulus)
     639                zmod_poly_scalar_mul(res, b, modulus)
     640                zmod_poly_zero(s)
     641                zmod_poly_zero(t)
     642                zmod_poly_set_coeff_ui(t, 0, leadinv)
     643            else:
     644                zmod_poly_set(res, b)
     645                zmod_poly_zero(s)
     646                zmod_poly_zero(t)
     647                zmod_poly_set_coeff_ui(t, 0, 1)
     648            return 0
     649    else:
     650        if celement_is_zero(b, n):
     651            celement_xgcd(res, t, s, b, a, n)
     652            return 0
     653        else:
     654            pass
     655   
     656    # General case
     657    #
    599658    zmod_poly_xgcd(res, s, t, a, b)
    600659
    601 
    602660cdef inline unsigned long zmod_poly_evaluate_horner(zmod_poly_t _poly, unsigned long c):
    603661    """
    604662    Evaluate _poly at c using Horner's rule.
  • sage/rings/polynomial/polynomial_template.pxi

    diff -r 70bbd9ed6103 -r 8525a44c1182 sage/rings/polynomial/polynomial_template.pxi
    a b  
    242242
    243243    def gcd(self, Polynomial_template other):
    244244        """
    245         Return the greatest common divisor of self and other.
    246 
     245        Return the (monic) greatest common divisor of ``self`` and ``other``.
     246       
     247        NOTES:
     248       
     249        There is no logic implemented at this level;  the call is simply
     250        forwarded to the *linkage* file.
     251       
    247252        EXAMPLE::
    248253
    249254            sage: P.<x> = GF(2)[]
     
    253258            sage: f.gcd(x^2)
    254259            x
    255260        """
    256         if(celement_is_zero(&self.x, get_cparent((<Polynomial_template>self)._parent))):
    257             return other
    258         if(celement_is_zero(&other.x, get_cparent((<Polynomial_template>self)._parent))):
    259             return self
    260 
    261261        cdef Polynomial_template r = <Polynomial_template>PY_NEW(self.__class__)
    262262        celement_construct(&r.x, get_cparent((<Polynomial_template>self)._parent))
    263263        r._parent = (<Polynomial_template>self)._parent
     
    267267
    268268    def xgcd(self, Polynomial_template other):
    269269        """
    270         Computes extended gcd of self and other.
    271 
     270        Computes the extended gcd of ``self`` and ``other``.
     271       
     272        It is expected that, if possible, the greatest common divisor is
     273        returned in monic form.
     274       
     275        NOTES:
     276       
     277        There is no logic implemented at this level:  the call is simply
     278        forwarded to the *linkage* file.
     279       
    272280        EXAMPLE::
    273281
    274282            sage: P.<x> = GF(7)[]
     
    278286            sage: f.xgcd(x^2)
    279287            (x, 1, 6)
    280288        """
    281         if(celement_is_zero(&self.x, get_cparent((<Polynomial_template>self)._parent))):
    282             return other, self._parent(0), self._parent(1)
    283         if(celement_is_zero(&other.x, get_cparent((<Polynomial_template>self)._parent))):
    284             return self, self._parent(1), self._parent(0)
    285 
    286289        cdef Polynomial_template r = <Polynomial_template>PY_NEW(self.__class__)
    287290        celement_construct(&r.x, get_cparent((<Polynomial_template>self)._parent))
    288291        r._parent = (<Polynomial_template>self)._parent
     
    300303        #assert(r._parent(rp) == r)
    301304        #assert(s._parent(sp) == s)
    302305        #assert(t._parent(tp) == t)
    303         return r,s,t
     306        return r, s, t
    304307
    305308    def __floordiv__(self, right):
    306309        """
  • sage/schemes/hyperelliptic_curves/jacobian_morphism.py

    diff -r 70bbd9ed6103 -r 8525a44c1182 sage/schemes/hyperelliptic_curves/jacobian_morphism.py
    a b  
    244244    a2, b2 = D2
    245245    if a1 == a2 and b1 == b2:
    246246        # Duplication law:
    247         d, h1, h3 = a1.xgcd(2*b1)
     247        # EDIT: Mimic the old XGCD behaviour
     248        tempdouble = 2*b1
     249        if a1 == 0:
     250            d  = tempdouble
     251            h1 = 0
     252            h3 = 1
     253        elif tempdouble == 0:
     254            d  = a1
     255            h1 = 1
     256            h3 = 0
     257        else:
     258            d, h1, h3 = a1.xgcd(tempdouble)
     259       
    248260        a = (a1 // d)**2
    249261        b = (b1 + h3*((f - b1**2) // d)) % (a)
    250262    else:
    251         d0, _, h2 = a1.xgcd(a2)
     263        # EDIT: Mimic the old XGCD behaviour
     264        if a1 == 0:
     265            d0 = a2
     266            h2 = 1
     267        elif a2 == 0:
     268            d0 = a1
     269            h2 = 0
     270        else:
     271            d0, _, h2 = a1.xgcd(a2)
     272       
    252273        if d0 == 1:
    253274            a = a1*a2
    254275            b = (b2 + h2*a2*(b1-b2)) % (a)
    255276        else:
    256             d, l, h3 = d0.xgcd(b1 + b2)
     277            # EDIT: Mimic the old XGCD behaviour
     278            tempsum = b1 + b2
     279            if d0 == 0:
     280                d  = tempsum
     281                l  = 0
     282                h3 = 1
     283            elif tempsum == 0:
     284                d  = d0
     285                l  = 1
     286                h3 = 0
     287            else:
     288                d, l, h3 = d0.xgcd(tempsum)
     289           
    257290            a = (a1*a2) // (d**2)
    258291            b = ((b2 + l*h2*(b1-b2)*(a2 // d)) + h3*((f - b2**2) // d)) % (a)
    259     a =a.monic()
     292    a = a.monic()
    260293    return (a, b)
    261294
    262295def cantor_composition(D1,D2,f,h,genus):
     
    308341    a2, b2 = D2
    309342    if a1 == a2 and b1 == b2:
    310343        # Duplication law:
    311         d, h1, h3 = a1.xgcd(2*b1 + h)
     344        # EDIT: Mimic the old XGCD behaviour
     345        tempsum = 2*b1 + h
     346        if a1 == 0:
     347            d = tempsum
     348            h1 = 0
     349            h3 = 1
     350        elif tempsum == 0:
     351            d  = a1
     352            h1 = 1
     353            h3 = 0
     354        else:
     355            d, h1, h3 = a1.xgcd(tempsum)
     356       
    312357        a = (a1 // d)**2;
    313358        b = (b1 + h3*((f-h*b1-b1**2) // d)) % (a)
    314359    else:
    315         d0, _, h2 = a1.xgcd(a2)
     360        # EDIT: Mimic the old XGCD behaviour
     361        if a1 == 0:
     362            d0 = a2
     363            h2 = 1
     364        elif a2 == 0:
     365            d0 = a1
     366            h2 = 0
     367        else:
     368            d0, _, h2 = a1.xgcd(a2)
     369       
    316370        if d0 == 1:
    317371            a = a1*a2;
    318372            b = (b2 + h2*a2*(b1-b2)) % (a)
     
    322376                a = (a1*a2) // (d0**2)
    323377                b = (b2 + h2*(b1-b2)*(a2 // d0)) % (a)
    324378            else:
    325                 d, l, h3 = d0.xgcd(e0)
     379                # EDIT: Mimic the old XGCD behaviour
     380                if d0 == 0:
     381                    d  = e0
     382                    l  = 0
     383                    h3 = 1
     384                elif e0 == 0:
     385                    d  = d0
     386                    l  = 1
     387                    h3 = 0
     388                else:
     389                    d, l, h3 = d0.xgcd(e0)
     390               
    326391                a = (a1*a2) // (d**2)
    327392                b = (b2 + l*h2*(b1-b2)*(a2 // d) + h3*((f-h*b2-b2**2) // d)) % (a)
    328393    a = a.monic()