Changeset 2ed3a3 in git for Singular/LIB/matrix.lib
 Timestamp:
 Nov 20, 2008, 6:17:33 PM (15 years ago)
 Branches:
 (u'spielwiese', 'ec94ef7a30b928574c0c3daf41f6804dff5f6b69')
 Children:
 c16f9922fe5d93bd72e4b2440cd6c7cd43889ca6
 Parents:
 2a03adf5a9e042150758e454f9153509b05e3633
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/matrix.lib
r2a03ad r2ed3a3 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: matrix.lib,v 1.4 3 20081001 15:29:23 SingularExp $";2 version="$Id: matrix.lib,v 1.44 20081120 17:17:33 motsak Exp $"; 3 3 category="Linear Algebra"; 4 4 info=" … … 34 34 rm_unitcol(A); remove unit columns and associated rows of A 35 35 headStand(A); A[ni+1,mj+1]:=A[i,j] 36 symmetricBasis(n,k ); ring, with a basis of a kth symmetric power of anndim v.space37 symmetricPower(A,k); module, kth symmetric power of a module/matrix A38 exteriorBasis(n,k); ring, with a basis of a kth exterior power of an ndim v.space39 exteriorPower(A,k); module,kth exterior power of a module/matrix A36 symmetricBasis(n,k[,s]); basis of kth symmetric power of ndim v.space 37 exteriorBasis(n,k[,s]); basis of kth exterior power of ndim v.space 38 symmetricPower(A,k); kth symmetric power of a module/matrix A 39 exteriorPower(A,k); kth exterior power of a module/matrix A 40 40 (parameters in square brackets [] are optional) 41 41 "; … … 1074 1074 ////////////////////////////////////////////////////////////////////////////// 1075 1075 1076 // Symmetric/Exterior powers, thanks to Oleksandr Iena for his persistence ;) 1077 1078 proc symmetricBasis(int n, int k) 1079 "USAGE: symmetricBasis(n, k); n int, k int 1080 RETURN: ring: the basis of the kth symmetric power of a ndim vector space 1081 NOTE: The output consists of a polynomial ring in n variables, together with the ideal called I. 1082 The ideal I is the basis of the kth symmetric power of a ndim vector space (ordered lexicographically). 1083 The char. of the returned ring is the same as of current basis ring or zero. 1076 // Symmetric/Exterior powers thanks to Oleksandr Iena for his persistence ;) 1077 1078 proc symmetricBasis(int n, int k, list #) 1079 "USAGE: symmetricBasis(n, k[,s]); n int, k int, s string 1080 RETURN: poynomial ring containing the ideal \"symBasis\", 1081 being a basis of the kth symmetric power of an ndim vector space. 1082 NOTE: The output polynomial ring has characteristics 0 and n variables 1083 named \"S(i)\", where the base variable name S is either given by the 1084 optional string argument(which must not contain brackets) or equal to 1085 "e" by default. 1084 1086 SEE ALSO: exteriorBasis 1085 1087 KEYWORDS: symmetric basis 1086 1088 EXAMPLE: example symmetricBasis; shows an example" 1087 1089 { 1088 int p = 0; 1089 if( nameof(basering) != "basering" ) 1090 // handle optional base variable name 1091 string S = "e"; 1092 if( size(#) > 0 ) 1090 1093 { 1091 p = char(basering); 1094 if( typeof(#[1]) != "string" ) 1095 { 1096 ERROR("Wrong optional argument: must be a string"); 1097 } 1098 S = #[1]; 1099 if( (find(S, "(") + find(S, ")")) > 0 ) 1100 { 1101 ERROR("Wrong optional argument: must be a string without brackets"); 1102 } 1092 1103 } 1093 1104 1094 ring T = (p), (e(1..n)), dp; 1095 ideal I = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 ); 1096 int N = ncols(I); 1097 I = sort(I)[1]; // lex 1098 I = I[N..1]; 1099 1100 export I; 1101 return(T); 1102 } 1103 example 1104 { "EXAMPLE:"; echo = 2; 1105 1106 def r = symmetricBasis(2, 3); setring r; r; 1107 I; // basis of 3rd sym. power of a 2dim v.s. 1108 kill r; 1109 1110 ring r = (0),(a, b, c, d), dp; r; 1111 def g = symmetricBasis(3, 2); setring g; g; I; 1112 1113 kill g, r; 1114 1115 ring r = (32003),(a, b, c, d), dp; r; 1116 def g = symmetricBasis(4, 2); setring g; g; I; 1105 // create ring container for symmetric power basis 1106 execute("ring @@@SYM_POWER_RING_NAME=(0),("+S+"(1.."+string(n)+")),dp;"); 1107 1108 // choose symmetric basis  1109 ideal symBasis = maxideal(k); 1110 1111 // export and return  1112 export symBasis; 1113 return(basering); 1114 } 1115 example 1116 { "EXAMPLE:"; echo = 2; 1117 1118 // basis of the 3rd symmetricPower of a 4dim vector space: 1119 def R = symmetricBasis(4, 3, "@e"); setring R; 1120 R; // container ring: 1121 symBasis; // symmetric basis: 1117 1122 } 1118 1123 1119 1124 ////////////////////////////////////////////////////////////////////////////// 1120 1125 1121 proc exteriorBasis(int n, int k) 1122 "USAGE: exteriorBasis(n, k); n int, k int 1123 RETURN: ring: the basis of the kth exterior power of a ndim vector space 1124 NOTE: The output consists of a polynomial ring in n variables, together with the ideal called I. 1125 The ideal I is the basis of the kth exterior power of a ndim vector space (ordered lexicographically). 1126 The char. of the returned ring is the same as of current basis ring or zero. 1127 1126 proc exteriorBasis(int n, int k, list #) 1127 "USAGE: exteriorBasis(n, k[,s]); n int, k int, s string 1128 RETURN: poynomial ring containing the ideal \"extBasis\", 1129 being a basis of the kth exterior power of an ndim vector space. 1130 NOTE: The output polynomial ring has characteristics 0 and n variables 1131 named \"S(i)\", where the base variable name S is either given by the 1132 optional string argument(which must not contain brackets) or equal to 1133 "e" by default. 1128 1134 SEE ALSO: symmetricBasis 1129 1135 KEYWORDS: exterior basis 1130 1136 EXAMPLE: example exteriorBasis; shows an example" 1131 1137 { 1132 int p = 0; 1133 if( nameof(basering) != "basering" ) 1138 // handle optional base variable name 1139 string S = "e"; 1140 if( size(#) > 0 ) 1134 1141 { 1135 p = char(basering); 1142 if( typeof(#[1]) != "string" ) 1143 { 1144 ERROR("Wrong optional argument: must be a string"); 1145 } 1146 S = #[1]; 1147 if( (find(S, "(") + find(S, ")")) > 0 ) 1148 { 1149 ERROR("Wrong optional argument: must be a string without brackets"); 1150 } 1136 1151 } 1137 1138 ring S = (p), (e(1..n)), dp; 1152 1153 // create ring container for symmetric power basis 1154 execute("ring @@@EXT_POWER_RING_NAME=(0),("+S+"(1.."+string(n)+")),dp;"); 1155 1156 // choose exterior basis  1139 1157 def T = SuperCommutative(); setring T; 1140 ideal I = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 ); 1141 int N = ncols(I); 1142 I = sort(I)[1]; 1143 I = I[N..1 ]; 1144 1145 export I; 1146 return(T); 1147 } 1148 example 1149 { "EXAMPLE:"; echo = 2; 1150 1151 def r = exteriorBasis(2, 3); setring r; r; 1152 "Basis: ", I; // basis of 3rd sym. power of a 2dim v.s. 1153 kill r; 1154 1155 ring r = (0),(a, b, c, d), dp; r; 1156 def g = exteriorBasis(3, 2); setring g; g; "Basis: ", I; 1157 1158 kill g, r; 1159 1160 ring r = (32003),(a, b, c, d), dp; r; 1161 def g = exteriorBasis(4, 2); setring g; g; "Basis: ", I; 1162 } 1163 1164 ////////////////////////////////////////////////////////////////////////////// 1165 1166 proc symmetricPower(module A, int k) 1167 "USAGE: symmetricPower(A, k); A module, k int 1168 RETURN: module: the kth symmetric power of A 1169 NOTE: the chosen bases (ordered lexicographically) and 1170 temporary data may be shown if printlevel is big enough 1171 SEE ALSO: exteriorPower 1172 KEYWORDS: symmetric power 1173 EXAMPLE: example symmetricPower; shows an example" 1174 { 1175 int p = printlevel  voice + 2; 1158 ideal extBasis = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 ); 1159 1160 // export and return  1161 export extBasis; 1162 return(basering); 1163 } 1164 example 1165 { "EXAMPLE:"; echo = 2; 1166 // basis of the 3rd symmetricPower of a 4dim vector space: 1167 def r = exteriorBasis(4, 3, "@e"); setring r; 1168 r; // container ring: 1169 extBasis; // exterior basis: 1170 } 1171 1172 1173 static proc chooseSafeVarName(string prefix, string suffix) 1174 "USAGE: give appropreate prefix for variable names 1175 RETURN: safe variable name (repeated prefix + suffix) 1176 " 1177 { 1178 string V = varstr(basering); 1179 string S = suffix; 1180 while( find(V, S) > 0 ) 1181 { 1182 S = prefix + S; 1183 } 1184 return(S); 1185 } 1186 1187 static proc mapPower(int p, module A, int k, def Tn, def Tm) 1188 "USAGE: by both symmetric and exteriorPower" 1189 NOTE: everything over the basering! 1190 module A (matrix of the map), int k (power) 1191 rings Tn is source and Tm is imagering with bases 1192 resp. Ink and Imk. 1193 M = max dim of Image, N  dim. of source 1194 SEE ALSO: symmetricPower, exteriorPower" 1195 { 1176 1196 def save = basering; 1177 1197 1198 int n = nvars(save); 1178 1199 int M = nrows(A); 1179 1200 int N = ncols(A); … … 1181 1202 int i, j; 1182 1203 1183 ///////////////////////////////////////////////////////////////// 1184 def Tn = symmetricBasis(N, k); setring Tn; 1185 ideal Ink = I; 1186 dbprint(p3, "Temporary Source Ring", basering); 1187 dbprint(p, "S^k(Source Basis)", Ink); 1188 1189 ///////////////////////////////////////////////////////////////// 1190 def Tm = symmetricBasis(M, k); setring Tm; 1191 ideal Imk = I; 1192 dbprint(p3, "Temporary Image Ring", basering); 1193 dbprint(p, "S^k(Image Basis)", Imk); 1194 1195 ///////////////////////////////////////////////////////////////// 1196 def Rm = save + Tm; 1197 setring Rm; 1204 // compute matrix of single images  1205 def Rm = save + Tm; setring Rm; 1198 1206 dbprint(p2, "Temporary Working Ring", Rm); 1199 1207 … … 1207 1215 for( j = M; j > 0; j ) 1208 1216 { 1209 t = t + A[i][j] * var( nvars(save) + j); // tensor product!!!1217 t = t + A[i][j] * var(n + j); 1210 1218 } 1211 1219 … … 1213 1221 } 1214 1222 1215 dbprint(p1, "Matrix of simgle images", B); 1216 1217 map TMap = Tn, B; ideal C = TMap(Ink); // apply S^k A to Source basis vectors... (Ink) 1218 C = NF(C, std(0)); 1219 1223 dbprint(p1, "Matrix of single images", B); 1224 1225 // compute image  1226 // apply S^k(A): Tn > Rm to Source basis vectors Ink: 1227 map TMap = Tn, B; 1228 1229 ideal C = NF(TMap(Ink), std(0)); 1220 1230 dbprint(p1, "Image Matrix: ", C); 1221 1231 1222 // and write it in Image basis (Imk) 1223 1232 1233 // write it in Image basis  1224 1234 ideal Imk = imap(Tm, Imk); 1225 1235 … … 1247 1257 1248 1258 D[i] = tt; 1249 // tt;1250 1259 } 1251 1260 1252 // D; print(D); 1253 1254 1261 // map it back and return  1255 1262 setring save; 1256 1257 // basering; 1258 1259 module AA = imap(Rm, D); 1260 1261 // nrows(AA) = Nk; // ???? 1262 // ncols(AA) = Mk; 1263 1264 // Nk, Mk; 1265 1266 // AA[Mk] = AA[Mk] * 1; 1267 1268 return(AA); 1269 } 1270 example 1271 { "EXAMPLE:"; echo = 2; 1272 int save = printlevel; printlevel = 1; 1263 return( imap(Rm, D) ); 1264 } 1265 1266 1267 1268 1269 ////////////////////////////////////////////////////////////////////////////// 1270 1271 proc symmetricPower(module A, int k) 1272 "USAGE: symmetricPower(A, k); A module, k int 1273 RETURN: module: the kth symmetric power of A 1274 NOTE: the chosen bases and most of intermediate data will be shown if 1275 printlevel is big enough 1276 SEE ALSO: exteriorPower 1277 KEYWORDS: symmetric power 1278 EXAMPLE: example symmetricPower; shows an example" 1279 { 1280 int p = printlevel  voice + 2; 1281 1282 def save = basering; 1283 1284 int M = nrows(A); 1285 int N = ncols(A); 1286 1287 string S = chooseSafeVarName("@", "@_e"); 1288 1289 // choose source basis  1290 def Tn = symmetricBasis(N, k, S); setring Tn; 1291 ideal Ink = symBasis; 1292 export Ink; 1293 dbprint(p3, "Temporary Source Ring", basering); 1294 dbprint(p, "S^k(Source Basis)", Ink); 1295 1296 // choose image basis  1297 def Tm = symmetricBasis(M, k, S); setring Tm; 1298 ideal Imk = symBasis; 1299 export Imk; 1300 dbprint(p3, "Temporary Image Ring", basering); 1301 dbprint(p, "S^k(Image Basis)", Imk); 1302 1303 // compute and return S^k(A) in chosen bases  1304 setring save; 1305 1306 return(mapPower(p, A, k, Tn, Tm)); 1307 } 1308 example 1309 { "EXAMPLE:"; echo = 2; 1273 1310 1274 1311 ring r = (0),(a, b, c, d), dp; r; 1275 1276 module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1); 1277 1278 print(symmetricPower(A, 2)); 1279 print(symmetricPower(A, 3)); 1280 1281 module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); 1282 1312 module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); print(B); 1313 1314 // symmetric power over a commutative Kalgebra: 1283 1315 print(symmetricPower(B, 2)); 1284 1316 print(symmetricPower(B, 3)); 1285 1317 1286 kill r; 1287 ring r = (0),(a, b, c, d), dp;def g = SuperCommutative();setring g; 1288 1289 module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1); 1290 1291 print(symmetricPower(A, 2)); 1292 print(symmetricPower(A, 3)); 1293 1294 module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); 1295 1296 print(symmetricPower(B, 2)); 1297 print(symmetricPower(B, 3)); 1298 1299 kill r; 1300 1301 printlevel = save; 1318 // symmetric power over an exterior algebra: 1319 def g = SuperCommutative(); setring g; g; 1320 1321 module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); print(B); 1322 1323 print(symmetricPower(B, 2)); // much smaller! 1324 print(symmetricPower(B, 3)); // zero! (over an exterior algebra!) 1325 1302 1326 } 1303 1327 … … 1307 1331 "USAGE: exteriorPower(A, k); A module, k int 1308 1332 RETURN: module: the kth exterior power of A 1309 NOTE: the chosen bases and temporary data 1310 may be shown if printlevel is big enough. 1311 Last rows may be invisible if zero. 1333 NOTE: the chosen bases and most of intermediate data will be shown if 1334 printlevel is big enough. Last rows will be invisible if zero. 1312 1335 SEE ALSO: symmetricPower 1313 1336 KEYWORDS: exterior power … … 1320 1343 int N = ncols(A); 1321 1344 1322 int i, j; 1323 1324 ///////////////////////////////////////////////////////////////// 1325 def Tn = exteriorBasis(N, k); setring Tn; 1326 ideal Ink = I; 1345 string S = chooseSafeVarName("@", "@_e"); 1346 1347 // choose source basis  1348 def Tn = exteriorBasis(N, k, S); setring Tn; 1349 ideal Ink = extBasis; 1350 export Ink; 1327 1351 dbprint(p3, "Temporary Source Ring", basering); 1328 1352 dbprint(p, "E^k(Source Basis)", Ink); 1329 1353 1330 ///////////////////////////////////////////////////////////////// 1331 def Tm = exteriorBasis(M, k); setring Tm; 1332 ideal Imk = I; 1354 // choose image basis  1355 def Tm = exteriorBasis(M, k, S); setring Tm; 1356 ideal Imk = extBasis; 1357 export Imk; 1333 1358 dbprint(p3, "Temporary Image Ring", basering); 1334 1359 dbprint(p, "E^k(Image Basis)", Imk); 1335 1360 1336 1337 def Rm = save + Tm; 1338 setring Rm; 1339 dbprint(p2, "Temporary Working Ring", Rm); 1340 1341 module A = imap(save, A); 1342 1343 ideal B; poly t; 1344 1345 for( i = N; i > 0; i ) 1346 { 1347 t = 0; 1348 for( j = M; j > 0; j ) 1349 { 1350 t = t + A[i][j] * var( nvars(save) + j); // tensor product!!! 1351 } 1352 1353 B[i] = t; 1354 } 1355 1356 dbprint(p1, "Matrix of simgle images", B); 1357 1358 map TMap = Tn, B; ideal C = TMap(Ink); // apply S^k A to Source basis vectors... (Ink) 1359 C = NF(C, std(0)); 1360 1361 dbprint(p1, "Image Matrix: ", C); 1362 1363 // and write it in Image basis (Imk) 1364 1365 ideal Imk = imap(Tm, Imk); 1366 1367 module D; poly lm; vector tt; 1368 1369 for( i = ncols(C); i > 0; i ) 1370 { 1371 t = C[i]; 1372 tt = 0; 1373 1374 while( t != 0 ) 1375 { 1376 lm = leadmonom(t); 1377 // lm; 1378 for( j = ncols(Imk); j > 0; j ) 1379 { 1380 if( lm / Imk[j] != 0 ) 1381 { 1382 tt = tt + (lead(t) / Imk[j]) * gen(j); 1383 break; 1384 } 1385 } 1386 t = t  lead(t); 1387 } 1388 1389 D[i] = tt; 1390 // tt; 1391 } 1392 1393 // D; print(D); 1394 1395 1361 // compute and return E^k(A) in chosen bases  1396 1362 setring save; 1397 1398 // basering; 1399 1400 module AA = imap(Rm, D); 1401 1402 // nrows(AA) = Nk; // ???? 1403 // ncols(AA) = Mk; 1404 1405 // Nk, Mk; 1406 1407 // AA[Mk] = AA[Mk] * 1; 1408 1409 return(AA); 1410 } 1411 example 1412 { "EXAMPLE:"; echo = 2; 1413 int save = printlevel; printlevel = 1; 1414 1415 ring r = (0),(a, b, c, d), dp; r; 1416 1417 module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1); 1363 return(mapPower(p, A, k, Tn, Tm)); 1364 } 1365 example 1366 { "EXAMPLE:"; echo = 2; 1367 ring r = (0),(a, b, c, d, e, f), dp; 1368 r; "base ring:"; 1369 1370 module B = a*gen(1) + c*gen(2) + e*gen(3), 1371 b*gen(1) + d*gen(2) + f*gen(3), 1372 e*gen(1) + f*gen(3); 1373 1374 print(B); 1375 print(exteriorPower(B, 2)); 1376 print(exteriorPower(B, 3)); 1377 1378 def g = SuperCommutative(); setring g; g; 1379 1380 module A = a*gen(1), b * gen(1), c*gen(2), d * gen(2); 1381 print(A); 1418 1382 1419 1383 print(exteriorPower(A, 2)); 1420 print(exteriorPower(A, 3)); 1421 1422 module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); 1384 1385 module B = a*gen(1) + c*gen(2) + e*gen(3), 1386 b*gen(1) + d*gen(2) + f*gen(3), 1387 e*gen(1) + f*gen(3); 1388 print(B); 1423 1389 1424 1390 print(exteriorPower(B, 2)); 1425 1391 print(exteriorPower(B, 3)); 1426 1392 1427 kill r;1428 ring r = (0),(a, b, c, d), dp;def g = SuperCommutative();setring g;1429 1430 module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);1431 1432 print(exteriorPower(A, 2));1433 print(exteriorPower(A, 3));1434 1435 module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);1436 1437 print(exteriorPower(B, 2));1438 print(exteriorPower(B, 3));1439 1440 kill r;1441 1442 printlevel = save;1443 1393 } 1444 1394
Note: See TracChangeset
for help on using the changeset viewer.