Changeset 3260:bafb6aef9dad


Ignore:
Timestamp:
03/03/07 11:36:34 (6 years ago)
Author:
William Stein <wstein@…>
Branch:
default
Message:

prelimnary partial implementation of decomposition

Location:
sage
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • sage/matrix/benchmark.py

    r3103 r3260  
    593593        A = random_matrix(ZZ, n, n, x=min, y=max+1).change_ring(QQ) 
    594594        t = cputime() 
    595         v = A**(-1) 
     595        v = ~A 
    596596        return cputime(t) 
    597597    elif system == 'magma': 

  • sage/matrix/matrix0.pyx

    r3254 r3260  
    11941194 
    11951195        INPUT: 
    1196             v -- list 
     1196            v -- list of length at most the number of rows of self. 
    11971197         
    11981198        EXAMPLES: 
    … …  
    12081208        R = self.rows() 
    12091209        n = len(R) 
    1210         if len(v) != n: 
    1211             raise ValueError, "length of v must equal number of rows." 
     1210        if len(v) > n: 
     1211            raise ValueError, "length of v (=%s) must be at most the number (=%s) of rows."%(len(v), n) 
    12121212        zero = self._base_ring(0) 
    12131213        s = None 
    1214         for i from 0 <= i < n: 
     1214        for i from 0 <= i < len(v): 
    12151215            if v[i] != zero: 
    12161216                a = v[i] * R[i] 
    … …  
    12181218                    s = a 
    12191219                else: 
    1220                     s = s + a 
     1220                    s += a 
    12211221        if s is None: 
    12221222            return self.parent().row_space()(0) 

  • sage/matrix/matrix2.pyx

    r3222 r3260  
    11381138        """ 
    11391139        M = self._row_ambient_module() 
    1140         return M.span(self.rows()) 
     1140        if self.fetch('in_echelon_form'): 
     1141            return M.span(self.rows(), already_echelonized=True) 
     1142        else: 
     1143            return M.span(self.rows(), already_echelonized=True) 
    11411144 
    11421145    def row_space(self): 
    … …  
    12161219        """ 
    12171220        Returns the decomposition of the free module on which this 
    1218         matrix acts from the right, along with whether this matrix 
    1219         acts irreducibly on each factor.  The factors are guaranteed 
    1220         to be sorted in the same way as the corresponding factors of 
    1221         the characteristic polynomial. 
     1221        matrix A acts from the right (i.e., the action is x goes to x 
     1222        A), along with whether this matrix acts irreducibly on each 
     1223        factor.  The factors are guaranteed to be sorted in the same 
     1224        way as the corresponding factors of the characteristic 
     1225        polynomial. 
    12221226 
    12231227        Let A be the matrix acting from the on the vector space V of 
    … …  
    13021306        F.sort() 
    13031307        for g, m in f.factor(): 
     1308            t = verbose('decomposition -- Computing g(self) for an irreducible factor g of degree %s'%g.degree(),level=2) 
    13041309            if is_diagonalizable: 
    13051310                B = g(self) 
    13061311            else: 
    1307                 B = g(self) ** m 
     1312                B = g(self) 
     1313                t2 = verbose('decomposition -- raising g(self) to the power %s'%m,level=2) 
     1314                B = B ** m 
     1315                verbose('done powering',t2) 
     1316            t = verbose('decomposition -- done computing g(self)', level=2, t=t) 
    13081317            E.append((B.kernel(), m==1)) 
     1318            t = verbose('decomposition -- time to compute kernel', level=2, t=t) 
    13091319            if dual: 
    13101320                Edual.append((B.transpose().kernel(), m==1)) 
     1321                verbose('decomposition -- time to compute dual kernel', level=2, t=t) 
    13111322        if dual: 
    13121323            return E, Edual 
     1324        return E 
     1325 
     1326    def _decomposition_devel(self, is_diagonalizable=False, dual=False, 
     1327                            echelon_algorithm='default', height_guess=None): 
     1328        """ 
     1329        Returns the decomposition of the free module on which this 
     1330        matrix A acts from the right (i.e., the action is x goes to x 
     1331        A), along with whether this matrix acts irreducibly on each 
     1332        factor.  The factors are guaranteed to be sorted in the same 
     1333        way as the corresponding factors of the characteristic 
     1334        polynomial. 
     1335 
     1336        Let A be the matrix acting from the on the vector space V of 
     1337        column vectors.  Assume that A is square.  This function 
     1338        computes maximal subspaces W_1, ..., W_n corresponding to 
     1339        Galois conjugacy classes of eigenvalues of A.  More precisely, 
     1340        let f(X) be the characteristic polynomial of A.  This function 
     1341        computes the subspace $W_i = ker(g_(A)^n)$, where g_i(X) is an 
     1342        irreducible factor of f(X) and g_i(X) exactly divides f(X). 
     1343        If the optional parameter is_diagonalizable is True, then we 
     1344        let W_i = ker(g(A)), since then we know that ker(g(A)) = 
     1345        $ker(g(A)^n)$. 
     1346 
     1347        If dual is True, also returns the corresponding decomposition 
     1348        of V under the action of the transpose of A.  The factors are 
     1349        guarenteed to correspond. 
     1350 
     1351        INPUT: 
     1352            self -- a matrix over a field 
     1353 
     1354        OUTPUT: 
     1355            Sequence -- list of pairs (V,t), where V is a vector spaces 
     1356                    and t is a bool, and t is True exactly when the 
     1357                    charpoly of self on V is irreducible. 
     1358 
     1359            (optional) list -- list of pairs (W,t), where W is a vector 
     1360                    space and t is a bool, and t is True exactly 
     1361                    when the charpoly of the transpose of self on W 
     1362                    is irreducible. 
     1363 
     1364        EXAMPLES: 
     1365            sage: MS1 = MatrixSpace(ZZ,4) 
     1366            sage: MS2 = MatrixSpace(QQ,6) 
     1367            sage: A = MS1.matrix([3,4,5,6,7,3,8,10,14,5,6,7,2,2,10,9]) 
     1368            sage: B = MS2(range(36)) 
     1369            sage: B*11   # random output 
     1370            [-11  22 -11 -11 -11 -11] 
     1371            [ 11 -22 -11 -22  11  11] 
     1372            [-11 -11 -11 -22 -22 -11] 
     1373            [-22  22 -22  22 -11  11] 
     1374            [ 22 -11  11 -22  11  22] 
     1375            [ 11  11  11 -22  22  22] 
     1376            sage: decomposition(A)  
     1377            [ 
     1378            (Ambient free module of rank 4 over the principal ideal domain Integer Ring, 1) 
     1379            ]            
     1380            sage: decomposition(B) 
     1381            [ 
     1382            (Vector space of degree 6 and dimension 2 over Rational Field 
     1383            Basis matrix: 
     1384            [ 1  0 -1 -2 -3 -4] 
     1385            [ 0  1  2  3  4  5], 1), 
     1386            (Vector space of degree 6 and dimension 4 over Rational Field 
     1387            Basis matrix: 
     1388            [ 1  0  0  0 -5  4] 
     1389            [ 0  1  0  0 -4  3] 
     1390            [ 0  0  1  0 -3  2] 
     1391            [ 0  0  0  1 -2  1], 0) 
     1392            ] 
     1393        """ 
     1394        if not self.is_square(): 
     1395            raise ArithmeticError, "self must be a square matrix" 
     1396 
     1397        if not self.base_ring().is_field(): 
     1398            raise TypeError, "self must be over a field." 
     1399 
     1400        if self.nrows() == 0: 
     1401            return decomp_seq([]) 
     1402 
     1403        f = self.charpoly('x') 
     1404        E = decomp_seq([]) 
     1405 
     1406        if dual: 
     1407            Edual = decomp_seq([]) 
     1408 
     1409        t = verbose('decomposition -- factoring the characteristic polynomial', level=2) 
     1410        F = f.factor() 
     1411        verbose('decomposition -- done factoring', t=t, level=2) 
     1412         
     1413        V = sage.modules.free_module.FreeModule( 
     1414                self.base_ring(), self.nrows(), sparse=self.is_sparse()) 
     1415         
     1416        if len(F) == 1: 
     1417            m = F[0][1] 
     1418            if dual: 
     1419                return decomp_seq([(V,m==1)]), decomp_seq([(V,m==1)]) 
     1420            else: 
     1421                return decomp_seq([(V,m==1)]) 
     1422 
     1423        v = V.random_element() 
     1424         
     1425        num_iterates = max([f.degree() - g.degree() for g, _ in F]) + 1 
     1426        t = verbose('decomposition -- computing %s iterates of a random vector.'%num_iterates) 
     1427         
     1428        S = self.iterates(v, num_iterates) 
     1429         
     1430        verbose('decomposition -- done computing iterates.', t = t) 
     1431             
     1432        F.sort() 
     1433        for i in range(len(F)): 
     1434            g, m = F[i] 
     1435 
     1436            # Compute the complementary factor. 
     1437            h = f // (g**m) 
     1438 
     1439            # Compute one element of the kernel of g(self)**m. 
     1440            t = verbose('decomposition -- compute element of kernel of g(self), for g of degree %s'%g.degree(),level=2) 
     1441            w = S.linear_combination_of_rows(h.list()) 
     1442            t = verbose('decomposition -- done computing element of kernel of g(self)', t=t,level=2) 
     1443 
     1444            # normalize?? 
     1445            t = verbose('decomposition -- normalize vector', level=2) 
     1446            w = w.normalize() 
     1447            verbose('decomposition -- normalize vector', t, level = 2) 
     1448 
     1449            # Get the rest of the kernel. 
     1450            t = verbose('decomposition -- fill out rest of kernel',level=2) 
     1451            W = self.iterates(w, g.degree()) 
     1452            t = verbose('decomposition -- finished filling out rest of kernel',level=2, t=t) 
     1453 
     1454            t = verbose('decomposition -- now computing row space', level=2) 
     1455            W.echelonize(algorithm = echelon_algorithm, height_guess=height_guess) 
     1456            E.append((W.row_space(), m==1)) 
     1457            verbose('decomposition -- computed row space', level=2,t=t) 
     1458             
    13131459        return E 
    13141460 

  • sage/matrix/matrix_integer_dense.pyx

    r3259 r3260  
    17061706        cdef Matrix_integer_dense M 
    17071707        cdef Integer D 
    1708          
     1708             
    17091709        if self._nrows != self._ncols: 
    17101710            raise ArithmeticError, "self must be square" 
    … …  
    18531853        A = self 
    18541854        # Step 1: Compute the rank 
    1855          
     1855 
     1856        t = verbose('computing rank', level=2, caller_name='p-adic echelon') 
    18561857        r = A.rank() 
     1858        verbose('done computing rank', level=2, t=t, caller_name='p-adic echelon') 
     1859         
    18571860        if r == self._nrows: 
    18581861            # The input matrix already has full rank.  
    … …  
    18951898            non_pivots = [i for i in range(B.ncols()) if not i in pivots_] 
    18961899            D = B.matrix_from_columns(non_pivots) 
     1900            t = verbose('calling IML solver', level=2, caller_name='p-adic echelon') 
    18971901            X, d = C._solve_iml(D, right=True) 
     1902            t = verbose('finished IML solver', level=2, caller_name='p-adic echelon', t=t) 
    18981903 
    18991904            # Step 6: Verify that we had chosen the correct pivot columns. 

  • sage/matrix/matrix_rational_dense.pyx

    r3259 r3260  
    5050  
    5151cimport sage.structure.element 
     52 
     53from sage.structure.sequence import Sequence 
    5254from sage.rings.rational cimport Rational 
    5355from matrix cimport Matrix 
    … …  
    10311033                                 height_guess=height_guess, proof=proof) 
    10321034     
     1035    def _decomposition_devel(self, is_diagonalizable=False, dual=False, 
     1036                            echelon_algorithm='default', height_guess=None): 
     1037        """ 
     1038        Returns the decomposition of the free module on which this 
     1039        matrix A acts from the right (i.e., the action is x goes to x 
     1040        A), along with whether this matrix acts irreducibly on each 
     1041        factor.  The factors are guaranteed to be sorted in the same 
     1042        way as the corresponding factors of the characteristic 
     1043        polynomial. 
     1044 
     1045        INPUT: 
     1046            self -- a matrix over a field 
     1047 
     1048        OUTPUT: 
     1049            Sequence -- list of pairs (V,t), where V is a vector spaces 
     1050                    and t is a bool, and t is True exactly when the 
     1051                    charpoly of self on V is irreducible. 
     1052 
     1053            (optional) list -- list of pairs (W,t), where W is a vector 
     1054                    space and t is a bool, and t is True exactly 
     1055                    when the charpoly of the transpose of self on W 
     1056                    is irreducible. 
     1057        """ 
     1058        if not self.is_square(): 
     1059            raise ArithmeticError, "self must be a square matrix" 
     1060 
     1061        if not self.base_ring().is_field(): 
     1062            raise TypeError, "self must be over a field." 
     1063 
     1064        if self.nrows() == 0: 
     1065            return decomp_seq([]) 
     1066 
     1067        f = self.charpoly('x') 
     1068        E = decomp_seq([]) 
     1069 
     1070        if dual: 
     1071            Edual = decomp_seq([]) 
     1072 
     1073        t = verbose('decomposition -- factoring the characteristic polynomial', level=2) 
     1074        F = f.factor() 
     1075        verbose('decomposition -- done factoring', t=t, level=2) 
     1076         
     1077        if len(F) == 1: 
     1078            V = QQ**self.nrows() 
     1079            m = F[0][1] 
     1080            if dual: 
     1081                return decomp_seq([(V,m==1)]), decomp_seq([(V,m==1)]) 
     1082            else: 
     1083                return decomp_seq([(V,m==1)]) 
     1084 
     1085 
     1086        V = ZZ**self.nrows() 
     1087        v = V.random_element() 
     1088         
     1089        num_iterates = max([f.degree() - g.degree() for g, _ in F]) + 1 
     1090 
     1091        A, _ = self._clear_denom() 
     1092         
     1093        t = verbose('decomposition -- computing %s iterates of a random vector.'%num_iterates) 
     1094        S = A.iterates(v, num_iterates) 
     1095        verbose('decomposition -- done computing iterates.', t = t) 
     1096             
     1097        F.sort() 
     1098        for i in range(len(F)): 
     1099            g, m = F[i] 
     1100 
     1101            # Compute the complementary factor. 
     1102            h = f // (g**m) 
     1103            h = h * h.denominator() 
     1104 
     1105            # Compute one element of the kernel of g(self)**m. 
     1106            t = verbose('decomposition -- compute element of kernel of g(self), for g of degree %s'%g.degree(),level=2) 
     1107            w = S.linear_combination_of_rows(h.list()) 
     1108            t = verbose('decomposition -- done computing element of kernel of g(self)', t=t,level=2) 
     1109 
     1110            # Get the rest of the kernel. 
     1111            t = verbose('decomposition -- fill out rest of kernel',level=2) 
     1112            W = A.iterates(w, g.degree()) 
     1113            t = verbose('decomposition -- finished filling out rest of kernel',level=2, t=t) 
     1114 
     1115            W = W.change_ring(QQ) 
     1116 
     1117            t = verbose('decomposition -- now computing row space', level=2) 
     1118            W.echelonize(algorithm = echelon_algorithm, height_guess=height_guess) 
     1119            E.append((W.row_space(), m==1)) 
     1120            verbose('decomposition -- computed row space', level=2,t=t) 
     1121             
     1122        return E 
     1123 
    10331124    ################################################################ 
    10341125    # second implementation of the above, usually over twice as fast 
    … …  
    13181409 
    13191410             
     1411 
     1412 
     1413cdef decomp_seq(v): 
     1414    return Sequence(v, universe=tuple, check=False, cr=True) 

  • sage/modules/free_module_element.pyx

    r3182 r3260  
    555555        self[i] = x 
    556556 
     557 
     558    def normalize(self): 
     559        """ 
     560        Return this vector divided through by the first nonzero entry of this vector. 
     561 
     562        EXAMPLES: 
     563            sage: v = vector(QQ,[0,4/3,5,1,2]) 
     564            sage: v.normalize() 
     565            (0, 1, 15/4, 3/4, 3/2) 
     566        """ 
     567        cdef Py_ssize_t i 
     568        for i from 0 <= i < self._degree: 
     569            if self[i] != 0: 
     570                return (~self[i]) * self 
     571        return self 
    557572         
    558573    def inner_product(self, right): 
Note: See TracChangeset for help on using the changeset viewer.