Changeset 7764:cd445dcb4bb7


Ignore:
Timestamp:
11/18/07 23:10:14 (6 years ago)
Author:
Robert L Miller <rlm@…>
Branch:
default
Message:

5th half - fixed many bugs, print a lot of crap: code 111100\001111 now gets correct aut group

Location:
sage/coding
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • sage/coding/binary_code.pxd

    r7763 r7764  
    5656    cdef int sat_225(self, int) 
    5757    cdef int min_cell_reps(self, int) 
     58    cdef void new_min_cell_reps(self, int, int *, int) 
    5859    cdef int fixed_cols(self, int, int) 
     60    cdef void fixed_vertices(self, int, int *, int *, int) 
    5961    cdef int first_smallest_nontrivial(self, int) 
    6062    cdef int new_first_smallest_nontrivial(self, int, int *, int) 

  • sage/coding/binary_code.pyx

    r7763 r7764  
    740740            for k from 0 <= k < nwords-1: 
    741741                wd_ents[k] = k 
    742                 wd_lvls[k] = nwords 
     742                wd_lvls[k] = 2*ncols 
    743743            for k from 0 <= k < ncols-1: 
    744744                col_ents[k] = k 
    745                 col_lvls[k] = ncols 
     745                col_lvls[k] = 2*ncols 
    746746            wd_ents[nwords-1] = nwords-1 
    747747            wd_lvls[nwords-1] = -1 
    … …  
    905905        cdef int i, j, k 
    906906        s = '' 
    907         for i from 0 <= i < self.nwords: 
    908             s += '({' 
    909             for j from 0 <= j < self.nwords: 
    910                 s += str(self.wd_ents[j]) 
    911                 if self.wd_lvls[j] <= i: 
    912                     s += '},{' 
    913                 else: 
    914                     s += ',' 
    915             s = s[:-2] + ')\n' 
    916         for i from 0 <= i < self.ncols: 
    917             s += '({' 
    918             for j from 0 <= j < self.ncols: 
    919                 s += str(self.col_ents[j]) 
    920                 if self.col_lvls[j] <= i: 
    921                     s += '},{' 
    922                 else: 
    923                     s += ',' 
    924             s = s[:-2] + ')\n' 
     907        for i from 0 <= i < 2*self.ncols: 
     908            if i == 0 or not self.is_discrete(i-1): 
     909                s += '({' 
     910                for j from 0 <= j < self.nwords: 
     911                    s += str(self.wd_ents[j]) 
     912                    if self.wd_lvls[j] <= i: 
     913                        s += '},{' 
     914                    else: 
     915                        s += ',' 
     916                s = s[:-2] + ')  ' 
     917                s += '({' 
     918                for j from 0 <= j < self.ncols: 
     919                    s += str(self.col_ents[j]) 
     920                    if self.col_lvls[j] <= i: 
     921                        s += '},{' 
     922                    else: 
     923                        s += ',' 
     924                s = s[:-2] + ')\n' 
    925925        return s 
    926926 
    … …  
    10981098        return reps 
    10991099 
     1100    cdef void new_min_cell_reps(self, int k, int *Omega, int start): 
     1101        cdef int i, j 
     1102        cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls 
     1103        cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents 
     1104        cdef int reps = (1 << self_col_ents[0]) 
     1105        cdef int radix = self.radix, nwords = self.nwords, ncols = self.ncols 
     1106        for i from 0 < i < ncols: 
     1107            if self_col_lvls[i-1] <= k: 
     1108                reps += (1 << self_col_ents[i]) 
     1109        Omega[start] = reps 
     1110        reps = 1 
     1111        for i from 0 < i < min(radix, nwords): 
     1112            if self_wd_lvls[i-1] <= k: 
     1113                reps += (1 << self_wd_ents[i]) 
     1114        Omega[start+1] = reps 
     1115        j = radix 
     1116        while j < nwords: 
     1117            reps = 0 
     1118            for i from 0 <= i < min(radix, nwords - j): 
     1119                if self_wd_lvls[j + i - 1] <= k: 
     1120                    reps += (1 << self_wd_ents[i]) 
     1121            Omega[start+1+j] = reps 
     1122            j += radix 
     1123 
    11001124    def _fixed_cols(self, mcrs, k): #TODO 
    11011125        """ 
    … …  
    11351159                fixed += (1 << i) 
    11361160        return fixed & mcrs 
     1161 
     1162    cdef void fixed_vertices(self, int k, int *Phi, int *Omega, int start): 
     1163        cdef int i, j, ell 
     1164        cdef int fixed = 0, ncols = self.ncols, nwords = self.nwords 
     1165        cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls 
     1166        cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents 
     1167        for i from 0 <= i < ncols: 
     1168            if self_col_lvls[i] <= k: 
     1169                fixed += (1 << self_col_ents[i]) 
     1170        Phi[start] = fixed & Omega[start] 
     1171        # zero out the rest of Phi 
     1172        ell = 1 + nwords/self.radix 
     1173        if nwords%self.radix: 
     1174            ell += 1 
     1175        for i from 0 < i < ell: 
     1176            Phi[start+i] = 0 
     1177        for i from 0 <= i < nwords: 
     1178            if self_wd_lvls[i] <= k: 
     1179                ell = self_wd_ents[i] 
     1180                j =   ell/self.radix 
     1181                ell = ell%self.radix 
     1182                if Omega[start+1+j]&(1 << ell): 
     1183                    Phi[start+1+j] ^= (1 << ell) 
    11371184 
    11381185    def _first_smallest_nontrivial(self, k): #TODO 
    … …  
    11841231            j += 1 
    11851232        # j now points to the last element of the cell 
    1186         print "fsnt:", location, j-location+1 
     1233#        print "fsnt:", location, j-location+1 
    11871234        i = self.radix - j - 1                 # the cell is represented in binary, reading from the right: 
    11881235        cell = (~0 << location) ^ (~0 << j+1)  # <-------            self.radix               -----> 
    … …  
    11901237 
    11911238    cdef int new_first_smallest_nontrivial(self, int k, int *W, int start): 
    1192         cdef int cell 
     1239        cdef int ell 
    11931240        cdef int i = 0, j = 0, location = 0, min = self.ncols, nwords = self.nwords 
    11941241        cdef int min_is_col = 1, radix = self.radix 
    11951242        cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls 
     1243        cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents 
    11961244        while True: 
    11971245            if self_col_lvls[i] <= k: 
    … …  
    12081256                    min = i - j + 1 
    12091257                    min_is_col = 0 
    1210                     location = k 
     1258                    location = j 
    12111259                j = i + 1 
    12121260            if self_wd_lvls[i] == -1: break 
    … …  
    12151263        # nontrivial cell 
    12161264        j = location 
     1265        #zero out this level of W: 
     1266        ell = 1 + nwords/radix 
     1267        if nwords%radix: 
     1268            ell += 1 
     1269        for i from 0 <= i < ell: 
     1270            W[start+i] = 0 
    12171271        if min_is_col: 
    12181272            while True: 
    … …  
    12201274                j += 1 
    12211275            # j now points to the last element of the cell 
    1222             W[start] = (~0 << location) ^ (~0 << (j + 1)) 
    1223             return self.col_ents[location] 
     1276            i = location 
     1277            while i <= j: 
     1278                W[start] ^= (1 << self_col_ents[i]) 
     1279                i += 1 
     1280            return self_col_ents[location] 
    12241281        else: 
    12251282            while True: 
    … …  
    12271284                j += 1 
    12281285            # j now points to the last element of the cell 
    1229             i = 0 
    1230             while location >= radix*i + radix: 
     1286            i = location 
     1287            while i <= j: 
     1288                ell = self_wd_ents[i] 
     1289                W[start+1+ell/radix] ^= (1 << ell%radix) 
    12311290                i += 1 
    1232             if j < radix*i + radix: 
    1233                 W[1 + i] = (~0 << (location - radix*i)) ^ (~0 << (j - radix*i + 1)) 
    1234             else: # it will span at most two ints: 
    1235                 W[1 + i] =      (~0 << (location - radix*i)) 
    1236                 W[1 + i + 1] = ~(~0 << (j - radix*(i+1) + 1)) 
    1237             return self.wd_ents[location] ^ self.flag 
     1291            return self_wd_ents[location] ^ self.flag 
    12381292 
    12391293    def _dangerous_dont_use_set_ents_lvls(self, col_ents, col_lvls, wd_ents, wd_lvls): 
    … …  
    14321486 
    14331487    cdef int split_vertex(self, int v, int k): 
    1434         cdef int i = 0, j 
     1488        cdef int i = 0, j, flag = self.flag 
    14351489        cdef int *ents 
    14361490        cdef int *lvls 
    1437         if v & self.flag: 
     1491        if v & flag: 
    14381492            ents = self.wd_ents 
    14391493            lvls = self.wd_lvls 
     1494            v = v ^ flag 
     1495            while ents[i] != v: i += 1 
     1496            v = v ^ flag 
    14401497        else: 
    14411498            ents = self.col_ents 
    14421499            lvls = self.col_lvls 
    1443         while ents[i] != v: i += 1 
     1500            while ents[i] != v: i += 1 
    14441501        j = i 
    14451502        while lvls[i] > k: i += 1 
    … …  
    14531510                ents[j] = ents[j-1] 
    14541511                j -= 1 
    1455             ents[j] = v 
     1512            if v & flag: 
     1513                ents[j] = v ^ flag 
     1514            else: 
     1515                ents[j] = v 
    14561516        lvls[j] = k 
    1457         return j 
     1517        if v & flag: 
     1518            return j ^ flag 
     1519        else: 
     1520            return j 
    14581521 
    14591522    def _col_degree(self, C, col, wd_ptr, k): 
    … …  
    14941557        cdef int i = 0 
    14951558        cdef int *self_wd_lvls = self.wd_lvls, *self_wd_ents = self.wd_ents 
    1496         col = self.col_ents[col] 
    14971559        while True: 
    14981560            if CG.is_one(self_wd_ents[wd_ptr], col): i += 1 
    … …  
    15401602    cdef int wd_degree(self, BinaryCode CG, int wd, int col_ptr, int k, int *ham_wts): 
    15411603 
    1542         cdef int mask = (1 << col_ptr) 
    1543         cdef int *self_col_lvls = self.col_lvls 
     1604        cdef int *self_col_lvls = self.col_lvls, *self_col_ents = self.col_ents 
     1605        cdef int mask = (1 << self_col_ents[col_ptr]) 
    15441606        while self_col_lvls[col_ptr] > k: 
    15451607            col_ptr += 1 
    1546             mask += (1 << col_ptr) 
    1547         mask &= CG.words[self.wd_ents[wd]] 
     1608            mask += (1 << self_col_ents[col_ptr]) 
     1609        mask &= CG.words[wd] 
    15481610        return ham_wts[mask & 65535] + ham_wts[(mask >> 16) & 65535] 
    15491611 
    … …  
    17891851        cdef int q, r, s, t, flag = self.flag, self_ncols = self.ncols 
    17901852        cdef int t_w, self_nwords = self.nwords, invariant = 0, i, j, m = 0 
    1791         cdef int *self_wd_degs = self.wd_degs, *self_wd_lvls = self.wd_lvls, *self_col_lvls = self.col_lvls 
    1792         cdef int *self_col_degs = self.col_degs 
     1853        cdef int *self_wd_degs = self.wd_degs, *self_wd_lvls = self.wd_lvls, *self_wd_ents = self.wd_ents 
     1854        cdef int *self_col_degs = self.col_degs, *self_col_lvls = self.col_lvls, *self_col_ents = self.col_ents 
    17931855        while not self.is_discrete(k) and m < alpha_length: 
    1794 #            print "m:", m 
    1795 #            print "alpha:", ','.join(['w'+str(alpha[i]^flag) if alpha[i]&flag else 'c'+str(alpha[i]) for i from 0 <= i < alpha_length]) 
     1856            print "m:", m 
     1857            print "alpha:", ','.join(['w'+str(alpha[i]^flag) if alpha[i]&flag else 'c'+str(alpha[i]) for i from 0 <= i < alpha_length]) 
     1858            print self 
    17961859            invariant += 1 
    17971860            j = 0 
    17981861            if alpha[m] & flag: 
     1862#                print 'word' 
    17991863                while j < self_ncols: 
     1864#                    print 'j', j 
     1865#                    print self 
    18001866                    i = j; s = 0 
    18011867                    invariant += 8 
    18021868                    while True: 
    1803                         self_col_degs[i-j] = self.col_degree(CG, i, alpha[m]^flag, k) 
     1869#                        print 'col_i', self_col_ents[i] 
     1870#                        print 'alpha[m]^flag', alpha[m]^flag 
     1871                        self_col_degs[i-j] = self.col_degree(CG, self_col_ents[i], alpha[m]^flag, k) 
     1872#                        print 'deg', self_col_degs[i-j] 
    18041873                        if s == 0 and self_col_degs[i-j] != self_col_degs[0]: s = 1 
    18051874                        i += 1 
    18061875                        if self_col_lvls[i-1] <= k: break 
    18071876                    if s: 
     1877#                        print 's' 
    18081878                        invariant += 8 
    18091879                        t = self.sort_cols(j, k) 
    … …  
    18261896                            j += 1 
    18271897                        j += 1 
     1898#                        print 'end s j:', j 
    18281899                        invariant += (i-j) 
    18291900                    else: j = i 
    18301901            else: 
     1902                print 'col' 
    18311903                while j < self.nwords: 
     1904                    print 'j', j 
     1905                    print self 
    18321906                    i = j; s = 0 
    18331907                    invariant += 64 
    18341908                    while True: 
    1835                         self_wd_degs[i-j] = self.wd_degree(CG, i, alpha[m], k, ham_wts) 
     1909                        print 'i', i 
     1910                        self_wd_degs[i-j] = self.wd_degree(CG, self_wd_ents[i], alpha[m], k, ham_wts) 
     1911                        print 'deg', self_wd_degs[i-j] 
    18361912                        if s == 0 and self_wd_degs[i-j] != self_wd_degs[0]: s = 1 
    18371913                        i += 1 
    … …  
    21802256        self.w_gamma =     <int *> sage_malloc( self.w_gamma_size              * sizeof(int) ) 
    21812257        self.alpha =       <int *> sage_malloc( self.alpha_size                * sizeof(int) ) 
    2182         self.new_Phi =     <int *> sage_malloc( self.Phi_size * self.L         * sizeof(int) ) 
    2183         self.new_Omega =   <int *> sage_malloc( self.Phi_size * self.L         * sizeof(int) ) 
    2184         self.new_W =       <int *> sage_malloc( self.Phi_size * self.radix * 2 * sizeof(int) ) 
     2258        self.Phi =     <int *> sage_malloc( self.Phi_size * (self.L+1)     * sizeof(int) ) 
     2259        self.Omega =   <int *> sage_malloc( self.Phi_size * self.L         * sizeof(int) ) 
     2260        self.W =       <int *> sage_malloc( self.Phi_size * self.radix * 2 * sizeof(int) ) 
    21852261 
    21862262        self.aut_gp_gens = <int *> sage_malloc( self.aut_gens_size             * sizeof(int) ) 
    … …  
    21932269        self.e =           <int *> sage_malloc( self.radix * 2                 * sizeof(int) ) 
    21942270 
    2195         # TODO: remove these 
    2196         self.Phi =         <int *> sage_malloc( self.L                         * sizeof(int) ) 
    2197         self.Omega =       <int *> sage_malloc( self.L                         * sizeof(int) ) 
    2198         self.W =           <int *> sage_malloc( self.radix                     * sizeof(int) ) 
    2199  
    22002271        if not (self.Phi and self.Omega and self.W and self.Lambda1 and self.Lambda2 and self.Lambda3 \ 
    22012272            and self.w_gamma and self.c_gamma and self.alpha and self.v and self.e and self.aut_gp_gens \ 
    2202             and self.labeling and self.new_Phi and self.new_Omega and self.new_W): 
     2273            and self.labeling): 
    22032274            if self.Phi:          sage_free(self.Phi) 
    22042275            if self.Omega:        sage_free(self.Omega) 
    22052276            if self.W:            sage_free(self.W) 
    2206             if self.new_Phi:      sage_free(self.new_Phi) 
    2207             if self.new_Omega:    sage_free(self.new_Omega) 
    2208             if self.new_W:        sage_free(self.new_W) 
    22092277            if self.Lambda1:      sage_free(self.Lambda1) 
    22102278            if self.Lambda2:      sage_free(self.Lambda2) 
    … …  
    22312299        if self.w_gamma:   sage_free(self.w_gamma) 
    22322300        if self.alpha:     sage_free(self.alpha) 
    2233         if self.new_Phi:   sage_free(self.new_Phi) 
    2234         if self.new_Omega: sage_free(self.new_Omega) 
    2235         if self.new_W:     sage_free(self.new_W) 
    22362301 
    22372302        if self.v:           sage_free(self.v) 
    … …  
    22692334 
    22702335        # declare variables: 
    2271         cdef int i, j, ii, jj # local variables 
     2336        cdef int i, j, ii, jj, iii, jjj # local variables 
    22722337 
    22732338        cdef PartitionStack nu, zeta, rho # nu is the current position in the tree, 
    … …  
    22912356        cdef OrbitPartition Theta # keeps track of which vertices have been discovered to be equivalent 
    22922357        cdef int *Phi   = self.Phi      # Phi stores the fixed point sets of each automorphism 
    2293         cdef int *new_Phi = self.new_Phi 
    22942358        cdef int *Omega = self.Omega    # Omega stores the minimal elements of each cell of the orbit partition 
    2295         cdef int *new_Omega = self.new_Omega 
    22962359        cdef int l = -1                 # current index for storing values in Phi and Omega- we start at -1 so that when 
    22972360                                        # we increment first, the first place we write to is 0. 
    … …  
    22992362                                # the current partition, at k. Phi and Omega are ultimately used to make the size of 
    23002363                                # W as small as possible 
    2301         cdef int *new_W = self.new_W 
    23022364        cdef int *e = self.e    # 0 or 1, whether or not we have used Omega and Phi to narrow down W[k] yet: see states 12 and 17 
    23032365 
    … …  
    23262388            self.alpha_size = self.w_gamma_size + self.radix 
    23272389            self.Phi_size = self.w_gamma_size/self.radix + 1 
    2328             self.w_gamma =     <int *> sage_realloc(self.w_gamma,   self.w_gamma_size              * sizeof(int) ) 
    2329             self.alpha =       <int *> sage_realloc(self.alpha,     self.alpha_size                * sizeof(int) ) 
    2330             self.new_Phi =     <int *> sage_realloc(self.new_Phi,   self.Phi_size * self.L         * sizeof(int) ) 
    2331             self.new_Omega =   <int *> sage_realloc(self.new_Omega, self.Phi_size * self.L         * sizeof(int) ) 
    2332             self.new_W =       <int *> sage_realloc(self.new_W,     self.Phi_size * self.radix * 2 * sizeof(int) ) 
    2333             if not (self.w_gamma and self.alpha and self.new_Phi and self.new_Omega and self.new_W): 
     2390            self.w_gamma = <int *> sage_realloc(self.w_gamma,   self.w_gamma_size              * sizeof(int) ) 
     2391            self.alpha =   <int *> sage_realloc(self.alpha,     self.alpha_size                * sizeof(int) ) 
     2392            self.Phi =     <int *> sage_realloc(self.new_Phi,   self.Phi_size * self.L         * sizeof(int) ) 
     2393            self.Omega =   <int *> sage_realloc(self.new_Omega, self.Phi_size * self.L         * sizeof(int) ) 
     2394            self.W =       <int *> sage_realloc(self.new_W,     self.Phi_size * self.radix * 2 * sizeof(int) ) 
     2395            if not (self.w_gamma and self.alpha and self.Phi and self.Omega and self.W): 
    23342396                if self.w_gamma: sage_free(self.w_gamma) 
    23352397                if self.alpha: sage_free(self.alpha) 
    2336                 if self.new_Phi: sage_free(self.new_Phi) 
    2337                 if self.new_Omega: sage_free(self.new_Omega) 
    2338                 if self.new_W: sage_free(self.new_W) 
     2398                if self.Phi: sage_free(self.Phi) 
     2399                if self.Omega: sage_free(self.Omega) 
     2400                if self.W: sage_free(self.W) 
    23392401                raise MemoryError("Memory.") 
    23402402        word_gamma = self.w_gamma 
    … …  
    23512413        state = 1 
    23522414        while state != -1: 
    2353             print "state:", state 
    2354             print "k:", k 
    2355             print "h:", h 
    2356             print "hzf:", hzf 
    2357             print "v[k]", v[k] 
    2358             print "tvc", tvc 
    2359             print "W[k]", Integer(W[k]).binary() 
    2360             print "aut_gp_index", self.aut_gp_index 
    2361             print nu 
    2362             badass += 1 #168-195 
    2363             #if badass > 168: break 
     2415            if True:#badass > 39: 
     2416                print '-----' 
     2417                print badass 
     2418                print "k:", k 
     2419                if k != -1: 
     2420                    if v[k]&nu.flag: 
     2421                        print "v[k]: word ", v[k]^nu.flag 
     2422                    else: 
     2423                        print "v[k]: col ", v[k] 
     2424                    if tvc&nu.flag: 
     2425                        print "tvc- wd", tvc^nu.flag 
     2426                    else: 
     2427                        print "tvc- col", tvc 
     2428                    if W[self.Phi_size * k]: 
     2429                        print "W[k]: cols", Integer(W[self.Phi_size * k]).binary() 
     2430                    else: 
     2431                        j = nwords/self.radix 
     2432                        if nwords%self.radix: 
     2433                            j += 1 
     2434                        print "W[k]: words", [Integer(W[self.Phi_size * k + 1 + i]).binary() for i from 0 <= i < j] 
     2435                print nu 
     2436                if hh != -1: print zeta 
     2437            if False: 
     2438                print "h:", h 
     2439                print "hzf:", hzf__h_zeta 
     2440                print "aut_gp_index", self.aut_gp_index 
     2441                print 'hh', hh 
     2442                print 'ht', ht 
     2443                print 'hzf__h_zeta', hzf__h_zeta 
     2444                print 'qzb', qzb 
     2445            if True: 
     2446                print "state:", state 
     2447                print '-----' 
     2448            badass += 1 
     2449            #if badass > 92: break 
     2450            #if badass > 2: break 
    23642451 
    23652452            if state == 1: # Entry point: once only 
    … …  
    23722459                # store the first smallest nontrivial cell in W[k], and set v[k] 
    23732460                # equal to its minimum element 
    2374                 v[k] = nu.new_first_smallest_nontrivial(k, new_W, 2*self.radix*k) 
    2375                 W[k] = nu.first_smallest_nontrivial(k) # TODO: remove 
    2376                 v[k] = nu.v # stored during first_smallest_nontrivial # TODO: remove 
     2461                v[k] = nu.new_first_smallest_nontrivial(k, W, self.Phi_size * k) 
    23772462 
    23782463                Lambda[k] = 0 
    … …  
    23892474 
    23902475                alpha[0] = nu.split_vertex(v[k-1], k) 
    2391                 alpha[0] = nu.split_column(v[k-1], k) # TODO: remove this 
    2392  
     2476                print nu 
    23932477                Lambda[k] = nu.refine(k, alpha, 1, C, ham_wts) # store the invariant to Lambda[k] 
    23942478                # only if this is the first time moving down the search tree: 
    … …  
    24272511                # store the first smallest nontrivial cell in W[k], and set v[k] 
    24282512                # equal to its minimum element 
    2429                 v[k] = nu.new_first_smallest_nontrivial(k, new_W, 2*self.radix*k) 
    2430                 W[k] = nu.first_smallest_nontrivial(k)  # TODO: remove 
    2431                 v[k] = nu.v # stored during first_smallest_nontrivial # TODO: remove 
     2513                v[k] = nu.new_first_smallest_nontrivial(k, W, self.Phi_size * k) 
    24322514                if not nu.sat_225(k): hh = k + 1 
    24332515                e[k] = 0 # see state 12 and 17 
    … …  
    24432525            elif state == 6: # at this stage, there is no reason to continue downward, so backtrack 
    24442526                j = k 
     2527#                print 'current k', j 
     2528#                print 'ht', ht 
     2529#                print 'hzb__h_rho', hzb__h_rho 
     2530#                print 'hh', hh 
    24452531                # return to the longest ancestor nu[i] of nu that could have a 
    24462532                # descendant equivalent to zeta or could improve on rho. 
    … …  
    24732559                # pruning, so we store its fixed set and a set of representatives of its cells. 
    24742560                if l < self.L-1: l += 1 
    2475                 Omega[l] = nu.min_cell_reps(hh) 
    2476                 Phi[l] = nu.fixed_cols(Omega[l], hh) 
     2561                nu.new_min_cell_reps(hh, Omega, self.Phi_size*l) 
     2562                nu.fixed_vertices(hh, Phi, Omega, self.Phi_size*l) 
     2563 
    24772564                state = 12 
    24782565 
    … …  
    24842571                # hzf is the extremal height of ancestors of both nu and zeta, so if k < hzf, nu is not 
    24852572                # equivalent to zeta, i.e. there is no automorphism to discover. 
    2486                 if k < hzf: state = 8; continue 
     2573                if k < hzf__h_zeta: state = 8; continue 
    24872574 
    24882575                nu.get_permutation(zeta, word_gamma, col_gamma, ham_wts) 
    … …  
    25112598 
    25122599                # if C(nu) == C(rho), get the automorphism and goto 10 
    2513                 nu.find_basis(ham_wts) 
    25142600                rho.get_permutation(nu, word_gamma, col_gamma, ham_wts) 
    25152601                print "gamma:", [word_gamma[i] for i from 0 <= i < nwords], [col_gamma[i] for i from 0 <= i < ncols] 
    … …  
    25292615                # increment l 
    25302616                if l < self.L-1: l += 1 
     2617 
    25312618                # store information about the automorphism to Omega and Phi 
     2619                ii = self.Phi_size*l 
     2620                Omega[ii] = ~(~0 << ncols) 
     2621                Phi[ii] = 0 
     2622                jj = 1 + nwords/self.radix 
     2623                if nwords%self.radix: 
     2624                    jj += 1 
     2625                for i from 0 < i < jj: 
     2626                    Omega[ii+i] = ~0 
     2627                    Phi[ii+i] = 0 
     2628                Omega[ii+jj-1] = ~(~0 << nwords%self.radix) 
    25322629                # Omega stores the minimum cell representatives 
    2533                 Omega[l] = ~0 
    25342630                i = 0 
    25352631                while i < ncols: 
    25362632                    j = col_gamma[i]         # i is a minimum 
    25372633                    while j != i:            # cell rep, 
    2538                         Omega[l] ^= (1<<j)   # so cancel 
     2634                        Omega[ii] ^= (1<<j)  # so cancel 
    25392635                        j = col_gamma[j]     # cellmates 
    25402636                    i += 1 
    2541                     while i < ncols and not Omega[l]&(1<<i): # find minimal element 
    2542                         i += 1                               # of next cell 
     2637                    while i < ncols and not Omega[ii]&(1<<i): # find minimal element 
     2638                        i += 1                                # of next cell 
     2639                i = 0 
     2640                jj = self.radix 
     2641                while i < nwords: 
     2642                    j = word_gamma[i] 
     2643                    while j != i: 
     2644                        Omega[ii+1+j/jj] ^= (1<<(j%jj)) 
     2645                        j = word_gamma[j] 
     2646                    i += 1 
     2647                    while i < nwords and not Omega[ii+1+i/jj]&(1<<(i%jj)): 
     2648                        i += 1 
    25432649                # Phi stores the columns fixed by the automorphism 
    2544                 Phi[l] = 0 
    25452650                for i from 0 <= i < ncols: 
    25462651                    if col_gamma[i] == i: 
    2547                         Phi[l] ^= (1 << i) 
     2652                        Phi[ii] ^= (1 << i) 
     2653                for i from 0 <= i < nwords: 
     2654                    if word_gamma[i] == i: 
     2655                        Phi[ii+1+i/jj] ^= (1<<(i%jj)) 
     2656 
    25482657                # Now incorporate the automorphism into Theta 
    25492658                j = Theta.merge_perm(col_gamma, word_gamma) 
    2550                 if not j: print "no j" 
     2659#                if not j: print "no j" 
    25512660                # j stores whether anything happened or not- if not, then the automorphism we have 
    25522661                # discovered is already in the subgroup spanned by the generators we have output 
    … …  
    25582667                # algorithm came up to meet zeta. At this point, we were considering the new 
    25592668                # possibilities for descending away from zeta at this level. 
    2560                 if Theta.col_min_cell_rep[Theta.col_find(tvc)] == tvc: 
    2561                     # if this is still a minimum cell representative of Theta, even in light 
    2562                     # of this new automorphism, then the current branch off of zeta hasn't been 
    2563                     # found equivalent to one already searched yet, so there may still be a 
    2564                     # better canonical label downward. 
    2565                     state = 11; continue 
     2669                # if this is still a minimum cell representative of Theta, even in light 
     2670                # of this new automorphism, then the current branch off of zeta hasn't been 
     2671                # found equivalent to one already searched yet, so there may still be a 
     2672                # better canonical label downward. 
     2673                if tvc & nu.flag: 
     2674                    i = tvc^nu.flag 
     2675                    if Theta.wd_min_cell_rep[Theta.wd_find(i)] == i: 
     2676                        state = 11; continue 
     2677                else: 
     2678                    if Theta.col_min_cell_rep[Theta.col_find(tvc)] == tvc: 
     2679                        state = 11; continue 
    25662680 
    25672681                # Otherwise, proceed to where zeta meets nu: 
    … …  
    25892703                    # before, so we have already intersected W[k] with the bulk of Omega and Phi, but 
    25902704                    # we should still catch up with the latest ones 
    2591                     W[k] &= Omega[l] 
     2705                    ii = self.Phi_size*l 
     2706                    jj = self.Phi_size*k 
     2707                    j = 1 + nwords/self.radix 
     2708                    if nwords%self.radix: 
     2709                        j += 1 
     2710                    W[jj] &= Omega[ii] 
     2711                    for i from 0 < i < j: 
     2712                        W[jj+i] &= Omega[ii+i] 
    25922713                state = 13 
    25932714 
    … …  
    26022723                h = k 
    26032724                tvc = 0 
    2604                 while not (1 << tvc) & W[k]: 
    2605                     tvc += 1 
     2725                jj = self.Phi_size*k 
     2726                if W[jj]: 
     2727#                    print 'W[jj]', W[jj] 
     2728#                    print tvc 
     2729                    while not (1 << tvc) & W[jj]: 
     2730                        tvc += 1 
     2731                else: 
     2732                    ii = 0 
     2733                    while not W[jj+1+ii]: 
     2734                        ii += 1 
     2735                    while not W[jj+1+ii] & (1 << tvc): 
     2736                        tvc += 1 
     2737                    tvc = (ii*self.radix + tvc) ^ nu.flag 
    26062738                # now tvc points to the minimal cell representative of W[k] 
    26072739                state = 14 
    26082740 
    26092741            elif state == 14: # see if there are any more splits to make from this level of zeta (see state 17) 
    2610                 if Theta.col_find(v[k]) == Theta.col_find(tvc): 
    2611                     index += 1 
    2612                     # keep tabs on how many elements are in the same cell of Theta as tvc 
     2742                if v[k]&nu.flag == tvc&nu.flag: 
     2743                    if tvc&nu.flag: 
     2744                        if Theta.wd_find(v[k]^nu.flag) == Theta.wd_find(tvc^nu.flag): 
     2745                            index += 1 
     2746                    else: 
     2747                        if Theta.col_find(v[k]) == Theta.col_find(tvc): 
     2748                            index += 1 
     2749                            # keep tabs on how many elements are in the same cell of Theta as tvc 
    26132750                # find the next split 
    2614                 i = v[k] + 1 
    2615                 while i < ncols and not (1 << i) & W[k]: 
    2616                     i += 1 
    2617                 if i < ncols: 
    2618                     v[k] = i 
     2751                jj = self.Phi_size*k 
     2752                if v[k]&nu.flag: 
     2753                    ii = self.radix 
     2754                    i = (v[k]^nu.flag) + 1 
     2755                    while i < nwords and not (1 << i%ii) & W[jj+1+i/ii]: 
     2756                        i += 1 
     2757                    if i < nwords: 
     2758                        v[k] = i^nu.flag 
     2759                    else: 
     2760                        # there is no new split at this level 
     2761                        state = 16; continue 
     2762                    # new split column better be a minimal representative in Theta, or wasted effort 
     2763                    if Theta.wd_min_cell_rep[Theta.wd_find(i)] == i: 
     2764                        state = 15 
     2765                    else: 
     2766                        state = 14 
    26192767                else: 
    2620                     # there is no new split at this level 
    2621                     v[k] = -1 
    2622                     state = 16; continue 
    2623  
    2624                 # new split column better be a minimal representative in Theta, or wasted effort 
    2625                 if Theta.col_min_cell_rep[Theta.col_find(tvc)] == tvc: 
    2626                     state = 15 
    2627                 else: 
    2628                     state = 14 
     2768                    i = v[k] + 1 
     2769                    while i < ncols and not (1 << i) & W[jj]: 
     2770                        i += 1 
     2771                    if i < ncols: 
     2772                        v[k] = i 
     2773                    else: 
     2774                        # there is no new split at this level 
     2775                        state = 16; continue 
     2776                    # new split column better be a minimal representative in Theta, or wasted effort 
     2777                    if Theta.col_min_cell_rep[Theta.col_find(v[k])] == v[k]: 
     2778                        state = 15 
     2779                    else: 
     2780                        state = 14 
    26292781 
    26302782            elif state == 15: # split out the column v[k] 
    … …  
    26342786                    hh = k + 1 
    26352787                # hzf is maximal such that indicators line up for nu and zeta 
    2636                 if k < hzf: 
    2637                     hzf = k 
     2788                if k < hzf__h_zeta: 
     2789                    hzf__h_zeta = k 
    26382790                # hb is longest common ancestor of nu and rho 
    26392791                if hb >= k: 
    … …  
    26432795 
    26442796            elif state == 16: # backtrack up zeta, updating information about stabilizer vector 
    2645                 i = W[k] 
    2646                 j = ham_wts[i & 65535] + ham_wts[(i >> 16) & 65535] 
     2797                jj = self.Phi_size*k 
     2798                if W[jj]: 
     2799                    i = W[jj] 
     2800                    j = ham_wts[i & 65535] + ham_wts[(i >> 16) & 65535] 
     2801                else: 
     2802                    i = 0; j = 0 
     2803                    ii = self.radix 
     2804                    while i*ii < nwords: 
     2805                        iii = W[jj+1+i] 
     2806                        j += ham_wts[iii & 65535] + ham_wts[(iii >> 16) & 65535] 
     2807                        i += 1 
    26472808                if j == index and ht == k + 1: ht = k 
    26482809                size = size*index 
    … …  
    26552816                if e[k] == 0: # now is the time to narrow down W[k] by Omega and Phi 
    26562817                    # intersect W[k] with each Omega[i] such that v[0]...v[k-1] is in Phi[i] 
    2657                     ii = 0 
    2658                     for i from 0 <= i < k: 
    2659                         ii ^= (1 << v[i]) 
     2818                    jjj = self.Phi_size*k 
     2819                    jj = self.Phi_size*self.L 
     2820                    iii = nwords/self.radix + 1 
     2821                    if nwords%self.radix: 
     2822                        iii += 1 
     2823                    for ii from 0 <= ii < iii: 
     2824                        Phi[jj+ii] = 0 
     2825                    for ii from 0 <= ii < k: 
     2826                        if v[ii]&nu.flag: 
     2827                            i = v[ii]^nu.flag 
     2828                            Phi[jj+1+i/self.radix] ^= (1 << i%self.radix) 
     2829                        else: 
     2830                            Phi[jj] ^= (1 << v[ii]) 
    26602831                    for i from 0 <= i <= l: 
    2661                         if Phi[i] & ii == ii: 
    2662                             W[k] &= Omega[i] 
     2832                        ii = self.Phi_size*i 
     2833                        for j from 0 <= j < iii: 
     2834                            if Phi[ii + j] & Phi[jj + j] == Phi[jj + j]: 
     2835                                W[jjj + j] &= Omega[ii + j] 
    26632836                e[k] = 1 
    26642837                # see if there is a vertex to split out 
    2665                 i = v[k] + 1 
    2666                 while i < ncols: 
    2667                     i += 1 
    2668                     if (1 << i) & W[k]: break 
    2669                 if i < ncols: 
    2670                     v[k] = i 
    2671                     state = 15; continue 
    2672  
    2673                 v[k] = -1 
     2838                if nu.flag&v[k]: 
     2839                    i = (v[k]^nu.flag) + 1 
     2840                    while i < nwords: 
     2841                        i += 1 
     2842                        if (1 << i%self.radix) & W[jjj+1+i/self.radix]: break 
     2843                    if i < nwords: 
     2844                        v[k] = i^nu.flag 
     2845                        state = 15; continue 
     2846                else: 
     2847                    i = v[k] + 1 
     2848                    while i < ncols: 
     2849                        i += 1 
     2850                        if (1 << i) & W[jjj]: break 
     2851                    if i < ncols: 
     2852                        v[k] = i 
     2853                        state = 15; continue 
     2854 
    26742855                k -= 1 
    26752856                state = 13 
    … …  
    26792860                h = k # zeta[h] == nu[h] 
    26802861                ht = k # nodes descended from zeta[ht] are all equivalent 
    2681                 hzf = k # max such that indicators for zeta and nu agree 
     2862                hzf__h_zeta = k # max such that indicators for zeta and nu agree 
    26822863                zeta = PartitionStack(nu) 
    26832864                # (POINT B) 
    … …  
    26922873 
    26932874        # end big while loop 
    2694         rho.find_basis(ham_wts) 
    2695         for i from 0 <= i < ncols: 
    2696             self.labeling[rho.col_ents[i]] = i 
    2697         for i from 0 <= i < nrows: 
    2698             self.labeling[i+ncols] = rho.basis_locations[i] 
    2699  
    2700  
    2701  
    2702  
    2703  
     2875#        rho.find_basis(ham_wts) 
     2876#        for i from 0 <= i < ncols: 
     2877#            self.labeling[rho.col_ents[i]] = i 
     2878#        for i from 0 <= i < nrows: 
     2879#            self.labeling[i+ncols] = rho.basis_locations[i] 
     2880 
     2881 
     2882 
     2883 
     2884 
Note: See TracChangeset for help on using the changeset viewer.