Ignore:
Files:
3 deleted
8 edited

Legend:

Unmodified
Added
Removed

  • sage/matrix/matrix_integer_dense.pxd

    r2653 r1769  
    11include "../ext/cdefs.pxi" 
    2  
    3 cdef extern from "../ext/multi_modular.h": 
    4     ctypedef unsigned long mod_int 
    52 
    63cimport matrix_dense 
    … …  
    129    cdef object _pivots 
    1310    cdef int mpz_height(self, mpz_t height) except -1 
    14     cdef _mod_int_c(self, mod_int modulus) 
    1511 
    1612    cdef void _zero_out_matrix(self) 

  • sage/matrix/matrix_integer_dense.pyx

    r2654 r2629  
    2020ctypedef unsigned int uint 
    2121 
    22 from sage.ext.multi_modular cimport MultiModularBasis 
    23 cdef MultiModularBasis mm 
    24 mm = MultiModularBasis() 
    25  
    2622from sage.rings.integer cimport Integer 
    2723from sage.rings.rational_field import QQ 
    2824from sage.rings.integer_ring import ZZ 
    29 from sage.rings.integer_mod_ring import IntegerModRing 
    3025from sage.rings.polynomial_ring import PolynomialRing 
    3126from sage.structure.element cimport ModuleElement 
    … …  
    4035import matrix_space 
    4136 
     37# for multi-modular methods 
     38# TODO: make this a const, adjust for machine word size, put somewhere more reasonable 
     39start_prime = 32771 
    4240 
    4341cdef class Matrix_integer_dense(matrix_dense.Matrix_dense):   # dense or sparse 
    … …  
    607605         
    608606    def _multiply_multi_modular(left, Matrix_integer_dense right): 
    609      
    610         cdef Integer h 
    611         cdef mod_int *moduli 
    612         cdef int i, n 
    613          
    614607        h = left.height() * right.height() 
    615         n = mm.moduli_list_c(&moduli, h.value) 
    616         res = [] 
    617         for i from 0 <= i < n: 
    618             res.append(left._mod_int_c(moduli[i]) * right._mod_int_c(moduli[i])) 
    619         sage_free(moduli) 
    620         return left._lift_crt(res) 
     608        product = start_prime 
     609#        moduli = [ start_prime ] 
     610        residues = [ left._mod_int(start_prime) * right._mod_int(start_prime) ] 
     611        cur_prime = start_prime 
     612        while product < height: 
     613            cur_prime = next_prime(cur_prime) 
     614            product *= cur_prime 
     615#            moduli.append(cur_prime) 
     616            res.append(left._mod_int(cur_prime) * right_res._mod_int(cur_prime)) 
     617        return left._lift_crt(residues) 
    621618             
    622619    def _mod_int(self, modulus): 
    623         return self._mod_int_c(modulus) 
    624  
    625     cdef _mod_int_c(self, mod_int p): 
     620        cdef int p 
    626621        cdef Py_ssize_t i, j 
     622        p = modulus 
    627623        cdef Matrix_modn_dense res 
    628624        cdef mpz_t* self_row 
    629         cdef mod_int* res_row 
    630         res = Matrix_modn_dense.__new__(Matrix_modn_dense, matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None) 
     625        cdef uint* res_row 
     626        res = Matrix_modn_dense.__new__(matrix_space.MatrixSpace(ZZ, self._nrows, self._ncols, sparse=False), None, None, None) 
    631627        for i from 0 <= i < self._nrows: 
    632628            self_row = self._matrix[i] 
    … …  
    636632        return res 
    637633         
    638     def _lift_crt(self, residues): 
    639  
    640         cdef size_t n, i, j, k 
    641         cdef Py_ssize_t nr, nc 
    642          
    643         n = len(residues) 
    644         nr = residues[0].nrows() 
    645         nc = residues[0].ncols() 
    646          
    647         for b in residues: 
    648             if not PY_TYPE_CHECK(b, Matrix_modn_dense): 
    649                 raise TypeError, "Can only perform CRT on list of type Matrix_modn_dense." 
    650         cdef PyObject** res 
    651         res = FAST_SEQ_UNSAFE(residues) 
    652  
    653         cdef mod_int **row_list 
    654         row_list = <mod_int**>sage_malloc(sizeof(mod_int*) * n) 
    655         if row_list == NULL: 
    656             raise MemoryError, "out of memory allocating multi-modular coefficent list" 
    657          
    658         cdef Matrix_integer_dense M 
    659         M = Matrix_integer_dense.__new__(Matrix_integer_dense, self.matrix_space(nr, nc), None, None, None) 
    660          
    661         _sig_on 
    662         for i from 0 <= i < nr: 
    663             for k from 0 <= k < n: 
    664                 row_list[k] = (<Matrix_modn_dense>res[k]).matrix[i] 
    665             mm.mpz_crt_vec(M._matrix[i], row_list, n, nc) 
    666         _sig_off 
    667          
    668         sage_free(row_list) 
    669         return M 
    670          
     634    def _lift_crt(residues): 
     635        raise NotImplementedError 
    671636 
    672637    def _echelon_in_place_classical(self): 

  • sage/matrix/matrix_modn_dense.pxd

    r2847 r1772  
    11cimport matrix_dense 
    22 
    3 cdef extern from "../ext/multi_modular.h": 
    4     ctypedef unsigned long mod_int 
    5     mod_int MOD_INT_MAX 
    6     mod_int MOD_INT_OVERFLOW 
     3ctypedef unsigned int uint 
    74 
    85cdef class Matrix_modn_dense(matrix_dense.Matrix_dense): 
    9     cdef mod_int **matrix 
    10     cdef mod_int *_entries 
    11     cdef mod_int p 
    12     cdef mod_int gather 
     6    cdef uint **matrix 
     7    cdef uint *_entries 
     8    cdef uint p 
     9    cdef uint gather 
    1310 
    14     #cdef set_matrix(Matrix_modn_dense self, mod_int **m) 
    15     #cdef mod_int **get_matrix(Matrix_modn_dense self) 
    16     #cdef mod_int entry(self, mod_int i, mod_int j) 
    17     cdef _rescale_row_c(self, Py_ssize_t row, mod_int multiple, Py_ssize_t start_col) 
    18     cdef _rescale_col_c(self, Py_ssize_t col, mod_int multiple, Py_ssize_t start_row)     
     11    #cdef set_matrix(Matrix_modn_dense self, uint **m) 
     12    #cdef uint **get_matrix(Matrix_modn_dense self) 
     13    #cdef uint entry(self, uint i, uint j) 
     14    cdef _rescale_row_c(self, Py_ssize_t row, uint multiple, Py_ssize_t start_col) 
     15    cdef _rescale_col_c(self, Py_ssize_t col, uint multiple, Py_ssize_t start_row)     
    1916    cdef _add_multiple_of_row_c(self,  Py_ssize_t row_to, Py_ssize_t row_from, 
    20                                 mod_int multiple, Py_ssize_t start_col) 
     17                                uint multiple, Py_ssize_t start_col) 
    2118    cdef _add_multiple_of_column_c(self, Py_ssize_t col_to, Py_ssize_t col_from, 
    22                                    mod_int multiple, Py_ssize_t start_row) 
     19                                   uint multiple, Py_ssize_t start_row) 
    2320     
    2421 

  • sage/matrix/matrix_modn_dense.pyx

    r2847 r2629  
    112112        matrix_dense.Matrix_dense.__init__(self, parent) 
    113113 
    114         cdef mod_int p 
     114        cdef uint p 
    115115        p = self._base_ring.characteristic() 
    116116        self.p = p 
    117         if p >= MOD_INT_MAX: 
    118             raise OverflowError, "p (=%s) must be < %s"%(p, MOD_INT_MAX) 
    119         self.gather = MOD_INT_OVERFLOW / (p*p) 
     117        if p >= 46340: 
     118            raise OverflowError, "p (=%s) must be < 46340"%p 
     119        self.gather = 2**32/(p*p) 
    120120             
    121         self._entries = <mod_int *> sage_malloc(sizeof(mod_int)*self._nrows*self._ncols) 
     121        self._entries = <uint *> sage_malloc(sizeof(uint)*self._nrows*self._ncols) 
    122122        if self._entries == NULL: 
    123123           raise MemoryError, "Error allocating matrix"         
    124124 
    125         self.matrix = <mod_int **> sage_malloc(sizeof(mod_int*)*self._nrows) 
     125        self.matrix = <uint **> sage_malloc(sizeof(uint*)*self._nrows) 
    126126        if self.matrix == NULL: 
    127127            sage_free(self._entries) 
    128128            raise MemoryError, "Error allocating memory" 
    129129 
    130         cdef mod_int k 
     130        cdef uint k 
    131131        cdef Py_ssize_t i 
    132132        k = 0 
    … …  
    143143    def __init__(self, parent, entries, copy, coerce): 
    144144         
    145         cdef mod_int e 
     145        cdef uint e 
    146146        cdef Py_ssize_t i, j, k 
    147         cdef mod_int *v 
     147        cdef uint *v 
    148148 
    149149        # scalar? 
    … …  
    164164         
    165165        k = 0 
    166         cdef mod_int n 
     166        cdef uint n 
    167167        R = self.base_ring() 
    168168 
    … …  
    229229        A = Matrix_modn_dense.__new__(Matrix_modn_dense, self._parent, 
    230230                                      0, 0, 0) 
    231         memcpy(A._entries, self._entries, sizeof(mod_int)*self._nrows*self._ncols) 
     231        memcpy(A._entries, self._entries, sizeof(uint)*self._nrows*self._ncols) 
    232232        A.p = self.p 
    233233        A.gather = self.gather 
    … …  
    257257 
    258258        cdef Py_ssize_t start_row, c, r, nr, nc, i 
    259         cdef mod_int p, a, a_inverse, b 
    260         cdef mod_int **m 
     259        cdef uint p, a, a_inverse, b 
     260        cdef uint **m 
    261261 
    262262        start_row = 0 
    … …  
    295295        self._rescale_row_c(row, multiple, start_col) 
    296296         
    297     cdef _rescale_row_c(self, Py_ssize_t row, mod_int multiple, Py_ssize_t start_col): 
    298         cdef mod_int r, p 
    299         cdef mod_int* v 
     297    cdef _rescale_row_c(self, Py_ssize_t row, uint multiple, Py_ssize_t start_col): 
     298        cdef uint r, p 
     299        cdef uint* v 
    300300        cdef Py_ssize_t i 
    301301        p = self.p 
    … …  
    307307        self._rescale_col_c(col, multiple, start_row) 
    308308         
    309     cdef _rescale_col_c(self, Py_ssize_t col, mod_int multiple, Py_ssize_t start_row): 
     309    cdef _rescale_col_c(self, Py_ssize_t col, uint multiple, Py_ssize_t start_row): 
    310310        """ 
    311311        EXAMPLES: 
    … …  
    339339            [ 2  5 34]             
    340340        """ 
    341         cdef mod_int r, p 
    342         cdef mod_int* v 
     341        cdef uint r, p 
     342        cdef uint* v 
    343343        cdef Py_ssize_t i 
    344344        p = self.p 
    … …  
    350350        self._add_multiple_of_row_c(row_to, row_from, multiple, start_col) 
    351351         
    352     cdef _add_multiple_of_row_c(self,  Py_ssize_t row_to, Py_ssize_t row_from, mod_int multiple, 
     352    cdef _add_multiple_of_row_c(self,  Py_ssize_t row_to, Py_ssize_t row_from, uint multiple, 
    353353                               Py_ssize_t start_col): 
    354         cdef mod_int p 
    355         cdef mod_int *v_from, *v_to 
     354        cdef uint p 
     355        cdef uint *v_from, *v_to 
    356356 
    357357        p = self.p 
    … …  
    368368        self._add_multiple_of_column_c(col_to, col_from, s, start_row) 
    369369 
    370     cdef _add_multiple_of_column_c(self, Py_ssize_t col_to, Py_ssize_t col_from, mod_int multiple, 
     370    cdef _add_multiple_of_column_c(self, Py_ssize_t col_to, Py_ssize_t col_from, uint multiple, 
    371371                                   Py_ssize_t start_row): 
    372         cdef mod_int  p 
    373         cdef mod_int **m 
     372        cdef uint  p 
     373        cdef uint **m 
    374374         
    375375        m = self.matrix 
    … …  
    382382 
    383383    cdef swap_rows_c(self, Py_ssize_t row1, Py_ssize_t row2): 
    384         cdef mod_int* temp 
     384        cdef uint* temp 
    385385        temp = self.matrix[row1] 
    386386        self.matrix[row1] = self.matrix[row2] 
    … …  
    389389    cdef swap_columns_c(self, Py_ssize_t col1, Py_ssize_t col2): 
    390390        cdef Py_ssize_t i, nr 
    391         cdef mod_int t 
    392         cdef mod_int **m 
     391        cdef uint t 
     392        cdef uint **m 
    393393        m = self.matrix 
    394394        nr = self._nrows 
    … …  
    412412        n = self._nrows 
    413413 
    414         cdef mod_int **h 
     414        cdef uint **h 
    415415        h = self.matrix 
    416416 
    417         cdef mod_int p, r, t, t_inv, u 
     417        cdef uint p, r, t, t_inv, u 
    418418        cdef Py_ssize_t i, j, m 
    419419        p = self.p 
    … …  
    477477        n = self._nrows 
    478478 
    479         cdef mod_int p, t 
     479        cdef uint p, t 
    480480        p = self.p 
    481481 
    … …  
    529529    #  A = matrix(Integers(389),4,range(16)); A._echelon_strassen(4) 
    530530    # *** glibc detected *** free(): invalid next size (fast): 0x0000000000fb15e0 *** 
    531     # error in set_to memcpy on 64-bit 
    532     def matrix_window(self, Py_ssize_t row=0, Py_ssize_t col=0, 
     531    def xxx_matrix_window(self, Py_ssize_t row=0, Py_ssize_t col=0, 
    533532                      Py_ssize_t nrows=-1, Py_ssize_t ncols=-1): 
    534533        """ 
    … …  
    551550            <type 'sage.matrix.matrix_modn_dense.Matrix_modn_dense'> 
    552551        """ 
    553         print "making matrix window" 
    554552        if nrows == -1: 
    555553            nrows = self._nrows - row 

  • sage/matrix/matrix_rational_dense.pyx

    r2654 r2628  
    439439        A, A_denom = left._clear_denom() 
    440440        B, B_denom = right._clear_denom() 
    441         AB = A._multiply_multi_modular(B) 
     441        AB = A*B # A._multiply_multi_modular(B) 
    442442        D = A_denom * B_denom 
    443443        res = Matrix_rational_dense.__new__(Matrix_rational_dense, left.matrix_space(AB._nrows, AB._ncols), 0, 0, 0) 

  • sage/matrix/matrix_window_modn_dense.pxd

    r2847 r1872  
    11from matrix_window cimport MatrixWindow 
    22 
    3 cdef extern from "../ext/multi_modular.h": 
    4     ctypedef unsigned long mod_int 
    5     mod_int MOD_INT_MAX 
     3ctypedef unsigned int uint 
    64 
    75cdef class MatrixWindow_modn_dense(MatrixWindow): 

  • sage/matrix/matrix_window_modn_dense.pyx

    r2847 r1993  
    88from sage.matrix.matrix_modn_dense cimport Matrix_modn_dense 
    99 
    10 # TODO: what if we used the actual pointers for our loop variables? 
    1110 
    1211cdef class MatrixWindow_modn_dense(matrix_window.MatrixWindow): 
    … …  
    2423        """ 
    2524        cdef Py_ssize_t i, j 
    26         cdef mod_int** self_rows 
    27         cdef mod_int** A_rows 
     25        cdef uint** self_rows 
     26        cdef uint** A_rows 
    2827        if self._nrows != A._nrows or self._ncols != A._ncols: 
    2928            raise ArithmeticError, "incompatible dimensions" 
    … …  
    3130        A_rows = ( <Matrix_modn_dense> A._matrix ).matrix 
    3231        for i from 0 <= i < self._nrows: 
    33             memcpy(self_rows[i+self._row] + self._col, A_rows[i+A._row] + A._col, self._ncols * sizeof(mod_int*)) 
     32            memcpy(self_rows[i+self._row] + self._col, A_rows[i+A._row] + A._col, self._ncols * sizeof(uint*)) 
    3433                 
    3534    cdef set_to_zero(MatrixWindow_modn_dense self): 
    3635        cdef Py_ssize_t i, j 
    37         cdef mod_int** rows 
     36        cdef uint** rows 
    3837        rows = ( <Matrix_modn_dense> self._matrix ).matrix 
    3938        for i from self._row <= i < self._row + self._nrows: 
    40             memset(rows[i] + self._col, 0, self._ncols * sizeof(mod_int*)) 
     39            memset(rows[i] + self._col, 0, self._ncols * sizeof(uint*)) 
    4140 
    4241    cdef add(self, MatrixWindow A): 
    4342        cdef Py_ssize_t i, j 
    44         cdef mod_int p 
    45         cdef mod_int* self_row 
    46         cdef mod_int* A_row 
     43        cdef uint p 
     44        cdef uint** self_matrix 
     45        cdef uint** A_matrix 
    4746        if self._nrows != A._nrows or self._ncols != A._ncols: 
    4847            raise ArithmeticError, "incompatible dimensions" 
     48        self_matrix = ( <Matrix_modn_dense> self._matrix ).matrix 
     49        A_matrix = ( <Matrix_modn_dense> A._matrix ).matrix 
    4950        p = ( <Matrix_modn_dense> self._matrix ).p 
    5051        for i from 0 <= i < self._nrows: 
    51             self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col 
    52             A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col 
    5352            for j from 0 <= j < self._ncols: 
    54                 self_row[j] += A_row[j] 
    55                 if self_row[j] >= p: 
    56                     self_row[j] -= p 
     53                # I really want to do in place operations here... 
     54                self_matrix[i+self._row][j+self._col] = self_matrix[i+self._row][j+self._col] + A_matrix[i+A._row][j+A._col] 
     55                if self_matrix[i+self._row][j+self._col] >= p: 
     56                    self_matrix[i+self._row][j+self._col] = self_matrix[i+self._row][j+self._col] - p 
    5757     
    5858    cdef subtract(self, MatrixWindow A): 
    5959        cdef Py_ssize_t i, j 
    60         cdef mod_int p 
    61         cdef mod_int* self_row 
    62         cdef mod_int* A_row 
     60        cdef uint p 
     61        cdef uint** self_matrix 
     62        cdef uint** A_matrix 
    6363        if self._nrows != A._nrows or self._ncols != A._ncols: 
    6464            raise ArithmeticError, "incompatible dimensions" 
     65        self_matrix = ( <Matrix_modn_dense> self._matrix ).matrix 
     66        A_matrix = ( <Matrix_modn_dense> A._matrix ).matrix 
    6567        p = ( <Matrix_modn_dense> self._matrix ).p 
    6668        for i from 0 <= i < self._nrows: 
    67             self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col 
    68             A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col 
    6969            for j from 0 <= j < self._ncols: 
    70                 if self_row[j] >= A_row[j]: 
    71                     self_row[j] -= A_row[j] 
    72                 else: 
    73                     self_row[j] += p - A_row[j] 
     70                # I really want to do in place operations here... 
     71                self_matrix[i+self._row][j+self._col] = p + self_matrix[i+self._row][j+self._col] - A_matrix[i+A._row][j+A._col] 
     72                if self_matrix[i+self._row][j+self._col] >= p: 
     73                    self_matrix[i+self._row][j+self._col] = self_matrix[i+self._row][j+self._col] - p 
    7474         
    7575    cdef set_to_sum(self, MatrixWindow A, MatrixWindow B): 
    7676        cdef Py_ssize_t i, j 
    77         cdef mod_int p 
    78         cdef mod_int* self_row 
    79         cdef mod_int* A_row 
    80         cdef mod_int* B_row 
     77        cdef uint p 
     78        cdef uint** self_matrix 
     79        cdef uint** A_matrix 
     80        cdef uint** B_matrix 
    8181        if self._nrows != A._nrows or self._ncols != A._ncols: 
    8282            raise ArithmeticError, "incompatible dimensions" 
    8383        if self._nrows != B._nrows or self._ncols != B._ncols: 
    8484            raise ArithmeticError, "incompatible dimensions" 
     85        self_matrix = ( <Matrix_modn_dense> self._matrix ).matrix 
     86        A_matrix = ( <Matrix_modn_dense> A._matrix ).matrix 
     87        B_matrix = ( <Matrix_modn_dense> B._matrix ).matrix 
    8588        p = ( <Matrix_modn_dense> self._matrix ).p 
    8689        for i from 0 <= i < self._nrows: 
    87             self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col 
    88             A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col 
    89             B_row    = ( <Matrix_modn_dense>    B._matrix ).matrix[i +    B._row] + B._col 
    9090            for j from 0 <= j < self._ncols: 
    91                 self_row[j] = A_row[j] + B_row[j] 
    92                 if self_row[j] >= p: 
    93                     self_row[j] -= p 
     91                self_matrix[i+self._row][j+self._col] = A_matrix[i+A._row][j+A._col] + B_matrix[i+B._row][j+B._col] 
     92                if self_matrix[i+self._row][j+self._col] >= p: 
     93                    # I really want to do in place operations here... 
     94                    self_matrix[i+self._row][j+self._col] = self_matrix[i+self._row][j+self._col] - p 
    9495 
    9596    cdef set_to_diff(self, MatrixWindow A, MatrixWindow B): 
    9697        cdef Py_ssize_t i, j 
    97         cdef mod_int p 
    98         cdef mod_int* self_row 
    99         cdef mod_int* A_row 
    100         cdef mod_int* B_row 
     98        cdef uint p 
     99        cdef uint** self_matrix 
     100        cdef uint** A_matrix 
     101        cdef uint** B_matrix 
    101102        if self._nrows != A._nrows or self._ncols != A._ncols: 
    102103            raise ArithmeticError, "incompatible dimensions" 
    103104        if self._nrows != B._nrows or self._ncols != B._ncols: 
    104105            raise ArithmeticError, "incompatible dimensions" 
     106        self_matrix = ( <Matrix_modn_dense> self._matrix ).matrix 
     107        A_matrix = ( <Matrix_modn_dense> A._matrix ).matrix 
     108        B_matrix = ( <Matrix_modn_dense> B._matrix ).matrix 
    105109        p = ( <Matrix_modn_dense> self._matrix ).p 
    106110        for i from 0 <= i < self._nrows: 
    107             self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col 
    108             A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col 
    109             B_row    = ( <Matrix_modn_dense>    B._matrix ).matrix[i +    B._row] + B._col 
    110111            for j from 0 <= j < self._ncols: 
    111                 if A_row[j] > B_row[j]: 
    112                     self_row[j] = A_row[j] - B_row[j] 
    113                 else: 
    114                     self_row[j] = p + A_row[j] - B_row[j] 
     112                self_matrix[i+self._row][j+self._col] =  p + A_matrix[i+A._row][j+A._col] - B_matrix[i+B._row][j+B._col] 
     113                if self_matrix[i+self._row][j+self._col] >= p: 
     114                    # I really want to do in place operations here... 
     115                    self_matrix[i+self._row][j+self._col] = self_matrix[i+self._row][j+self._col] - p 
    115116                 
    116117    cdef set_to_prod(self, MatrixWindow A, MatrixWindow B): 
    117         cdef Py_ssize_t i, j, k, gather, top, A_ncols 
    118         cdef mod_int p, s 
    119         cdef mod_int* self_row 
    120         cdef mod_int* A_row 
    121         cdef mod_int** B_matrix_off 
     118        # TODO: gather 
     119        cdef Py_ssize_t i, j, k 
     120        cdef uint p, s 
     121        cdef uint** self_matrix 
     122        cdef uint** A_matrix 
     123        cdef uint** B_matrix 
    122124        if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols: 
    123125            raise ArithmeticError, "incompatible dimensions" 
    124         B_matrix_off = ( <Matrix_modn_dense> B._matrix ).matrix + B._row 
     126        self_matrix = ( <Matrix_modn_dense> self._matrix ).matrix 
     127        A_matrix = ( <Matrix_modn_dense> A._matrix ).matrix 
     128        B_matrix = ( <Matrix_modn_dense> B._matrix ).matrix 
    125129        p = ( <Matrix_modn_dense> self._matrix ).p 
    126         gather = ( <Matrix_modn_dense> self._matrix ).gather 
    127         A_ncols = A._ncols 
    128          
    129         if gather <= 1: 
    130             for i from 0 <= i < A._nrows: 
    131                 self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col 
    132                 A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col 
    133                 for j from 0 <= j < B._ncols: 
    134                     s = 0 
    135                     for k from 0 <= k < A._ncols: 
    136                         s = (s + A_row[k] * B_matrix_off[k][j+B._col]) % p 
    137                     self_row[j] = s 
    138         else: 
    139             for i from 0 <= i < A._nrows: 
    140                 self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col 
    141                 A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col 
    142                 for j from 0 <= j < B._ncols: 
    143                     s = 0 
    144                     k = 0 
    145                     while k < A_ncols: 
    146                         top = k + gather 
    147                         if top > A_ncols: 
    148                             top = A_ncols 
    149                         for k from k <= k < top: # = min(k+gather, A._ncols) 
    150                             s += A_row[k] * B_matrix_off[k][j+B._col] 
    151                         s %= p 
    152                     self_row[j] = s 
     130        for i from 0 <= i < A._nrows: 
     131            for j from 0 <= j < B._ncols: 
     132                s = 0 
     133                for k from 0 <= k < A._ncols: 
     134                    s = (s + A_matrix[i+A._row][k+A._col] * B_matrix[k+B._row][j+B._col]) % p 
     135                self_matrix[i+self._row][j+self._col] = s % p 
    153136         
    154137    cdef add_prod(self, MatrixWindow A, MatrixWindow B): 
    155138        # TODO: gather 
    156         cdef Py_ssize_t i, j, k, A_ncols 
    157         cdef mod_int p, s 
    158         cdef mod_int* self_row 
    159         cdef mod_int* A_row 
    160         cdef mod_int* B_row 
     139        cdef Py_ssize_t i, j, k 
     140        cdef uint p, s 
     141        cdef uint** self_matrix 
     142        cdef uint** A_matrix 
     143        cdef uint** B_matrix 
    161144        if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols: 
    162145            raise ArithmeticError, "incompatible dimensions" 
     146        self_matrix = ( <Matrix_modn_dense> self._matrix ).matrix 
     147        A_matrix = ( <Matrix_modn_dense> A._matrix ).matrix 
     148        B_matrix = ( <Matrix_modn_dense> B._matrix ).matrix 
    163149        p = ( <Matrix_modn_dense> self._matrix ).p 
    164         gather = ( <Matrix_modn_dense> self._matrix ).gather 
    165         A_ncols = A._ncols 
    166          
    167         if gather <= 1: 
    168             for i from 0 <= i < A._nrows: 
    169                 self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col 
    170                 A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col 
    171                 B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col 
    172                 for j from 0 <= j < B._ncols: 
    173                     s = self_row[j] 
    174                     for k from 0 <= k < A._ncols: 
    175                         s = ( s + A_row[k] * B_row[j] ) % p 
    176                     self_row[j] = s 
    177                      
    178         else: 
    179             for i from 0 <= i < A._nrows: 
    180                 self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col 
    181                 A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col 
    182                 B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col 
    183                 for j from 0 <= j < B._ncols: 
    184                     s = self_row[j] 
    185                     k = 0 
    186                     while k < A_ncols: 
    187                         top = k + gather 
    188                         if top > A_ncols: 
    189                             top = A_ncols 
    190                         for k from k <= k < top: # = min(k+gather, A._ncols) 
    191                             s += A_row[k] * B_matrix_off[k][j+B._col] 
    192                         s %= p 
    193                     self_row[j] = s 
    194                      
     150        for i from 0 <= i < A._nrows: 
     151            for j from 0 <= j < B._ncols: 
     152                s = self_matrix[i+self._row][j+self._col] 
     153                for k from 0 <= k < A._ncols: 
     154                    s = ( s + A_matrix[i+A._row][k+A._col] * B_matrix[k+B._row][j+B._col] ) % p 
     155                self_matrix[i+self._row][j+self._col] = s % p 
    195156                 
    196157    cdef subtract_prod(self, MatrixWindow A, MatrixWindow B): 
    197158        # TODO: gather 
    198         cdef Py_ssize_t i, j, k, A_ncols 
    199         cdef mod_int p, s 
    200         cdef mod_int* self_row 
    201         cdef mod_int* A_row 
    202         cdef mod_int* B_row 
     159        cdef Py_ssize_t i, j, k 
     160        cdef uint p, s 
     161        cdef uint** self_matrix 
     162        cdef uint** A_matrix 
     163        cdef uint** B_matrix 
    203164        if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols: 
    204165            raise ArithmeticError, "incompatible dimensions" 
     166        self_matrix = ( <Matrix_modn_dense> self._matrix ).matrix 
     167        A_matrix = ( <Matrix_modn_dense> A._matrix ).matrix 
     168        B_matrix = ( <Matrix_modn_dense> B._matrix ).matrix 
    205169        p = ( <Matrix_modn_dense> self._matrix ).p 
    206         A_ncols = A._ncols 
    207          
    208         if gather <= 1: 
    209             for i from 0 <= i < A._nrows: 
    210                 self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col 
    211                 A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col 
    212                 B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col 
    213                 for j from 0 <= j < B._ncols: 
    214                     s = self_row[j] 
    215                     for k from 0 <= k < A._ncols: 
    216                         s = s + p - ( A_row[k] * B_row[j] ) % p 
    217                     self_row[j] = s 
    218                      
    219         else: 
    220             for i from 0 <= i < A._nrows: 
    221                 self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col 
    222                 A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col 
    223                 B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col 
    224                 for j from 0 <= j < B._ncols: 
    225                     s = self_row[j] + p * gather # because they're unsigned 
    226                     k = 0 
    227                     while k < A_ncols: 
    228                         top = k + gather 
    229                         if top > A_ncols: 
    230                             top = A_ncols 
    231                         for k from k <= k < top: # = min(k+gather, A._ncols) 
    232                             s -= A_row[k] * B_matrix_off[k][j+B._col] 
    233                         s = s % p  +  p * gather # because they're unsigne 
    234                     self_row[j] = s % p 
     170        for i from 0 <= i < A._nrows: 
     171            for j from 0 <= j < B._ncols: 
     172                s = self_matrix[i+self._row][j+self._col] 
     173                for k from 0 <= k < A._ncols: 
     174                    s = s + p - ( A_matrix[i+A._row][k+A._col] * B_matrix[k+B._row][j+B._col] ) % p 
     175                self_matrix[i+self._row][j+self._col] = s % p 
    235176 
    236177    cdef int element_is_zero(self, Py_ssize_t i, Py_ssize_t j): 

  • setup.py

    r2847 r2632  
    314314    Extension('sage.ext.arith_gmp', 
    315315              sources = ['sage/ext/arith_gmp.pyx'], 
    316               libraries=['gmp']), \ 
    317  
    318     Extension('sage.ext.multi_modular', 
    319               sources = ['sage/ext/multi_modular.pyx'], 
    320316              libraries=['gmp']), \ 
    321317 
Note: See TracChangeset for help on using the changeset viewer.