Ticket #4302: 4302_speedup1.patch

File 4302_speedup1.patch, 2.2 KB (added by zimmerma, 11 years ago)
  • sage/rings/polynomial/polynomial_gf2x.pyx

    # HG changeset patch
    # User Paul Zimmermann <zimmerma@loria.fr>
    # Date 1224348127 -7200
    # Node ID a4f95b4f41101f203e14ebc4a6a3bc069e9763c8
    # Parent  3be92d723a616f423472372f4ee1a5846da93821
    patch for #4302: replace half mults by squarings for computing G matrix
    
    diff -r 3be92d723a61 -r a4f95b4f4110 sage/rings/polynomial/polynomial_gf2x.pyx
    a b  
    122122
    123123        cdef GF2X_c _f = (<Polynomial_GF2X>self).x
    124124        cdef GF2X_c _g = (<Polynomial_GF2X>g).x
    125         cdef GF2X_c gpow, tt
     125        cdef GF2X_c gpow, g2, tt
    126126        GF2X_conv_long(gpow, 1)
    127127
    128128        maxlength = GF2X_NumBits(_f)
     
    138138        G = <Matrix_mod2_dense>Matrix(GF(2), k, n)
    139139
    140140        # first compute g^j mod h, 2 <= j < k
    141         for j in range(0, k):
     141        # first deal with j=0
     142        for i from 0 <= i < GF2X_NumBits(gpow):
     143            mzd_write_bit(G._entries, 0, i, GF2_conv_to_long(GF2X_coeff(gpow, i)))
     144        # precompute g^2
     145        GF2X_SqrMod_pre(g2, _g, modulus)
     146        gpow = _g
     147        for j in range(1, k, 2):
     148            if j > 1:
     149                GF2X_MulMod_pre(gpow, gpow, g2, modulus) # gpow = g^j
    142150            for i from 0 <= i < GF2X_NumBits(gpow):
    143151                mzd_write_bit(G._entries, j, i, GF2_conv_to_long(GF2X_coeff(gpow, i)))
    144             #gpow = (gpow * g) % h # we'll need g^k below
     152            # we now process 2j, 4j, 8j, ... by squaring each time
     153            if 2*j < k:
     154                tt = gpow
     155                jj = j
     156                while 2*jj < k:
     157                   GF2X_SqrMod_pre(tt, tt, modulus)
     158                   jj = 2*jj
     159                   for i from 0 <= i < GF2X_NumBits(tt):
     160                       mzd_write_bit(G._entries, jj, i, GF2_conv_to_long(GF2X_coeff(tt, i)))
     161        # we need that gpow = g^k at the end
     162        if k % 2 == 1: # k is odd, last j is k-2
     163            GF2X_MulMod_pre(gpow, gpow, g2, modulus)
     164        else:          # k is even, last j is k-1
    145165            GF2X_MulMod_pre(gpow, gpow, _g, modulus)
    146166        verbose("G %d x %d %5.3f s"%(G.nrows(), G.ncols(),cputime(t)),level=1)
    147167