Ticket #6768: 12538.patch

File 12538.patch, 16.7 KB (added by wakep, 13 years ago)
  • new file sage/rings/polynomial/char_two_quad.py

    # HG changeset patch
    # User Preston Wake <preston.wake@gmail.com>
    # Date 1248373201 14400
    # Node ID 28635b3d0b70e166aab75e0492f397f051df35ce
    # Parent  2e793d2a0e123293b73eed40715e43185fd9ccfe
    added speciallized code to factor quadratic polynomials over GF(2^m)
    
    diff -r 2e793d2a0e12 -r 28635b3d0b70 sage/rings/polynomial/char_two_quad.py
    - +  
     1def quad(f):
     2    r"""
     3    Returns the factorization of f, for f a quadratic polynomial over GF(2^m).
     4    Taken from CHIN-LONG CHEN's 1982 paper Formulas for the Solutions of Quadratic Equations over GF (2^m).
     5    By Preston Wake and Wouter Castryck, 2009
     6
     7    EXAMPLES:
     8
     9    ::
     10        sage: from sage.rings.polynomial.char_two_quad import quad
     11        sage: F.<x> = GF(4)
     12        sage: PP = PolynomialRing(F,'y')
     13        sage: f = PP([1,1,1]);f
     14        y^2 + y + 1
     15        sage: a = quad(f); a
     16        ([(y + x, 1), (y + x + 1, 1)], 1)
     17        sage: r=a[0][0][0];r
     18        y + x
     19        sage: f(r(0))
     20        0
     21        sage: f/r
     22        y + x + 1
     23   
     24    One where f is a square
     25
     26    ::
     27        sage: f = PP([1,0,1]);f
     28        y^2 + 1
     29        sage: a = quad(f); a
     30        ([(y + 1, 2)], 1)
     31   
     32    One where f is irreducible
     33   
     34    ::
     35        sage: f = PP([x,1,1]);f
     36        y^2 + y + x
     37        sage: a = quad(f);a
     38        ([(y^2 + y + x, 1)], 1)
     39           
     40    """
     41    PP=f.parent()
     42    F=PP.base_ring()
     43    m=F.degree();
     44    a=f[2];b=f[1];c=f[0];
     45    if a==0:
     46        raise ValueError, 'not a quad'
     47    if b==0:
     48        y=gftwosqrt(c/a)
     49        f1=PP([y,1])
     50        return ([(f1,2)],a)
     51    k=a*c/b**2
     52    if not(k.trace()==0):
     53        return ([(f*(a**(-1)),1)],a)
     54    else:
     55        y=doquad(k,m,F)
     56    y1=y+1
     57    beta=b/a
     58    f1=PP([y*beta,1]);f2=PP([y1*beta,1])
     59    return ([(f1,1),(f2,1)],a)
     60
     61
     62def gftwosqrt(r):
     63    """
     64    Quickly finds the squareroot of an element in GF(2^m)
     65    """
     66    F=r.parent()
     67    c=F.cardinality()
     68    return r**(c/2)
     69
     70def qodd(k,m):
     71    """
     72    Does the compution for m odd
     73    """
     74    y=k;to_add=k
     75    for i in xrange((m-1)//2):
     76        to_add=to_add**4
     77        y+=to_add
     78    return y
     79
     80def T4_calc(k,m):
     81    """
     82    Computes Trace_F/GF(4)_(k), as well as the solution for m=2mod4, and other numbers that we will need for the computation with m=0mod4
     83    """
     84    if m.mod(4)==0:
     85        T4=k;to_add=k; #print '0mod4'
     86        for i in xrange((m-2)//2):
     87            to_add=to_add**4
     88            if 2*(i+1)==m-2:
     89                squ=to_add
     90            if i==(m//4)-1:
     91                hfwy=to_add
     92            T4+=to_add
     93        to_add=to_add**2
     94        return [T4,squ,hfwy]
     95    else: #m.mod(4)==2
     96        T4=k;sol=k;to_add=k; #print '2mod4'
     97        for i in xrange((m-6)//4):
     98            to_add=to_add**4
     99            T4+=to_add
     100            to_add=to_add**4
     101            T4+=to_add
     102            sol+=to_add
     103        for i in [1,2]:
     104            to_add=to_add**4
     105            T4+=to_add
     106        return [T4,sol]
     107
     108def alp(F):
     109    """
     110    Computes an alp satisfying: alp^2+alp=1, for use in twomod_t1. This part is slow.
     111    """
     112    e=(F.order()-1)/3
     113    alp=(F.gen())**e
     114    while not(alp**2+alp==1):    #probably <3 times
     115        a=F.random_element()
     116        alp=a**e
     117    return alp
     118
     119def zmod_t1(k,m,squ,hfwy):
     120    """
     121    Does the compution for m=0mod4, Tr_4(k)=1. The case m=0mod4, Tr_4(k)=0 is reduces to this one by doquad.
     122    """
     123    m4=m//4
     124    y=hfwy;to_add=hfwy     #;print 'zmod_t1'
     125    for i in xrange(m4-2):
     126        to_add=to_add**4
     127        y+=to_add
     128    s=k*y**2
     129    y+=to_add**4
     130    to_times=k
     131    to_add_start=hfwy**2
     132    for j in xrange(1,m4-1):
     133        to_times=to_times**4
     134        to_add_start=to_add_start**4
     135        to_add=to_add_start
     136        s_temp=to_add
     137        for i in xrange(j,m4-2):
     138            to_add=to_add**4
     139            s_temp+=to_add
     140        s+=s_temp*to_times
     141    y=s+s**2+squ**2*(1+y) #squ from T4
     142    return y
     143
     144def doquad(k,m,F):
     145    """
     146    Divides the task of factoring into cases, depending on whether m is even or odd, and on the trace down to GF(4).
     147    """
     148    if m%2==1:
     149        return qodd(k,m)     #;print 'odd'
     150    else: #m is even
     151        if m.mod(4)==2:
     152            T4,sol=T4_calc(k,m)
     153            if T4==0:
     154                return sol**4+sol**8     #; print 't4=0,2mod4'
     155            else: #T4==1
     156                return (sol**4+sol**8)+alp(F)   #  ; print 't4=1,2mod4'
     157        else: #m = 0 mod4
     158            T4,squ,hfwy=T4_calc(k,m)
     159            if T4==1:
     160                if m==4:
     161                    return (k**8+k**12)
     162                else:
     163                    return zmod_t1(k,m,squ,hfwy)
     164            else: #T4=0:
     165                t=F.random_element()    #;print 't4=0,0mod4'
     166                while t.trace()==0: #we want trace=1, about 50%
     167                    t=F.random_element()
     168                k1=t+t**2
     169                y1=doquad(k+k1,m,F) #we've reduced to another case
     170                return (t+y1)
     171
  • sage/rings/polynomial/polynomial_element.pyx

    diff -r 2e793d2a0e12 -r 28635b3d0b70 sage/rings/polynomial/polynomial_element.pyx
    a b  
    24072407            Traceback (most recent call last):
    24082408            ...
    24092409            NotImplementedError: factorization of polynomials over rings with composite characteristic is not implemented
     2410
     2411        For quadratrics over GF(2^m) ::
     2412
     2413            sage: F.<x> = GF(4)
     2414            sage: PP = PolynomialRing(F,'y')
     2415            sage: f = PP([1,1,1]);f
     2416            y^2 + y + 1
     2417            sage: f.factor()
     2418            (y + x) * (y + x + 1)
    24102419        """
    24112420
    24122421        # PERFORMANCE NOTE:
     
    24582467        from sage.rings.integer_mod_ring import is_IntegerModRing
    24592468        from sage.rings.integer_ring import is_IntegerRing
    24602469        from sage.rings.rational_field import is_RationalField
     2470        from sage.rings.polynomial.char_two_quad import quad
    24612471
    24622472        n = None
    24632473        if is_IntegerModRing(R) or is_IntegerRing(R) or is_RationalField(R):
     
    24782488            return Factorization(v, from_M(F.unit()))
    24792489
    24802490        elif is_FiniteField(R):
    2481             v = [x._pari_("a") for x in self.list()]
    2482             f = pari(v).Polrev()
    2483             G = list(f.factor())
     2491            if ch==2 and self.degree()==2:
     2492                G=quad(self)
     2493                return Factorization(G[0],G[1])
     2494            else:
     2495                v = [x._pari_("a") for x in self.list()]
     2496                f = pari(v).Polrev()
     2497                G = list(f.factor())
    24842498
    24852499           
    24862500        elif is_NumberField(R):
  • new file sage/rings/polynomial/char_two_quad.py

    # HG changeset patch
    # User Preston Wake <preston.wake@gmail.com>
    # Date 1248373201 14400
    # Node ID 28635b3d0b70e166aab75e0492f397f051df35ce
    # Parent  2e793d2a0e123293b73eed40715e43185fd9ccfe
    added speciallized code to factor quadratic polynomials over GF(2^m)
    
    diff -r 2e793d2a0e12 -r 28635b3d0b70 sage/rings/polynomial/char_two_quad.py
    - +  
     1def quad(f):
     2    r"""
     3    Returns the factorization of f, for f a quadratic polynomial over GF(2^m).
     4    Taken from CHIN-LONG CHEN's 1982 paper Formulas for the Solutions of Quadratic Equations over GF (2^m).
     5    By Preston Wake and Wouter Castryck, 2009
     6
     7    EXAMPLES:
     8
     9    ::
     10        sage: from sage.rings.polynomial.char_two_quad import quad
     11        sage: F.<x> = GF(4)
     12        sage: PP = PolynomialRing(F,'y')
     13        sage: f = PP([1,1,1]);f
     14        y^2 + y + 1
     15        sage: a = quad(f); a
     16        ([(y + x, 1), (y + x + 1, 1)], 1)
     17        sage: r=a[0][0][0];r
     18        y + x
     19        sage: f(r(0))
     20        0
     21        sage: f/r
     22        y + x + 1
     23   
     24    One where f is a square
     25
     26    ::
     27        sage: f = PP([1,0,1]);f
     28        y^2 + 1
     29        sage: a = quad(f); a
     30        ([(y + 1, 2)], 1)
     31   
     32    One where f is irreducible
     33   
     34    ::
     35        sage: f = PP([x,1,1]);f
     36        y^2 + y + x
     37        sage: a = quad(f);a
     38        ([(y^2 + y + x, 1)], 1)
     39           
     40    """
     41    PP=f.parent()
     42    F=PP.base_ring()
     43    m=F.degree();
     44    a=f[2];b=f[1];c=f[0];
     45    if a==0:
     46        raise ValueError, 'not a quad'
     47    if b==0:
     48        y=gftwosqrt(c/a)
     49        f1=PP([y,1])
     50        return ([(f1,2)],a)
     51    k=a*c/b**2
     52    if not(k.trace()==0):
     53        return ([(f*(a**(-1)),1)],a)
     54    else:
     55        y=doquad(k,m,F)
     56    y1=y+1
     57    beta=b/a
     58    f1=PP([y*beta,1]);f2=PP([y1*beta,1])
     59    return ([(f1,1),(f2,1)],a)
     60
     61
     62def gftwosqrt(r):
     63    """
     64    Quickly finds the squareroot of an element in GF(2^m)
     65    """
     66    F=r.parent()
     67    c=F.cardinality()
     68    return r**(c/2)
     69
     70def qodd(k,m):
     71    """
     72    Does the compution for m odd
     73    """
     74    y=k;to_add=k
     75    for i in xrange((m-1)//2):
     76        to_add=to_add**4
     77        y+=to_add
     78    return y
     79
     80def T4_calc(k,m):
     81    """
     82    Computes Trace_F/GF(4)_(k), as well as the solution for m=2mod4, and other numbers that we will need for the computation with m=0mod4
     83    """
     84    if m.mod(4)==0:
     85        T4=k;to_add=k; #print '0mod4'
     86        for i in xrange((m-2)//2):
     87            to_add=to_add**4
     88            if 2*(i+1)==m-2:
     89                squ=to_add
     90            if i==(m//4)-1:
     91                hfwy=to_add
     92            T4+=to_add
     93        to_add=to_add**2
     94        return [T4,squ,hfwy]
     95    else: #m.mod(4)==2
     96        T4=k;sol=k;to_add=k; #print '2mod4'
     97        for i in xrange((m-6)//4):
     98            to_add=to_add**4
     99            T4+=to_add
     100            to_add=to_add**4
     101            T4+=to_add
     102            sol+=to_add
     103        for i in [1,2]:
     104            to_add=to_add**4
     105            T4+=to_add
     106        return [T4,sol]
     107
     108def alp(F):
     109    """
     110    Computes an alp satisfying: alp^2+alp=1, for use in twomod_t1. This part is slow.
     111    """
     112    e=(F.order()-1)/3
     113    alp=(F.gen())**e
     114    while not(alp**2+alp==1):    #probably <3 times
     115        a=F.random_element()
     116        alp=a**e
     117    return alp
     118
     119def zmod_t1(k,m,squ,hfwy):
     120    """
     121    Does the compution for m=0mod4, Tr_4(k)=1. The case m=0mod4, Tr_4(k)=0 is reduces to this one by doquad.
     122    """
     123    m4=m//4
     124    y=hfwy;to_add=hfwy     #;print 'zmod_t1'
     125    for i in xrange(m4-2):
     126        to_add=to_add**4
     127        y+=to_add
     128    s=k*y**2
     129    y+=to_add**4
     130    to_times=k
     131    to_add_start=hfwy**2
     132    for j in xrange(1,m4-1):
     133        to_times=to_times**4
     134        to_add_start=to_add_start**4
     135        to_add=to_add_start
     136        s_temp=to_add
     137        for i in xrange(j,m4-2):
     138            to_add=to_add**4
     139            s_temp+=to_add
     140        s+=s_temp*to_times
     141    y=s+s**2+squ**2*(1+y) #squ from T4
     142    return y
     143
     144def doquad(k,m,F):
     145    """
     146    Divides the task of factoring into cases, depending on whether m is even or odd, and on the trace down to GF(4).
     147    """
     148    if m%2==1:
     149        return qodd(k,m)     #;print 'odd'
     150    else: #m is even
     151        if m.mod(4)==2:
     152            T4,sol=T4_calc(k,m)
     153            if T4==0:
     154                return sol**4+sol**8     #; print 't4=0,2mod4'
     155            else: #T4==1
     156                return (sol**4+sol**8)+alp(F)   #  ; print 't4=1,2mod4'
     157        else: #m = 0 mod4
     158            T4,squ,hfwy=T4_calc(k,m)
     159            if T4==1:
     160                if m==4:
     161                    return (k**8+k**12)
     162                else:
     163                    return zmod_t1(k,m,squ,hfwy)
     164            else: #T4=0:
     165                t=F.random_element()    #;print 't4=0,0mod4'
     166                while t.trace()==0: #we want trace=1, about 50%
     167                    t=F.random_element()
     168                k1=t+t**2
     169                y1=doquad(k+k1,m,F) #we've reduced to another case
     170                return (t+y1)
     171
  • sage/rings/polynomial/polynomial_element.pyx

    diff -r 2e793d2a0e12 -r 28635b3d0b70 sage/rings/polynomial/polynomial_element.pyx
    a b  
    24072407            Traceback (most recent call last):
    24082408            ...
    24092409            NotImplementedError: factorization of polynomials over rings with composite characteristic is not implemented
     2410
     2411        For quadratrics over GF(2^m) ::
     2412
     2413            sage: F.<x> = GF(4)
     2414            sage: PP = PolynomialRing(F,'y')
     2415            sage: f = PP([1,1,1]);f
     2416            y^2 + y + 1
     2417            sage: f.factor()
     2418            (y + x) * (y + x + 1)
    24102419        """
    24112420
    24122421        # PERFORMANCE NOTE:
     
    24582467        from sage.rings.integer_mod_ring import is_IntegerModRing
    24592468        from sage.rings.integer_ring import is_IntegerRing
    24602469        from sage.rings.rational_field import is_RationalField
     2470        from sage.rings.polynomial.char_two_quad import quad
    24612471
    24622472        n = None
    24632473        if is_IntegerModRing(R) or is_IntegerRing(R) or is_RationalField(R):
     
    24782488            return Factorization(v, from_M(F.unit()))
    24792489
    24802490        elif is_FiniteField(R):
    2481             v = [x._pari_("a") for x in self.list()]
    2482             f = pari(v).Polrev()
    2483             G = list(f.factor())
     2491            if ch==2 and self.degree()==2:
     2492                G=quad(self)
     2493                return Factorization(G[0],G[1])
     2494            else:
     2495                v = [x._pari_("a") for x in self.list()]
     2496                f = pari(v).Polrev()
     2497                G = list(f.factor())
    24842498
    24852499           
    24862500        elif is_NumberField(R):
  • sage/modular/hecke/algebra.py

    # HG changeset patch
    # User Preston Wake <preston.wake@gmail.com>
    # Date 1250471218 14400
    # Node ID 0fce5ebf9305c04d7c997245e2608ed25c761f0d
    # Parent  2e793d2a0e123293b73eed40715e43185fd9ccfe
    Added code for bases of Hecke Algebras
    
    diff -r 2e793d2a0e12 -r 0fce5ebf9305 sage/modular/hecke/algebra.py
    a b  
    3131import math
    3232import weakref
    3333
     34import sage.rings.all as rings
    3435import sage.rings.arith as arith
    3536import sage.rings.infinity
    3637import sage.misc.latex as latex
     
    4041import sage.rings.commutative_algebra
    4142from sage.misc.misc import verbose
    4243from sage.matrix.constructor import matrix
     44from sage.rings.arith import lcm
     45from sage.matrix.matrix_space import MatrixSpace
    4346
    4447def is_HeckeAlgebra(x):
    4548    r"""
     
    142145    return T
    143146
    144147
     148def _heckebasis(M):
     149    r"""
     150    Gives a basis of the hecke algebra of M as a ZZ-module
     151
     152    INPUT:
     153
     154    - ``M`` - a hecke module
     155
     156    OUTPUT:
     157
     158    - a list of hecke algebra elements represented as matrices
     159
     160    EXAMPLES::
     161   
     162        sage: M = ModularSymbols(11,2,1)
     163        sage: sage.modular.hecke.algebra._heckebasis(M)
     164        [Hecke operator on Modular Symbols space of dimension 2 for Gamma_0(11) of weight 2 with sign 1 over Rational Field defined by:
     165        [1 0]
     166        [0 1],
     167        Hecke operator on Modular Symbols space of dimension 2 for Gamma_0(11) of weight 2 with sign 1 over Rational Field defined by:
     168        [0 1]
     169        [0 5]]
     170    """
     171    QQ = rings.QQ
     172    ZZ = rings.ZZ
     173    d = M.rank()
     174    VV = QQ**(d**2)
     175    WW = ZZ**(d**2)
     176    MM = MatrixSpace(QQ,d)
     177    MMZ = MatrixSpace(ZZ,d)
     178    S = []; Denom = []; B = []; B1 = []
     179    for i in xrange(1, M.hecke_bound() + 1):
     180        v = M.hecke_operator(i).matrix()
     181        den = v.denominator()
     182        Denom.append(den)
     183        S.append(v)
     184    den = lcm(Denom)
     185    for m in S:
     186        B.append(WW((den*m).list()))
     187    UU = WW.submodule(B)
     188    B = UU.basis()
     189    for u in B:
     190        u1 = u.list()
     191        m1 = M.hecke_algebra()(MM(u1), check=False)
     192        #m1 = MM(u1)
     193        B1.append((1/den)*m1)
     194    return B1
     195
     196
    145197class HeckeAlgebra_base(sage.rings.commutative_algebra.CommutativeAlgebra):
    146198    """
    147199    Base class for algebras of Hecke operators on a fixed Hecke module.
     
    404456    def basis(self):
    405457        r"""
    406458        Return a basis for this Hecke algebra as a free module over
    407         its base ring. Not implemented at present.
     459        its base ring.
    408460
    409461        EXAMPLE::
    410462
    411463            sage: ModularSymbols(Gamma1(3), 3).hecke_algebra().basis()
    412             Traceback (most recent call last):
    413             ...
    414             NotImplementedError
     464            [Hecke operator on Modular Symbols space of dimension 2 for Gamma_1(3) of weight 3 with sign 0 and over Rational Field defined by:
     465            [1 0]
     466            [0 1],
     467            Hecke operator on Modular Symbols space of dimension 2 for Gamma_1(3) of weight 3 with sign 0 and over Rational Field defined by:
     468            [0 0]
     469            [0 2]]
    415470        """
    416         raise NotImplementedError
     471        try:
     472            return self.__basis_cache
     473        except AttributeError:
     474            self.__basis_cache=_heckebasis(self.__M)
     475            return self.__basis_cache
    417476
    418477    def discriminant(self):
    419478        r"""