Changeset fa85ef in git
- Timestamp:
- Oct 8, 2010, 9:44:27 PM (14 years ago)
- Branches:
- (u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '38077648e7239f98078663eb941c3c979511150a')
- Children:
- 1e33d5df4b2c03f3879e5b15e65274116be26b6b
- Parents:
- cc6f7a79d5d975a00b14db3835c4b6c0f13c5553
- File:
-
- 1 edited
-
Singular/LIB/ncfactor.lib (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/ncfactor.lib
rcc6f7a rfa85ef 3 3 category="Noncommutative"; 4 4 info=" 5 LIBRARY: ncfactor.lib Tools for factorization in the first Weyl algebra5 LIBRARY: ncfactor.lib Tools for factorization in some noncommutative algebras 6 6 AUTHORS: Albert Heinle, albert.heinle@rwth-aachen.de 7 @* Viktor Levandovskyy, levandov@math.rwth-aachen.de 8 OVERVIEW: In this library, new methods for factorization on polynomials 9 10 are implemented for two algebras over a field K. Recall, that 11 12 the first Weyl algebra over K is generated by x,d obeying the relation d*x=x*d+1. 13 14 The first shift algebra over K is generated by x,s obeying the relation s*x=x*s+s. 7 15 8 16 PROCEDURES: … … 10 18 testNCfac(l[,h]); tests factorizations from a given list for correctness 11 19 facSubWeyl(h,X,D); factorization in the first Weyl algebra as a subalgebra 20 facFirstShift(h); factorization in the first shift algebra 12 21 "; 13 22 … … 17 26 LIB "freegb.lib"; // for isVar 18 27 19 /* ring R = 0,(x,y),Ws(-1,1); */ 20 /* def r = nc_algebra(1,1); */ 21 /* setring(r); */ 28 proc tst_ncfactor() 29 { 30 example facFirstWeyl; 31 example facFirstShift; 32 example facSubWeyl; 33 example testNCfac; 34 } 22 35 23 36 ///////////////////////////////////////////////////// … … 27 40 static proc delete_dublicates_noteval(list l) 28 41 {//proc delete_dublicates_noteval 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 42 list result= l; 43 int j; int k; int i; 44 int deleted = 0; 45 int is_equal; 46 for (i = 1; i<= size(l); i++) 47 {//Iterate over the different factorizations 48 for (j = i+1; j<= size(l); j++) 49 {//Compare the i'th factorization to the j'th 50 if (size(l[i])!= size(l[j])) 51 {//different sizes => not equal 52 j++; 53 continue; 54 }//different sizes => not equal 55 is_equal = 1; 56 for (k = 1; k <= size(l[i]);k++) 57 {//Compare every entry 58 if (l[i][k]!=l[j][k]) 59 { 60 is_equal = 0; 61 break; 62 } 63 }//Compare every entry 64 if (is_equal == 1) 65 {//Delete this entry, because there is another equal one int the list 66 result = delete(result, i-deleted); 67 deleted = deleted+1; 68 break; 69 }//Delete this entry, because there is another equal one int the list 70 }//Compare the i'th factorization to the j'th 71 }//Iterate over the different factorizations 72 return(result); 60 73 }//proc delete_dublicates_noteval 61 74 … … 65 78 static proc delete_dublicates_eval(list l) 66 79 {//proc delete_dublicates_eval 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 80 list result=l; 81 int j; int k; int i; 82 int deleted = 0; 83 int is_equal; 84 for (i = 1; i<= size(result); i++) 85 {//Iterating over all elements in result 86 for (j = i+1; j<= size(result); j++) 87 {//comparing with the other elements 88 if (product(result[i]) == product(result[j])) 89 {//There are two equal results; throw away that one with the smaller size 90 if (size(result[i])>=size(result[j])) 91 {//result[i] has more entries 92 result = delete(result,j); 93 continue; 94 }//result[i] has more entries 95 else 96 {//result[j] has more entries 97 result = delete(result,i); 98 i--; 99 break; 100 }//result[j] has more entries 101 }//There are two equal results; throw away that one with the smaller size 102 }//comparing with the other elements 103 }//Iterating over all elements in result 104 return(result); 92 105 }//proc delete_dublicates_eval 93 106 … … 99 112 static proc combinekfinlf(list g, int nof, intvec limits) //nof stands for "number of factors" 100 113 {//Procedure combinekfinlf 101 list result; 102 int i; int j; int k; //iteration variables 103 list fc; //fc stands for "factors combined" 104 list temp; //a temporary store for factors 105 def nofgl = size(g); //nofgl stands for "number of factors of the given list" 106 if (nofgl == 0) 107 {//g was the empty list 108 return(result); 109 }//g was the empty list 110 if (nof <= 0) 111 {//The user wants to recieve a negative number or no element as a result 112 return(result); 113 }//The user wants to recieve a negative number or no element as a result 114 if (nofgl == nof) 115 {//There are no factors to combine 116 if (limitcheck(g,limits)) 117 { 118 result = result + list(g); 119 } 120 return(result); 121 }//There are no factors to combine 122 if (nof == 1) 123 {//User wants to get just one factor 124 if (limitcheck(list(product(g)),limits)) 125 { 126 result = result + list(list(product(g))); 127 } 128 return(result); 129 }//User wants to get just one factor 130 for (i = nof; i > 1; i--) 131 {//computing the possibilities that have at least one original factor from g 132 for (j = i; j>=1; j--) 133 {//shifting the window of combinable factors to the left 134 fc = combinekfinlf(list(g[(j)..(j+nofgl - i)]),nof - i + 1,limits); //fc stands for "factors combined" 135 for (k = 1; k<=size(fc); k++) 136 {//iterating over the different solutions of the smaller problem 137 if (j>1) 138 {//There are g_i before the combination 139 if (j+nofgl -i < nofgl) 140 {//There are g_i after the combination 141 temp = list(g[1..(j-1)]) + fc[k] + list(g[(j+nofgl-i+1)..nofgl]); 142 }//There are g_i after the combination 143 else 144 {//There are no g_i after the combination 145 temp = list(g[1..(j-1)]) + fc[k]; 146 }//There are no g_i after the combination 147 }//There are g_i before the combination 148 if (j==1) 149 {//There are no g_i before the combination 150 if (j+ nofgl -i <nofgl) 151 {//There are g_i after the combination 152 temp = fc[k]+ list(g[(j + nofgl - i +1)..nofgl]); 153 }//There are g_i after the combination 154 }//There are no g_i before the combination 155 if (limitcheck(temp,limits)) 156 { 157 result = result + list(temp); 158 } 159 }//iterating over the different solutions of the smaller problem 160 }//shifting the window of combinable factors to the left 161 }//computing the possibilities that have at least one original factor from g 162 for (i = 2; i<=nofgl/nof;i++) 163 {//getting the other possible results 164 result = result + combinekfinlf(list(product(list(g[1..i])))+list(g[(i+1)..nofgl]),nof,limits); 165 }//getting the other possible results 166 result = delete_dublicates_noteval(result); 167 return(result); 114 list result; 115 int i; int j; int k; //iteration variables 116 list fc; //fc stands for "factors combined" 117 list temp; //a temporary store for factors 118 def nofgl = size(g); //nofgl stands for "number of factors of the given list" 119 if (nofgl == 0) 120 {//g was the empty list 121 return(result); 122 }//g was the empty list 123 if (nof <= 0) 124 {//The user wants to recieve a negative number or no element as a result 125 return(result); 126 }//The user wants to recieve a negative number or no element as a result 127 if (nofgl == nof) 128 {//There are no factors to combine 129 if (limitcheck(g,limits)) 130 { 131 result = result + list(g); 132 } 133 return(result); 134 }//There are no factors to combine 135 if (nof == 1) 136 {//User wants to get just one factor 137 if (limitcheck(list(product(g)),limits)) 138 { 139 result = result + list(list(product(g))); 140 } 141 return(result); 142 }//User wants to get just one factor 143 for (i = nof; i > 1; i--) 144 {//computing the possibilities that have at least one original factor from g 145 for (j = i; j>=1; j--) 146 {//shifting the window of combinable factors to the left 147 //fc below stands for "factors combined" 148 fc = combinekfinlf(list(g[(j)..(j+nofgl - i)]),nof - i + 1,limits); 149 for (k = 1; k<=size(fc); k++) 150 {//iterating over the different solutions of the smaller problem 151 if (j>1) 152 {//There are g_i before the combination 153 if (j+nofgl -i < nofgl) 154 {//There are g_i after the combination 155 temp = list(g[1..(j-1)]) + fc[k] + list(g[(j+nofgl-i+1)..nofgl]); 156 }//There are g_i after the combination 157 else 158 {//There are no g_i after the combination 159 temp = list(g[1..(j-1)]) + fc[k]; 160 }//There are no g_i after the combination 161 }//There are g_i before the combination 162 if (j==1) 163 {//There are no g_i before the combination 164 if (j+ nofgl -i <nofgl) 165 {//There are g_i after the combination 166 temp = fc[k]+ list(g[(j + nofgl - i +1)..nofgl]); 167 }//There are g_i after the combination 168 }//There are no g_i before the combination 169 if (limitcheck(temp,limits)) 170 { 171 result = result + list(temp); 172 } 173 }//iterating over the different solutions of the smaller problem 174 }//shifting the window of combinable factors to the left 175 }//computing the possibilities that have at least one original factor from g 176 for (i = 2; i<=nofgl/nof;i++) 177 {//getting the other possible results 178 result = result + combinekfinlf(list(product(list(g[1..i])))+list(g[(i+1)..nofgl]),nof,limits); 179 }//getting the other possible results 180 result = delete_dublicates_noteval(result); 181 return(result); 168 182 }//Procedure combinekfinlf 169 183 … … 174 188 static proc merge_icf(list l1, list l2, intvec limits) 175 189 {//proc merge_icf 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 190 list g; 191 list f; 192 int i; int j; 193 if (size(l1)==0) 194 { 195 return(list()); 196 } 197 if (size(l2)==0) 198 { 199 return(list()); 200 } 201 if (size(l2)<=size(l1)) 202 {//l1 will be our g, l2 our f 203 g = l1; 204 f = l2; 205 }//l1 will be our g, l2 our f 206 else 207 {//l1 will be our f, l2 our g 208 g = l2; 209 f = l1; 210 }//l1 will be our f, l2 our g 211 def result = combinekfinlf(g,size(f),limits); 212 for (i = 1 ; i<= size(result); i++) 213 {//Adding the factors of f to every possibility listed in temp 214 for (j = 1; j<= size(f); j++) 215 { 216 result[i][j] = result[i][j]+f[j]; 217 } 218 if(!limitcheck(result[i],limits)) 219 { 220 result = delete(result,i); 221 i--; 222 } 223 }//Adding the factors of f to every possibility listed in temp 224 return(result); 211 225 }//proc merge_icf 212 226 … … 216 230 static proc merge_cf(list l1, list l2, intvec limits) 217 231 {//proc merge_cf 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 232 list g; 233 list f; 234 int i; int j; 235 list pre; 236 list post; 237 list candidate; 238 list temp; 239 int temppos; 240 if (size(l1)==0) 241 {//the first list is empty 242 return(list()); 243 }//the first list is empty 244 if(size(l2)==0) 245 {//the second list is empty 246 return(list()); 247 }//the second list is empty 248 if (size(l2)<=size(l1)) 249 {//l1 will be our g, l2 our f 250 g = l1; 251 f = l2; 252 }//l1 will be our g, l2 our f 253 else 254 {//l1 will be our f, l2 our g 255 g = l2; 256 f = l1; 257 }//l1 will be our f, l2 our g 258 list M; 259 for (i = 2; i<size(f); i++) 260 {//finding common factors of f and g... 247 261 for (j=2; j<size(g);j++) 248 262 {//... with g 249 250 251 252 263 if (f[i] == g[j]) 264 {//we have an equal pair 265 M = M + list(list(i,j)); 266 }//we have an equal pair 253 267 }//... with g 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 268 }//finding common factors of f and g... 269 if (g[1]==f[1]) 270 {//Checking for the first elements to be equal 271 M = M + list(list(1,1)); 272 }//Checking for the first elements to be equal 273 if (g[size(g)]==f[size(f)]) 274 {//Checking for the last elements to be equal 275 M = M + list(list(size(f),size(g))); 276 }//Checking for the last elements to be equal 277 list result;//= list(list()); 278 while(size(M)>0) 279 {//set of equal pairs is not empty 280 temp = M[1]; 281 temppos = 1; 282 for (i = 2; i<=size(M); i++) 283 {//finding the minimal element of M 284 if (M[i][1]<=temp[1]) 285 {//a possible candidate that is smaller than temp could have been found 286 if (M[i][1]==temp[1]) 287 {//In this case we must look at the second number 288 if (M[i][2]< temp[2]) 289 {//the candidate is smaller 290 temp = M[i]; 291 temppos = i; 292 }//the candidate is smaller 293 }//In this case we must look at the second number 294 else 295 {//The candidate is definately smaller 296 temp = M[i]; 297 temppos = i; 298 }//The candidate is definately smaller 299 }//a possible candidate that is smaller than temp could have been found 300 }//finding the minimal element of M 301 M = delete(M, temppos); 302 if(temp[1]>1) 303 {//There are factors to combine before the equal factor 304 if (temp[1]<size(f)) 305 {//The most common case 306 //first the combinations ignoring common factors 307 pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits); 308 post = merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits); 309 for (i = 1; i <= size(pre); i++) 310 {//all possible pre's... 311 for (j = 1; j<= size(post); j++) 312 {//...combined with all possible post's 313 candidate = pre[i]+list(f[temp[1]])+post[j]; 314 if (limitcheck(candidate,limits)) 315 { 316 result = result + list(candidate); 317 } 318 }//...combined with all possible post's 319 }//all possible pre's... 320 //Now the combinations with respect to common factors 321 post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits); 322 if (size(post)>0) 323 {//There are factors to combine 324 for (i = 1; i <= size(pre); i++) 325 {//all possible pre's... 326 for (j = 1; j<= size(post); j++) 327 {//...combined with all possible post's 328 candidate= pre[i]+list(f[temp[1]])+post[j]; 329 if (limitcheck(candidate,limits)) 330 { 331 result = result + list(candidate); 332 } 333 }//...combined with all possible post's 334 }//all possible pre's... 335 }//There are factors to combine 336 }//The most common case 337 else 338 {//the last factor is the common one 339 pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits); 340 for (i = 1; i<= size(pre); i++) 341 {//iterating over the possible pre-factors 342 candidate = pre[i]+list(f[temp[1]]); 343 if (limitcheck(candidate,limits)) 344 { 345 result = result + list(candidate); 346 } 347 }//iterating over the possible pre-factors 348 }//the last factor is the common one 349 }//There are factors to combine before the equal factor 350 else 351 {//There are no factors to combine before the equal factor 352 if (temp[1]<size(f)) 353 {//Just a check for security 354 //first without common factors 355 post=merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits); 356 for (i = 1; i<=size(post); i++) 357 { 358 candidate = list(f[temp[1]])+post[i]; 359 if (limitcheck(candidate,limits)) 360 { 361 result = result + list(candidate); 362 } 363 } 364 //Now with common factors 365 post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits); 366 if(size(post)>0) 367 {//we could find other combinations 368 for (i = 1; i<=size(post); i++) 369 { 370 candidate = list(f[temp[1]])+post[i]; 371 if (limitcheck(candidate,limits)) 372 { 373 result = result + list(candidate); 374 } 375 } 376 }//we could find other combinations 377 }//Just a check for security 378 }//There are no factors to combine before the equal factor 379 }//set of equal pairs is not empty 380 return(result); 367 381 }//proc merge_cf 368 382 … … 373 387 static proc mergence(list l1, list l2, intvec limits) 374 388 {//Procedure mergence 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 389 list g; 390 list f; 391 int l; int k; 392 list F; 393 if (size(l2)<=size(l1)) 394 {//l1 will be our g, l2 our f 395 g = l1; 396 f = l2; 397 }//l1 will be our g, l2 our f 398 else 399 {//l1 will be our f, l2 our g 400 g = l2; 401 f = l1; 402 }//l1 will be our f, l2 our g 403 list result; 404 for (l = size(f); l>=1; l--) 405 {//all possibilities to combine the factors of f 406 F = combinekfinlf(f,l,limits); 407 for (k = 1; k<= size(F);k++) 408 {//for all possibilities of combinations of the factors of f 409 result = result + merge_cf(F[k],g,limits); 410 result = result + merge_icf(F[k],g,limits); 411 }//for all possibilities of combinations of the factors of f 412 }//all possibilities to combine the factors of f 413 return(result); 400 414 }//Procedure mergence 401 415 … … 405 419 static proc limitcheck(list g, intvec limits) 406 420 {//proc limitcheck 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 421 int i; 422 if (size(limits)!=3) 423 {//check the input 424 return(0); 425 }//check the input 426 if(size(g)==0) 427 { 428 return(0); 429 } 430 def prod = product(g); 431 def limg = intvec(deg(prod,intvec(1,1)) ,deg(prod,intvec(1,0)),deg(prod,intvec(0,1))); 432 for (i = 1; i<=size(limg);i++) 433 {//the final check 434 if(limg[i]>limits[i]) 435 { 436 return(0); 437 } 438 }//the final check 439 return(1); 426 440 }//proc limitcheck 427 441 … … 430 444 //one factorization of a homogeneous polynomial 431 445 //in the first Weyl Algebra 432 static proc homogfac (poly h)433 "USAGE: homogfac (h); h is a homogeneous polynomial in the446 static proc homogfacFirstWeyl(poly h) 447 "USAGE: homogfacFirstWeyl(h); h is a homogeneous polynomial in the 434 448 first Weyl algebra with respect to the weight vector [-1,1] 435 449 RETURN: list 436 450 PURPOSE: Computes a factorization of a homogeneous polynomial h with 437 451 respect to the weight vector [-1,1] in the first Weyl algebra 438 THEORY: @code{homogfac } returns a list with a factorization of the given,452 THEORY: @code{homogfacFirstWeyl} returns a list with a factorization of the given, 439 453 [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with 440 454 k positive, the last k entries in the output list are the second 441 455 variable. If k is positive, the last k entries will be x. The other 442 456 entries will be irreducible polynomials of degree zero or 1 resp. -1. 443 SEE ALSO: homogfac _all444 "{//proc homogfac 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 //beginning to transform x^i*y^i in teta(teta-1)...(teta-i+1)482 483 484 485 486 487 ring tempRing = 0,(x,y,teta),dp;488 489 map tetamap = r,x,y;490 list mons = tetamap(mons);491 492 493 494 495 496 497 entry = entry * (teta-j);498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 }//proc homogfac 457 SEE ALSO: homogfacFirstWeyl_all 458 "{//proc homogfacFirstWeyl 459 def r = basering; 460 poly hath; 461 int i; int j; 462 if (!homogwithorder(h,intvec(-1,1))) 463 {//The given polynomial is not homogeneous 464 return(list()); 465 }//The given polynomial is not homogeneous 466 if (h==0) 467 { 468 return(list(0)); 469 } 470 list result; 471 int m = deg(h,intvec(-1,1)); 472 if (m!=0) 473 {//The degree is not zero 474 if (m <0) 475 {//There are more x than y 476 hath = lift(var(1)^(-m),h)[1,1]; 477 for (i = 1; i<=-m; i++) 478 { 479 result = result + list(var(1)); 480 } 481 }//There are more x than y 482 else 483 {//There are more y than x 484 hath = lift(var(2)^m,h)[1,1]; 485 for (i = 1; i<=m;i++) 486 { 487 result = result + list(var(2)); 488 } 489 }//There are more y than x 490 }//The degree is not zero 491 else 492 {//The degree is zero 493 hath = h; 494 }//The degree is zero 495 //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1) 496 list mons; 497 for(i = 1; i<=size(hath);i++) 498 {//Putting the monomials in a list 499 mons = mons+list(hath[i]); 500 }//Putting the monomials in a list 501 ring tempRing = 0,(x,y,theta),dp; 502 setring tempRing; 503 map thetamap = r,x,y; 504 list mons = thetamap(mons); 505 poly entry; 506 for (i = 1; i<=size(mons);i++) 507 {//transforming the monomials as monomials in theta 508 entry = leadcoef(mons[i]); 509 for (j = 0; j<leadexp(mons[i])[2];j++) 510 { 511 entry = entry * (theta-j); 512 } 513 mons[i] = entry; 514 }//transforming the monomials as monomials in theta 515 list azeroresult = factorize(sum(mons)); 516 list azeroresult_return_form; 517 for (i = 1; i<=size(azeroresult[1]);i++) 518 {//rewrite the result of the commutative factorization 519 for (j = 1; j <= azeroresult[2][i];j++) 520 { 521 azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]); 522 } 523 }//rewrite the result of the commutative factorization 524 setring(r); 525 map finalmap = tempRing,var(1),var(2),var(1)*var(2); 526 list tempresult = finalmap(azeroresult_return_form); 527 for (i = 1; i<=size(tempresult);i++) 528 {//factorizations of theta resp. theta +1 529 if(tempresult[i]==var(1)*var(2)) 530 { 531 tempresult = insert(tempresult,var(1),i-1); 532 i++; 533 tempresult[i]=var(2); 534 } 535 if(tempresult[i]==var(2)*var(1)) 536 { 537 tempresult = insert(tempresult,var(2),i-1); 538 i++; 539 tempresult[i]=var(1); 540 } 541 }//factorizations of theta resp. theta +1 542 result = tempresult+result; 543 return(result); 544 }//proc homogfacFirstWeyl 531 545 /* example */ 532 546 /* { */ … … 536 550 /* setring(r); */ 537 551 /* poly h = (x^2*y^2+1)*(x^4); */ 538 /* homogfac (h); */552 /* homogfacFirstWeyl(h); */ 539 553 /* } */ 540 554 541 555 //================================================== 542 556 //Computes all possible homogeneous factorizations 543 static proc homogfac _all(poly h)544 "USAGE: homogfac _all(h); h is a homogeneous polynomial in the first Weyl algebra557 static proc homogfacFirstWeyl_all(poly h) 558 "USAGE: homogfacFirstWeyl_all(h); h is a homogeneous polynomial in the first Weyl algebra 545 559 with respect to the weight vector [-1,1] 546 560 RETURN: list 547 561 PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect 548 562 to the weight vector [-1,1] in the first Weyl algebra 549 THEORY: @code{homogfac } returns a list with all factorization of the given,550 homogeneous polynomial. It uses the output of homogfac and permutes563 THEORY: @code{homogfacFirstWeyl} returns a list with all factorization of the given, 564 homogeneous polynomial. It uses the output of homogfacFirstWeyl and permutes 551 565 its entries with respect to the commutation rule. Furthermore, if a 552 566 factor of degree zero is irreducible in K[\theta], but reducible in 553 567 the first Weyl algebra, the permutations of this element with the other 554 568 entries will also be computed. 555 SEE ALSO: homogfac 556 "{//proc HomogFacAll 557 if (deg(h,intvec(1,1)) <= 0 ) 558 {//h is a constant 559 return(list(list(h))); 560 }//h is a constant 561 def r = basering; 562 list one_hom_fac; //stands for one homogeneous factorization 563 int i; int j; int k; 564 //Compute again a homogeneous factorization 565 one_hom_fac = homogfac(h); 566 if (size(one_hom_fac) == 0) 567 {//there is no homogeneous factorization or the polynomial was not homogeneous 568 return(list()); 569 }//there is no homogeneous factorization or the polynomial was not homogeneous 570 //divide list in A0-Part and a list of x's resp. y's 571 list list_not_azero = list(); 572 list list_azero; 573 list k_factor; 574 int is_list_not_azero_empty = 1; 575 int is_list_azero_empty = 1; 576 k_factor = list(one_hom_fac[1]); 577 if (absValue(deg(h,intvec(-1,1)))<size(one_hom_fac)-1) 578 {//There is a nontrivial A_0-part 579 list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,intvec(-1,1))))]; 580 is_list_azero_empty = 0; 581 }//There is a nontrivial A_0 part 582 for (i = 1; i<=size(list_azero)-1;i++) 583 {//in homogfac, we factorized theta, and this will be made undone 584 if (list_azero[i] == var(1)) 585 { 586 if (list_azero[i+1]==var(2)) 587 { 588 list_azero[i] = var(1)*var(2); 589 list_azero = delete(list_azero,i+1); 590 } 591 } 592 if (list_azero[i] == var(2)) 593 { 594 if (list_azero[i+1]==var(1)) 595 { 596 list_azero[i] = var(2)*var(1); 597 list_azero = delete(list_azero,i+1); 598 } 599 } 600 }//in homogfac, we factorized theta, and this will be made undone 601 if(deg(h,intvec(-1,1))!=0) 602 {//list_not_azero is not empty 603 list_not_azero = one_hom_fac[(size(one_hom_fac)-absValue(deg(h,intvec(-1,1)))+1)..size(one_hom_fac)]; 604 is_list_not_azero_empty = 0; 605 }//list_not_azero is not empty 606 //Map list_azero in K[theta] 607 ring tempRing = 0,(x,y,theta), dp; 608 setring(tempRing); 609 poly entry; 610 map thetamap = r,x,y; 611 if(!is_list_not_azero_empty) 612 {//Mapping in Singular is only possible, if the list before contained at least one element of the other ring 613 list list_not_azero = thetamap(list_not_azero); 614 }//Mapping in Singular is only possible, if the list before contained at least one element of the other ring 615 if(!is_list_azero_empty) 616 {//Mapping in Singular is only possible, if the list before contained at least one element of the other ring 617 list list_azero= thetamap(list_azero); 618 }//Mapping in Singular is only possible, if the list before contained at least one element of the other ring 619 list k_factor = thetamap(k_factor); 620 list tempmons; 621 for(i = 1; i<=size(list_azero);i++) 622 {//rewrite the polynomials in A1 as polynomials in K[theta] 623 tempmons = list(); 624 for (j = 1; j<=size(list_azero[i]);j++) 625 { 626 tempmons = tempmons + list(list_azero[i][j]); 627 } 628 for (j = 1 ; j<=size(tempmons);j++) 629 { 630 entry = leadcoef(tempmons[j]); 631 for (k = 0; k < leadexp(tempmons[j])[2];k++) 632 { 633 entry = entry*(theta-k); 634 } 635 tempmons[j] = entry; 636 } 637 list_azero[i] = sum(tempmons); 638 }//rewrite the polynomials in A1 as polynomials in K[theta] 639 //Compute all permutations of the A0-part 640 list result; 641 int shift_sign; 642 int shift; 643 poly shiftvar; 644 if (size(list_not_azero)!=0) 645 {//Compute all possibilities to permute the x's resp. the y's in the list 646 if (list_not_azero[1] == x) 647 {//h had a negative weighted degree 648 shift_sign = 1; 649 shiftvar = x; 650 }//h had a negative weighted degree 651 else 652 {//h had a positive weighted degree 653 shift_sign = -1; 654 shiftvar = y; 655 }//h had a positive weighted degree 656 result = permpp(list_azero + list_not_azero); 657 for (i = 1; i<= size(result); i++) 658 {//adjust the a_0-parts 659 shift = 0; 660 for (j=1; j<=size(result[i]);j++) 661 { 662 if (result[i][j]==shiftvar) 663 { 664 shift = shift + shift_sign; 665 } 666 else 667 { 668 result[i][j] = subst(result[i][j],theta,theta + shift); 669 } 670 } 671 }//adjust the a_0-parts 672 }//Compute all possibilities to permute the x's resp. the y's in the list 673 else 674 {//The result is just all the permutations of the a_0-part 675 result = permpp(list_azero); 676 }//The result is just all the permutations of the a_0 part 677 if (size(result)==0) 678 { 679 return(result); 680 } 681 //Now we are going deeper and search for theta resp. theta + 1, substitute them by xy resp. yx and go on permuting 682 int found_theta; 683 int thetapos; 684 list leftpart; 685 list rightpart; 686 list lparts; 687 list rparts; 688 list tempadd; 689 for (i = 1; i<=size(result) ; i++) 690 {//checking every entry of result for theta or theta +1 691 found_theta = 0; 692 for(j=1;j<=size(result[i]);j++) 693 { 694 if (result[i][j]==theta) 695 {//the jth entry is theta and can be written as x*y 696 thetapos = j; 697 result[i]= insert(result[i],x,j-1); 698 j++; 699 result[i][j] = y; 700 found_theta = 1; 701 break; 702 }//the jth entry is theta and can be written as x*y 703 if(result[i][j] == theta +1) 704 { 705 thetapos = j; 706 result[i] = insert(result[i],y,j-1); 707 j++; 708 result[i][j] = x; 709 found_theta = 1; 710 break; 711 } 712 } 713 if (found_theta) 714 {//One entry was theta resp. theta +1 715 leftpart = result[i]; 716 leftpart = leftpart[1..thetapos]; 717 rightpart = result[i]; 718 rightpart = rightpart[(thetapos+1)..size(rightpart)]; 719 lparts = list(leftpart); 720 rparts = list(rightpart); 721 //first deal with the left part 722 if (leftpart[thetapos] == x) 723 { 724 shift_sign = 1; 725 shiftvar = x; 726 } 727 else 728 { 729 shift_sign = -1; 730 shiftvar = y; 731 } 732 for (j = size(leftpart); j>1;j--) 733 {//drip x resp. y 734 if (leftpart[j-1]==shiftvar) 735 {//commutative 736 j--; 737 continue; 738 }//commutative 739 if (deg(leftpart[j-1],intvec(-1,1,0))!=0) 740 {//stop here 741 break; 742 }//stop here 743 //Here, we can only have a a0- part 744 leftpart[j] = subst(leftpart[j-1],theta, theta + shift_sign); 745 leftpart[j-1] = shiftvar; 746 lparts = lparts + list(leftpart); 747 }//drip x resp. y 748 //and now deal with the right part 749 if (rightpart[1] == x) 750 { 751 shift_sign = 1; 752 shiftvar = x; 753 } 754 else 755 { 756 shift_sign = -1; 757 shiftvar = y; 758 } 759 for (j = 1 ; j < size(rightpart); j++) 760 { 761 if (rightpart[j+1] == shiftvar) 762 { 763 j++; 764 continue; 765 } 766 if (deg(rightpart[j+1],intvec(-1,1,0))!=0) 767 { 768 break; 769 } 770 rightpart[j] = subst(rightpart[j+1], theta, theta - shift_sign); 771 rightpart[j+1] = shiftvar; 772 rparts = rparts + list(rightpart); 773 } 774 //And now, we put all possibilities together 775 tempadd = list(); 776 for (j = 1; j<=size(lparts); j++) 777 { 778 for (k = 1; k<=size(rparts);k++) 779 { 780 tempadd = tempadd + list(lparts[j]+rparts[k]); 781 } 782 } 783 tempadd = delete(tempadd,1); // The first entry is already in the list 784 result = result + tempadd; 785 continue; //We can may be not be done already with the ith entry 786 }//One entry was theta resp. theta +1 787 }//checking every entry of result for theta or theta +1 788 //map back to the basering 789 setring(r); 790 map finalmap = tempRing, var(1), var(2),var(1)*var(2); 791 list result = finalmap(result); 792 for (i=1; i<=size(result);i++) 793 {//adding the K factor 794 result[i] = k_factor + result[i]; 795 }//adding the k-factor 796 result = delete_dublicates_noteval(result); 797 return(result); 798 }//proc HomogFacAll 569 SEE ALSO: homogfacFirstWeyl 570 "{//proc HomogfacFirstWeylAll 571 if (deg(h,intvec(1,1)) <= 0 ) 572 {//h is a constant 573 return(list(list(h))); 574 }//h is a constant 575 def r = basering; 576 list one_hom_fac; //stands for one homogeneous factorization 577 int i; int j; int k; 578 //Compute again a homogeneous factorization 579 one_hom_fac = homogfacFirstWeyl(h); 580 if (size(one_hom_fac) == 0) 581 {//there is no homogeneous factorization or the polynomial was not homogeneous 582 return(list()); 583 }//there is no homogeneous factorization or the polynomial was not homogeneous 584 //divide list in A0-Part and a list of x's resp. y's 585 list list_not_azero = list(); 586 list list_azero; 587 list k_factor; 588 int is_list_not_azero_empty = 1; 589 int is_list_azero_empty = 1; 590 k_factor = list(one_hom_fac[1]); 591 if (absValue(deg(h,intvec(-1,1)))<size(one_hom_fac)-1) 592 {//There is a nontrivial A_0-part 593 list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,intvec(-1,1))))]; 594 is_list_azero_empty = 0; 595 }//There is a nontrivial A_0 part 596 for (i = 1; i<=size(list_azero)-1;i++) 597 {//in homogfacFirstWeyl, we factorized theta, and this will be made undone 598 if (list_azero[i] == var(1)) 599 { 600 if (list_azero[i+1]==var(2)) 601 { 602 list_azero[i] = var(1)*var(2); 603 list_azero = delete(list_azero,i+1); 604 } 605 } 606 if (list_azero[i] == var(2)) 607 { 608 if (list_azero[i+1]==var(1)) 609 { 610 list_azero[i] = var(2)*var(1); 611 list_azero = delete(list_azero,i+1); 612 } 613 } 614 }//in homogfacFirstWeyl, we factorized theta, and this will be made undone 615 if(deg(h,intvec(-1,1))!=0) 616 {//list_not_azero is not empty 617 list_not_azero = 618 one_hom_fac[(size(one_hom_fac)-absValue(deg(h,intvec(-1,1)))+1)..size(one_hom_fac)]; 619 is_list_not_azero_empty = 0; 620 }//list_not_azero is not empty 621 //Map list_azero in K[theta] 622 ring tempRing = 0,(x,y,theta), dp; 623 setring(tempRing); 624 poly entry; 625 map thetamap = r,x,y; 626 if(!is_list_not_azero_empty) 627 {//Mapping in Singular is only possible, if the list before 628 //contained at least one element of the other ring 629 list list_not_azero = thetamap(list_not_azero); 630 }//Mapping in Singular is only possible, if the list before 631 //contained at least one element of the other ring 632 if(!is_list_azero_empty) 633 {//Mapping in Singular is only possible, if the list before 634 //contained at least one element of the other ring 635 list list_azero= thetamap(list_azero); 636 }//Mapping in Singular is only possible, if the list before 637 //contained at least one element of the other ring 638 list k_factor = thetamap(k_factor); 639 list tempmons; 640 for(i = 1; i<=size(list_azero);i++) 641 {//rewrite the polynomials in A1 as polynomials in K[theta] 642 tempmons = list(); 643 for (j = 1; j<=size(list_azero[i]);j++) 644 { 645 tempmons = tempmons + list(list_azero[i][j]); 646 } 647 for (j = 1 ; j<=size(tempmons);j++) 648 { 649 entry = leadcoef(tempmons[j]); 650 for (k = 0; k < leadexp(tempmons[j])[2];k++) 651 { 652 entry = entry*(theta-k); 653 } 654 tempmons[j] = entry; 655 } 656 list_azero[i] = sum(tempmons); 657 }//rewrite the polynomials in A1 as polynomials in K[theta] 658 //Compute all permutations of the A0-part 659 list result; 660 int shift_sign; 661 int shift; 662 poly shiftvar; 663 if (size(list_not_azero)!=0) 664 {//Compute all possibilities to permute the x's resp. the y's in the list 665 if (list_not_azero[1] == x) 666 {//h had a negative weighted degree 667 shift_sign = 1; 668 shiftvar = x; 669 }//h had a negative weighted degree 670 else 671 {//h had a positive weighted degree 672 shift_sign = -1; 673 shiftvar = y; 674 }//h had a positive weighted degree 675 result = permpp(list_azero + list_not_azero); 676 for (i = 1; i<= size(result); i++) 677 {//adjust the a_0-parts 678 shift = 0; 679 for (j=1; j<=size(result[i]);j++) 680 { 681 if (result[i][j]==shiftvar) 682 { 683 shift = shift + shift_sign; 684 } 685 else 686 { 687 result[i][j] = subst(result[i][j],theta,theta + shift); 688 } 689 } 690 }//adjust the a_0-parts 691 }//Compute all possibilities to permute the x's resp. the y's in the list 692 else 693 {//The result is just all the permutations of the a_0-part 694 result = permpp(list_azero); 695 }//The result is just all the permutations of the a_0 part 696 if (size(result)==0) 697 { 698 return(result); 699 } 700 //Now we are going deeper and search for theta resp. theta + 1, substitute 701 //them by xy resp. yx and go on permuting 702 int found_theta; 703 int thetapos; 704 list leftpart; 705 list rightpart; 706 list lparts; 707 list rparts; 708 list tempadd; 709 for (i = 1; i<=size(result) ; i++) 710 {//checking every entry of result for theta or theta +1 711 found_theta = 0; 712 for(j=1;j<=size(result[i]);j++) 713 { 714 if (result[i][j]==theta) 715 {//the jth entry is theta and can be written as x*y 716 thetapos = j; 717 result[i]= insert(result[i],x,j-1); 718 j++; 719 result[i][j] = y; 720 found_theta = 1; 721 break; 722 }//the jth entry is theta and can be written as x*y 723 if(result[i][j] == theta +1) 724 { 725 thetapos = j; 726 result[i] = insert(result[i],y,j-1); 727 j++; 728 result[i][j] = x; 729 found_theta = 1; 730 break; 731 } 732 } 733 if (found_theta) 734 {//One entry was theta resp. theta +1 735 leftpart = result[i]; 736 leftpart = leftpart[1..thetapos]; 737 rightpart = result[i]; 738 rightpart = rightpart[(thetapos+1)..size(rightpart)]; 739 lparts = list(leftpart); 740 rparts = list(rightpart); 741 //first deal with the left part 742 if (leftpart[thetapos] == x) 743 { 744 shift_sign = 1; 745 shiftvar = x; 746 } 747 else 748 { 749 shift_sign = -1; 750 shiftvar = y; 751 } 752 for (j = size(leftpart); j>1;j--) 753 {//drip x resp. y 754 if (leftpart[j-1]==shiftvar) 755 {//commutative 756 j--; 757 continue; 758 }//commutative 759 if (deg(leftpart[j-1],intvec(-1,1,0))!=0) 760 {//stop here 761 break; 762 }//stop here 763 //Here, we can only have a a0- part 764 leftpart[j] = subst(leftpart[j-1],theta, theta + shift_sign); 765 leftpart[j-1] = shiftvar; 766 lparts = lparts + list(leftpart); 767 }//drip x resp. y 768 //and now deal with the right part 769 if (rightpart[1] == x) 770 { 771 shift_sign = 1; 772 shiftvar = x; 773 } 774 else 775 { 776 shift_sign = -1; 777 shiftvar = y; 778 } 779 for (j = 1 ; j < size(rightpart); j++) 780 { 781 if (rightpart[j+1] == shiftvar) 782 { 783 j++; 784 continue; 785 } 786 if (deg(rightpart[j+1],intvec(-1,1,0))!=0) 787 { 788 break; 789 } 790 rightpart[j] = subst(rightpart[j+1], theta, theta - shift_sign); 791 rightpart[j+1] = shiftvar; 792 rparts = rparts + list(rightpart); 793 } 794 //And now, we put all possibilities together 795 tempadd = list(); 796 for (j = 1; j<=size(lparts); j++) 797 { 798 for (k = 1; k<=size(rparts);k++) 799 { 800 tempadd = tempadd + list(lparts[j]+rparts[k]); 801 } 802 } 803 tempadd = delete(tempadd,1); // The first entry is already in the list 804 result = result + tempadd; 805 continue; //We can may be not be done already with the ith entry 806 }//One entry was theta resp. theta +1 807 }//checking every entry of result for theta or theta +1 808 //map back to the basering 809 setring(r); 810 map finalmap = tempRing, var(1), var(2),var(1)*var(2); 811 list result = finalmap(result); 812 for (i=1; i<=size(result);i++) 813 {//adding the K factor 814 result[i] = k_factor + result[i]; 815 }//adding the k-factor 816 result = delete_dublicates_noteval(result); 817 return(result); 818 }//proc HomogfacFirstWeylAll 799 819 /* example */ 800 820 /* { */ … … 804 824 /* setring(r); */ 805 825 /* poly h = (x^2*y^2+1)*(x^4); */ 806 /* homogfac _all(h); */826 /* homogfacFirstWeyl_all(h); */ 807 827 /* } */ 808 828 … … 811 831 static proc perm(list l) 812 832 {//proc perm 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 833 int i; int j; 834 list tempresult; 835 list result; 836 if (size(l)==0) 837 { 838 return(list()); 839 } 840 if (size(l)==1) 841 { 842 return(list(l)); 843 } 844 for (i = 1; i<=size(l); i++ ) 845 { 846 tempresult = perm(delete(l,i)); 847 for (j = 1; j<=size(tempresult);j++) 848 { 849 tempresult[j] = list(l[i])+tempresult[j]; 850 } 851 result = result+tempresult; 852 } 853 return(result); 834 854 }//proc perm 835 855 … … 839 859 static proc permpp(list l) 840 860 {//proc permpp 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 861 int i; int j; 862 list tempresult; 863 list l_without_double; 864 list l_without_double_pos; 865 int double_entry; 866 list result; 867 if (size(l)==0) 868 { 869 return(list()); 870 } 871 if (size(l)==1) 872 { 873 return(list(l)); 874 } 875 for (i = 1; i<=size(l);i++) 876 {//Filling the list with unique entries 877 double_entry = 0; 878 for (j = 1; j<=size(l_without_double);j++) 879 { 880 if (l_without_double[j] == l[i]) 881 { 882 double_entry = 1; 883 break; 884 } 885 } 886 if (!double_entry) 887 { 888 l_without_double = l_without_double + list(l[i]); 889 l_without_double_pos = l_without_double_pos + list(i); 890 } 891 }//Filling the list with unique entries 892 for (i = 1; i<=size(l_without_double); i++ ) 893 { 894 tempresult = permpp(delete(l,l_without_double_pos[i])); 895 for (j = 1; j<=size(tempresult);j++) 896 { 897 tempresult[j] = list(l_without_double[i])+tempresult[j]; 898 } 899 result = result+tempresult; 900 } 901 return(result); 882 902 }//proc permpp 883 903 … … 890 910 //variables in a certain order. 891 911 proc facFirstWeyl(poly h) 892 "USAGE: facFirstWeyl(h); h a poly 912 "USAGE: facFirstWeyl(h); h a polynomial in the first Weyl algebra 893 913 RETURN: list 894 914 PURPOSE: compute all factorizations of a polynomial in the first Weyl algebra 895 915 THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle 896 916 ASSUME: basering in the first Weyl algebra 917 NOTE: Every entry of the output list is a list with factors for one possible factorization. 918 The first factor is always a constant (1, if no nontrivial constant could be excluded). 897 919 EXAMPLE: example facFirstWeyl; shows examples 898 SEE ALSO: facSubWeyl, testNCfac 920 SEE ALSO: facSubWeyl, testNCfac, facFirstShift 899 921 "{//proc facFirstWeyl 900 //Redefine the ring in my standard form 901 if (!isWeyl()) 902 {//Our basering is not the Weyl algebra 903 return(list()); 904 }//Our basering is not the Weyl algebra 905 if(nvars(basering)!=2) 906 {//Our basering is the Weyl algebra, but not the first 907 return(list()); 908 }//Our basering is the Weyl algebra, but not the first 909 if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why 910 { 911 def r = basering; 912 ring tempRing = ringlist(r)[1][1],(x,y),Ws(-1,1); 913 def tempRingnc = nc_algebra(1,1); 914 setring(tempRingnc); 915 map transf = r, var(2), var(1); 916 list result = sfacwa(transf(h)); 917 setring(r); 918 map transfback = tempRingnc, var(2),var(1); 919 return(transfback(result)); 920 } 921 else { return(sfacwa(h));} 922 //Redefine the ring in my standard form 923 if (!isWeyl()) 924 {//Our basering is not the Weyl algebra 925 return(list()); 926 }//Our basering is not the Weyl algebra 927 if(nvars(basering)!=2) 928 {//Our basering is the Weyl algebra, but not the first 929 return(list()); 930 }//Our basering is the Weyl algebra, but not the first 931 list result = list(); 932 int i;int j; int k; int l; //counter 933 if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why 934 { 935 def r = basering; 936 ring tempRing = ringlist(r)[1][1],(x,y),Ws(-1,1); // very strange: 937 // setting Wp(-1,1) leads to SegFault; to clarify why!!! 938 def NTR = nc_algebra(1,1); 939 setring NTR ; 940 map transf = r, var(2), var(1); 941 list resulttemp = sfacwa(transf(h)); 942 setring(r); 943 map transfback = NTR, var(2),var(1); 944 result = transfback(resulttemp); 945 } 946 else { result = sfacwa(h);} 947 list recursivetemp; 948 for(i = 1; i<=size(result);i++) 949 {//recursively factorize factors 950 if(size(result[i])>2) 951 {//Nontrivial factorization 952 for (j=2;j<=size(result[i]);j++) 953 {//Factorize every factor 954 recursivetemp = facFirstWeyl(result[i][j]); 955 if(size(recursivetemp)>1) 956 {//we have a nontrivial factorization 957 for(k=1; k<=size(recursivetemp);k++) 958 {//insert factorized factors 959 if(size(recursivetemp[k])>2) 960 {//nontrivial 961 result = insert(result,result[i],i); 962 for(l = size(recursivetemp[k]);l>=2;l--) 963 { 964 result[i+1] = insert(result[i+1],recursivetemp[k][l],j); 965 } 966 result[i+1] = delete(result[i+1],j); 967 }//nontrivial 968 }//insert factorized factors 969 }//we have a nontrivial factorization 970 }//Factorize every factor 971 }//Nontrivial factorization 972 }//recursively factorize factors 973 if (size(result)==0) 974 {//only the trivial factorization could be found 975 result = list(list(1,h)); 976 }//only the trivial factorization could be found 977 return(result); 922 978 }//proc facFirstWeyl 923 979 example 924 980 { 925 926 927 928 929 930 981 "EXAMPLE:";echo=2; 982 ring R = 0,(x,y),dp; 983 def r = nc_algebra(1,1); 984 setring(r); 985 poly h = (x^2*y^2+x)*(x+1); 986 facFirstWeyl(h); 931 987 } 932 988 … … 940 996 homogeneous part and those of the lowest will be merged. If during this 941 997 procedure a factorization of the polynomial occurs, it will be added to 942 the output list. For a more detail led description visit943 @url{http:// math.rwth-aachen.de/\~Albert.Heinle}944 SEE ALSO: homogfac _all, homogfac998 the output list. For a more detailed description visit 999 @url{http://www.math.rwth-aachen.de/\~Albert.Heinle} 1000 SEE ALSO: homogfacFirstWeyl_all, homogfacFirstWeyl 945 1001 "{//proc sfacwa 946 947 948 return(homogfac_all(h));949 950 951 952 953 954 955 956 957 958 959 960 961 962 list f1 = homogfac_all(maxh);963 list f2 = homogfac_all(minh);964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 f1 = homogfac_all(maxh);1050 f2 = homogfac_all(minh);1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1002 if(homogwithorder(h,intvec(-1,1))) 1003 { 1004 return(homogfacFirstWeyl_all(h)); 1005 } 1006 def r = basering; 1007 map invo = basering,-var(1),var(2); 1008 int i; int j; int k; 1009 intvec limits = deg(h,intvec(1,1)) ,deg(h,intvec(1,0)),deg(h,intvec(0,1)); 1010 def prod; 1011 //end finding the limits 1012 poly maxh = jet(h,deg(h,intvec(-1,1)),intvec(-1,1))-jet(h,deg(h,intvec(-1,1))-1,intvec(-1,1)); 1013 poly minh = jet(h,deg(h,intvec(1,-1)),intvec(1,-1))-jet(h,deg(h,intvec(1,-1))-1,intvec(1,-1)); 1014 list result; 1015 list temp; 1016 list homogtemp; 1017 list M; list hatM; 1018 list f1 = homogfacFirstWeyl_all(maxh); 1019 list f2 = homogfacFirstWeyl_all(minh); 1020 int is_equal; 1021 poly hath; 1022 for (i = 1; i<=size(f1);i++) 1023 {//Merge all combinations 1024 for (j = 1; j<=size(f2); j++) 1025 { 1026 M = M + mergence(f1[i],f2[j],limits); 1027 } 1028 }//Merge all combinations 1029 for (i = 1 ; i<= size(M); i++) 1030 {//filter valid combinations 1031 if (product(M[i]) == h) 1032 {//We have one factorization 1033 result = result + list(M[i]); 1034 M = delete(M,i); 1035 continue; 1036 }//We have one factorization 1037 else 1038 { 1039 if (deg(h,intvec(-1,1))<=deg(h-product(M[i]),intvec(-1,1))) 1040 { 1041 M = delete(M,i); 1042 continue; 1043 } 1044 if (deg(h,intvec(1,-1))<=deg(h-product(M[i]),intvec(1,-1))) 1045 { 1046 M = delete(M,i); 1047 continue; 1048 } 1049 } 1050 }//filter valid combinations 1051 M = delete_dublicates_eval(M); 1052 while(size(M)>0) 1053 {//work on the elements of M 1054 hatM = list(); 1055 for(i = 1; i<=size(M); i++) 1056 {//iterate over all elements of M 1057 hath = h-product(M[i]); 1058 temp = list(); 1059 //First check for common inhomogeneous factors between hath and h 1060 if (involution(NF(involution(hath,invo), std(involution(ideal(M[i][1]),invo))),invo)==0) 1061 {//hath and h have a common factor on the left 1062 j = 1; 1063 f1 = M[i]; 1064 if (j+1<=size(f1)) 1065 {//Checking for more than one common factor 1066 while(involution(NF(involution(hath,invo),std(involution(ideal(product(f1[1..(j+1)])),invo))),invo)==0) 1067 { 1068 if (j+1<size(f1)) 1069 { 1070 j++; 1071 } 1072 else 1073 { 1074 break; 1075 } 1076 } 1077 }//Checking for more than one common factor 1078 f2 = list(f1[1..j])+list(involution(lift(involution(product(f1[1..j]),invo),involution(hath,invo))[1,1],invo)); 1079 temp = temp + merge_cf(f2,f1,limits); 1080 }//hath and h have a common factor on the left 1081 if (reduce(hath, std(ideal(M[i][size(M[i])])))==0) 1082 {//hath and h have a common factor on the right 1083 j = size(M[i]); 1084 f1 = M[i]; 1085 if (j-1>0) 1086 {//Checking for more than one factor 1087 while(reduce(hath,std(ideal(product(f1[(j-1)..size(f1)]))))==0) 1088 { 1089 if (j-1>1) 1090 { 1091 j--; 1092 } 1093 else 1094 { 1095 break; 1096 } 1097 } 1098 }//Checking for more than one factor 1099 f2 = list(lift(product(f1[j..size(f1)]),hath)[1,1])+list(f1[j..size(f1)]); 1100 temp = temp + merge_cf(f2,M[i],limits); 1101 }//hath and h have a common factor on the right 1102 //and now the homogeneous 1103 maxh = jet(hath,deg(hath,intvec(-1,1)),intvec(-1,1))-jet(hath,deg(hath,intvec(-1,1))-1,intvec(-1,1)); 1104 minh = jet(hath,deg(hath,intvec(1,-1)),intvec(1,-1))-jet(hath,deg(hath,intvec(1,-1))-1,intvec(1,-1)); 1105 f1 = homogfacFirstWeyl_all(maxh); 1106 f2 = homogfacFirstWeyl_all(minh); 1107 for (j = 1; j<=size(f1);j++) 1108 { 1109 for (k=1; k<=size(f2);k++) 1110 { 1111 homogtemp = mergence(f1[j],f2[k],limits); 1112 } 1113 } 1114 for (j = 1; j<= size(homogtemp); j++) 1115 { 1116 temp = temp + mergence(homogtemp[j],M[i],limits); 1117 } 1118 for (j = 1; j<=size(temp); j++) 1119 {//filtering invalid entries in temp 1120 if(product(temp[j])==h) 1121 {//This is already a result 1122 result = result + list(temp[j]); 1123 temp = delete(temp,j); 1124 continue; 1125 }//This is already a result 1126 if (deg(hath,intvec(-1,1))<=deg(hath-product(temp[j]),intvec(-1,1))) 1127 { 1128 temp = delete(temp,j); 1129 continue; 1130 } 1131 }//filtering invalid entries in temp 1132 hatM = hatM + temp; 1133 }//iterate over all elements of M 1134 M = hatM; 1135 for (i = 1; i<=size(M);i++) 1136 {//checking for complete factorizations 1137 if (h == product(M[i])) 1138 { 1139 result = result + list(M[i]); 1140 M = delete(M,i); 1141 continue; 1142 } 1143 }//checking for complete factorizations 1144 M = delete_dublicates_eval(M); 1145 }//work on the elements of M 1146 //In the case, that there is none, write a constant factor before the factor of interest. 1147 for (i = 1 ; i<=size(result);i++) 1148 {//add a constant factor 1149 if (deg(result[i][1],intvec(1,1))!=0) 1150 { 1151 result[i] = insert(result[i],1); 1152 } 1153 }//add a constant factor 1154 result = delete_dublicates_noteval(result); 1155 return(result); 1100 1156 }//proc sfacwa 1101 1157 … … 1103 1159 //================================================== 1104 1160 /*Singular has no way implemented to test polynomials 1105 for homogenity with respect to a weight vector.1106 The following procedure does exactly this*/1161 for homogenity with respect to a weight vector. 1162 The following procedure does exactly this*/ 1107 1163 static proc homogwithorder(poly h, intvec weights) 1108 1164 {//proc homogwithorder 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1165 if(size(weights) != nvars(basering)) 1166 {//The user does not know how many variables the current ring has 1167 return(0); 1168 }//The user does not know how many variables the current ring has 1169 int i; 1170 int dofp = deg(h,weights); //degree of polynomial 1171 for (i = 1; i<=size(h);i++) 1172 { 1173 if (deg(h[i],weights)!=dofp) 1174 { 1175 return(0); 1176 } 1177 } 1178 return(1); 1123 1179 }//proc homogwithorder 1124 1180 … … 1154 1210 by the multiplied candidate for its factorization. 1155 1211 EXAMPLE: example testNCfac; shows examples 1156 SEE ALSO: facFirstWeyl, facSubWeyl 1212 SEE ALSO: facFirstWeyl, facSubWeyl, facFirstShift 1157 1213 "{//proc testfac 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1214 if (size(l)==0) 1215 {//The empty list is given 1216 return(list()); 1217 }//The empty list is given 1218 if (size(#)>1) 1219 {//We want max. one optional argument 1220 return(list()); 1221 }//We want max. one optional argument 1222 list result; 1223 int i; int j; 1224 if (size(#)==0) 1225 {//No optional argument is given 1226 int valid = 1; 1227 for (i = size(l);i>=1;i--) 1228 {//iterate over the elements of the given list 1229 if (size(result)>0) 1230 { 1231 if (product(l[i])!=result[size(l)-i]) 1232 { 1233 valid = 0; 1234 break; 1235 } 1236 } 1237 result = insert(result, product(l[i])); 1238 }//iterate over the elements of the given list 1239 if (valid) 1240 { 1241 return(result); 1242 } 1243 return(list()); 1244 }//No optional argument is given 1245 else 1246 { 1247 int valid = 1; 1248 for (i = size(l);i>=1;i--) 1249 {//iterate over the elements of the given list 1250 if (product(l[i])!=#[1]) 1251 { 1252 valid = 0; 1253 } 1254 result = insert(result, product(l[i])); 1255 }//iterate over the elements of the given list 1256 if (valid) 1257 { 1258 return(result); 1259 } 1260 for (i = 1 ; i<= size(result);i++) 1261 {//subtract p from every entry in result 1262 result[i] = result[i] - #[1]; 1263 }//subtract p from every entry in result 1264 return(result); 1265 } 1210 1266 }//proc testfac 1211 1267 example 1212 1268 { 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1269 "EXAMPLE:";echo=2; 1270 ring r = 0,(x,y),dp; 1271 def R = nc_algebra(1,1); 1272 setring R; 1273 poly h = (x^2*y^2+1)*(x^2); 1274 def t1 = facFirstWeyl(h); 1275 //fist a correct list 1276 testNCfac(t1); 1277 //now a correct list with the factorized polynomial 1278 testNCfac(t1,h); 1279 //now we put in an incorrect list without a polynomial 1280 t1[3][3] = y; 1281 testNCfac(t1); 1282 //and last but not least we take h as additional input 1283 testNCfac(t1,h); 1228 1284 } 1229 1285 //================================================== … … 1237 1293 "USAGE: facSubWeyl(h,x,y); h, X, D polynomials 1238 1294 RETURN: list 1239 ASSUME: X and D are variables of a basering, which satisfy DX = XD +1. 1295 ASSUME: X and D are variables of a basering, which satisfy DX = XD +1. 1240 1296 @* That is, they generate the copy of the first Weyl algebra in a basering. 1241 1297 @* Moreover, h is a polynomial in X and D only. 1242 1298 PURPOSE: compute factorizations of the polynomial, which depends on X and D. 1243 1299 EXAMPLE: example facSubWeyl; shows examples 1244 SEE ALSO: facFirstWeyl, testNCfac 1300 SEE ALSO: facFirstWeyl, testNCfac, facFirstShift 1245 1301 "{ 1246 //proc facSubWeyl 1247 // basering can be anything having a Weyl algebra as subalgebra 1248 def r = basering; 1249 //We begin to check the input for assumptions 1250 // which are: X,D are vars of the basering, 1251 if ( (isVar(X)!=1) || (isVar(D)!=1) || (size(X)>1) || (size(D)>1) || 1252 (leadcoef(X) != number(1)) || (leadcoef(D) != number(1)) ) 1253 { 1254 ERROR("expected pure variables as generators of a subalgebra"); 1255 } 1256 // Weyl algebra: 1257 poly w = D*X-X*D-1; // [D,X]=1 1258 poly u = D*X-X*D+1; // [X,D]=1 1259 if (u*w!=0) 1260 { 1261 // that is no combination gives Weyl 1262 ERROR("2nd and 3rd argument do not generate a Weyl algebra"); 1263 } 1264 // one of two is correct 1265 if (u==0) 1266 { 1267 // hence w != 0, swap X<->D 1268 w = D; D=X; X=w; 1269 } 1270 // DONE with assumptions 1271 //Input successfully checked 1272 intvec lexpofX = leadexp(X); 1273 intvec lexpofD = leadexp(D); 1274 int varnumX=1; 1275 int varnumD=1; 1276 while(lexpofX[varnumX] != 1) 1277 { 1278 varnumX++; 1279 } 1280 while(lexpofD[varnumD] != 1) 1281 { 1282 varnumD++; 1283 } 1284 /* VL: it's not clear what you do here!!! I comment out */ 1285 // lexpofX = lexpofX +1; 1286 // lexpofX[varnumX] = 0; 1287 // lexpofX[varnumD] = 0; 1288 // if (deg(h,lexpofX)!=0) 1289 // { 1290 // return(list()); 1291 // } 1292 /* VL : to add printlevel stuff */ 1293 1294 ring firstweyl = 0,(var(varnumX),var(varnumD)),dp; 1295 def Firstweyl = nc_algebra(1,1); 1296 setring Firstweyl; 1297 poly h = imap(r,h); 1298 list result = facFirstWeyl(h); 1299 setring r; 1300 list result = imap(Firstweyl,result); 1301 return(result); 1302 //proc facSubWeyl 1303 // basering can be anything having a Weyl algebra as subalgebra 1304 def @r = basering; 1305 //We begin to check the input for assumptions 1306 // which are: X,D are vars of the basering, 1307 if ( (isVar(X)!=1) || (isVar(D)!=1) || (size(X)>1) || (size(D)>1) || 1308 (leadcoef(X) != number(1)) || (leadcoef(D) != number(1)) ) 1309 { 1310 ERROR("expected pure variables as generators of a subalgebra"); 1311 } 1312 // Weyl algebra: 1313 poly w = D*X-X*D-1; // [D,X]=1 1314 poly u = D*X-X*D+1; // [X,D]=1 1315 if (u*w!=0) 1316 { 1317 // that is no combination gives Weyl 1318 ERROR("2nd and 3rd argument do not generate a Weyl algebra"); 1319 } 1320 // one of two is correct 1321 int isReverted = 0; // Reverted Weyl if dx=xd-1 holds 1322 if (u==0) 1323 { 1324 // hence w != 0, swap X<->D later 1325 isReverted = 1; 1326 // w = D; D=X; X=w; 1327 } 1328 // else: do nothing 1329 // DONE with assumptions 1330 //Input successfuy checked 1331 intvec lexpofX = leadexp(X); 1332 intvec lexpofD = leadexp(D); 1333 int varnumX=1; 1334 int varnumD=1; 1335 while(lexpofX[varnumX] != 1) 1336 { 1337 varnumX++; 1338 } 1339 while(lexpofD[varnumD] != 1) 1340 { 1341 varnumD++; 1342 } 1343 /* VL: it's not clear what you do here!!! I comment out */ 1344 // lexpofX = lexpofX +1; 1345 // lexpofX[varnumX] = 0; 1346 // lexpofX[varnumD] = 0; 1347 // if (deg(h,lexpofX)!=0) 1348 // { 1349 // return(list()); 1350 // } 1351 /* VL : to add printlevel stuff */ 1352 1353 if (isReverted) 1354 { 1355 ring firstweyl = 0,(var(varnumD),var(varnumX)),dp; 1356 def Firstweyl = nc_algebra(1,1); 1357 setring Firstweyl; 1358 ideal M = 0:nvars(@r); 1359 M[varnumX]=var(2); 1360 M[varnumD]=var(1); 1361 map Q = @r,M; 1362 poly h= Q(h); 1363 } 1364 else 1365 { // that is unReverted 1366 ring firstweyl = 0,(var(varnumX),var(varnumD)),dp; 1367 def Firstweyl = nc_algebra(1,1); 1368 setring Firstweyl; 1369 poly h= imap(@r,h); 1370 } 1371 list result = facFirstWeyl(h); 1372 setring @r; 1373 list result; 1374 if (isReverted) 1375 { 1376 // map swap back 1377 ideal M; M[1] = var(varnumD); M[2] = var(varnumX); 1378 map S = Firstweyl, M; 1379 result = S(result); 1380 } 1381 else 1382 { 1383 // that is unReverted 1384 result = imap(Firstweyl,result); 1385 } 1386 return(result); 1302 1387 }//proc facSubWeyl 1303 1388 example 1304 1389 { 1305 "EXAMPLE:";echo=2; 1306 ring r = 0,(x,y,z),dp; 1307 matrix D[3][3]; D[1,3]=-1; 1308 def R = nc_algebra(1,D); // x,z generate Weyl subalgebra 1309 setring R; 1310 poly h = (x^2*z^2+x)*x; 1311 list fact1 = facSubWeyl(h,x,z); fact1[2]; 1312 poly test1 = fact1[2][2]; // for testing 1313 // compare with facFirstWeyl: 1314 ring s = 0,(z,x),dp; 1315 def S = nc_algebra(1,1); setring S; 1316 poly h = (x^2*z^2+x)*x; 1317 list fact2 = facFirstWeyl(h); fact2[2]; 1318 poly test2 = fact2[2][2]; // for testing 1319 map F = R,x,z; 1320 F(test1) - test2; // they are identical 1390 "EXAMPLE:";echo=2; 1391 ring r = 0,(x,y,z),dp; 1392 matrix D[3][3]; D[1,3]=-1; 1393 def R = nc_algebra(1,D); // x,z generate Weyl subalgebra 1394 setring R; 1395 poly h = (x^2*z^2+x)*x; 1396 list fact1 = facSubWeyl(h,x,z); 1397 // compare with facFirstWeyl: 1398 ring s = 0,(z,x),dp; 1399 def S = nc_algebra(1,1); setring S; 1400 poly h = (x^2*z^2+x)*x; 1401 list fact2 = facFirstWeyl(h); 1402 map F = R,x,0,z; 1403 list fact1 = F(fact1); // it is identical to list fact2 1404 testNCfac(fact1); // check the correctness again 1321 1405 } 1322 1406 //================================================== 1407 1408 //================================================== 1409 //************From here: Shift-Algebra************** 1410 //================================================== 1411 1412 //==================================================* 1413 //one factorization of a homogeneous polynomial 1414 //in the first Shift Algebra 1415 static proc homogfacFirstShift(poly h) 1416 {//proc homogfacFirstShift 1417 def r = basering; 1418 poly hath; 1419 int i; int j; 1420 if (!homogwithorder(h,intvec(0,1))) 1421 {//The given polynomial is not homogeneous 1422 return(list()); 1423 }//The given polynomial is not homogeneous 1424 if (h==0) 1425 { 1426 return(list(0)); 1427 } 1428 list result; 1429 int m = deg(h,intvec(0,1)); 1430 if (m>0) 1431 {//The degree is not zero 1432 hath = lift(var(2)^m,h)[1,1]; 1433 for (i = 1; i<=m;i++) 1434 { 1435 result = result + list(var(2)); 1436 } 1437 }//The degree is not zero 1438 else 1439 {//The degree is zero 1440 hath = h; 1441 }//The degree is zero 1442 ring tempRing = 0,(x),dp; 1443 setring tempRing; 1444 map thetamap = r,x,1; 1445 poly hath = thetamap(hath); 1446 list azeroresult = factorize(hath); 1447 list azeroresult_return_form; 1448 for (i = 1; i<=size(azeroresult[1]);i++) 1449 {//rewrite the result of the commutative factorization 1450 for (j = 1; j <= azeroresult[2][i];j++) 1451 { 1452 azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]); 1453 } 1454 }//rewrite the result of the commutative factorization 1455 setring(r); 1456 map finalmap = tempRing,var(1); 1457 list tempresult = finalmap(azeroresult_return_form); 1458 result = tempresult+result; 1459 return(result); 1460 }//proc homogfacFirstShift 1461 1462 //================================================== 1463 //Computes all possible homogeneous factorizations 1464 static proc homogfacFirstShift_all(poly h) 1465 {//proc HomogfacFirstShiftAll 1466 if (deg(h,intvec(1,1)) <= 0 ) 1467 {//h is a constant 1468 return(list(list(h))); 1469 }//h is a constant 1470 def r = basering; 1471 list one_hom_fac; //stands for one homogeneous factorization 1472 int i; int j; int k; 1473 int shiftcounter; 1474 //Compute again a homogeneous factorization 1475 one_hom_fac = homogfacFirstShift(h); 1476 one_hom_fac = delete(one_hom_fac,1); 1477 if (size(one_hom_fac) == 0) 1478 {//there is no homogeneous factorization or the polynomial was not homogeneous 1479 return(list()); 1480 }//there is no homogeneous factorization or the polynomial was not homogeneous 1481 list result = permpp(one_hom_fac); 1482 for (i = 1; i<=size(result);i++) 1483 { 1484 shiftcounter = 0; 1485 for (j = 1; j<=size(result[i]); j++) 1486 { 1487 if (result[i][j]==var(2)) 1488 { 1489 shiftcounter++; 1490 } 1491 else 1492 { 1493 result[i][j] = subst(result[i][j], var(1), var(1)-shiftcounter); 1494 } 1495 } 1496 result[i] = insert(result[i],1); 1497 } 1498 result = delete_dublicates_noteval(result); 1499 return(result); 1500 }//proc HomogfacFirstShiftAll 1501 1502 //================================================== 1503 //factorization of the first Shift Algebra 1504 proc facFirstShift(poly h) 1505 "USAGE: facFirstShift(h); h a polynomial in the first shift algebra 1506 RETURN: list 1507 PURPOSE: compute all factorizations of a polynomial in the first shift algebra 1508 THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle 1509 ASSUME: basering in the first shift algebra 1510 NOTE: Every entry of the output list is a list with factors for one possible factorization. 1511 EXAMPLE: example facFirstShift; shows examples 1512 SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl 1513 "{//facFirstShift 1514 if(nvars(basering)!=2) 1515 {//Our basering is the Shift algebra, but not the first 1516 ERROR("Basering is not the first shift algebra"); 1517 return(list()); 1518 }//Our basering is the Shift algebra, but not the first 1519 def r = basering; 1520 setring r; 1521 list LR = ringlist(r); 1522 number @n = leadcoef(LR[5][1,2]); 1523 poly @p = LR[6][1,2]; 1524 if ( @n!=number(1) ) 1525 { 1526 ERROR("Basering is not the first shift algebra"); 1527 return(list()); 1528 } 1529 list result = list(); 1530 int i;int j; int k; int l; //counter 1531 // create a ring with the ordering which makes shift algebra 1532 // graded 1533 // def r = basering; // done before 1534 ring tempRing = LR[1][1],(x,s),(a(0,1),Dp); 1535 def tempRingnc = nc_algebra(1,s); 1536 setring r; 1537 // information on relations 1538 if (@p == -var(1)) // reverted shift algebra 1539 { 1540 setring(tempRingnc); 1541 map transf = r, var(2), var(1); 1542 setring(r); 1543 map transfback = tempRingnc, var(2),var(1); 1544 // result = transfback(resulttemp); 1545 } 1546 else 1547 { 1548 if ( @p == var(2)) // usual shift algebra 1549 { 1550 setring(tempRingnc); 1551 map transf = r, var(1), var(2); 1552 // result = facshift(h); 1553 setring(r); 1554 map transfback = tempRingnc, var(1),var(2); 1555 } 1556 else 1557 { 1558 ERROR("Basering is not the first shift algebra"); 1559 return(list()); 1560 } 1561 } 1562 // main calls 1563 setring(tempRingnc); 1564 list resulttemp = facshift(transf(h)); 1565 setring(r); 1566 result = transfback(resulttemp); 1567 1568 list recursivetemp; 1569 for(i = 1; i<=size(result);i++) 1570 {//recursively factorize factors 1571 if(size(result[i])>2) 1572 {//Nontrivial factorization 1573 for (j=2;j<=size(result[i]);j++) 1574 {//Factorize every factor 1575 recursivetemp = facFirstShift(result[i][j]); 1576 if(size(recursivetemp)>1) 1577 {//we have a nontrivial factorization 1578 for(k=1; k<=size(recursivetemp);k++) 1579 {//insert factorized factors 1580 if(size(recursivetemp[k])>2) 1581 {//nontrivial 1582 result = insert(result,result[i],i); 1583 for(l = size(recursivetemp[k]);l>=2;l--) 1584 { 1585 result[i+1] = insert(result[i+1],recursivetemp[k][l],j); 1586 } 1587 result[i+1] = delete(result[i+1],j); 1588 }//nontrivial 1589 }//insert factorized factors 1590 }//we have a nontrivial factorization 1591 }//Factorize every factor 1592 }//Nontrivial factorization 1593 }//recursively factorize factors 1594 return(result); 1595 }//facFirstShift 1596 example 1597 { 1598 "EXAMPLE:";echo=2; 1599 ring R = 0,(x,s),dp; 1600 def r = nc_algebra(1,s); 1601 setring(r); 1602 poly h = (s^2*x+x)*s; 1603 facFirstShift(h); 1604 } 1605 1606 static proc facshift(poly h) 1607 "USAGE: facshift(h); h is a polynomial in the first Shift algebra 1608 RETURN: list 1609 PURPOSE: Computes a factorization of a polynomial h in the first Shift algebra 1610 THEORY: @code{facshift} returns a list with some factorizations of the given 1611 polynomial. The possibilities of the factorization of the highest 1612 homogeneous part and those of the lowest will be merged. If during this 1613 procedure a factorization of the polynomial occurs, it will be added to 1614 the output list. For a more detailled description visit 1615 @url{http://math.rwth-aachen.de/\~Albert.Heinle} 1616 SEE ALSO: homogfacFirstShift_all, homogfacFirstShift 1617 "{//proc facshift 1618 if(homogwithorder(h,intvec(0,1))) 1619 { 1620 return(homogfacFirstShift_all(h)); 1621 } 1622 def r = basering; 1623 map invo = basering,-var(1),-var(2); 1624 int i; int j; int k; 1625 intvec limits = deg(h,intvec(1,1)) ,deg(h,intvec(1,0)),deg(h,intvec(0,1)); 1626 def prod; 1627 //end finding the limits 1628 poly maxh = jet(h,deg(h,intvec(0,1)),intvec(0,1))-jet(h,deg(h,intvec(0,1))-1,intvec(0,1)); 1629 poly minh = jet(h,deg(h,intvec(0,-1)),intvec(0,-1))-jet(h,deg(h,intvec(0,-1))-1,intvec(0,-1)); 1630 list result; 1631 list temp; 1632 list homogtemp; 1633 list M; list hatM; 1634 list f1 = homogfacFirstShift_all(maxh); 1635 list f2 = homogfacFirstShift_all(minh); 1636 int is_equal; 1637 poly hath; 1638 for (i = 1; i<=size(f1);i++) 1639 {//Merge all combinations 1640 for (j = 1; j<=size(f2); j++) 1641 { 1642 M = M + mergence(f1[i],f2[j],limits); 1643 } 1644 }//Merge all combinations 1645 for (i = 1 ; i<= size(M); i++) 1646 {//filter valid combinations 1647 if (product(M[i]) == h) 1648 {//We have one factorization 1649 result = result + list(M[i]); 1650 M = delete(M,i); 1651 continue; 1652 }//We have one factorization 1653 else 1654 { 1655 if (deg(h,intvec(0,1))<=deg(h-product(M[i]),intvec(0,1))) 1656 { 1657 M = delete(M,i); 1658 continue; 1659 } 1660 if (deg(h,intvec(0,-1))<=deg(h-product(M[i]),intvec(0,-1))) 1661 { 1662 M = delete(M,i); 1663 continue; 1664 } 1665 } 1666 }//filter valid combinations 1667 M = delete_dublicates_eval(M); 1668 while(size(M)>0) 1669 {//work on the elements of M 1670 hatM = list(); 1671 for(i = 1; i<=size(M); i++) 1672 {//iterate over all elements of M 1673 hath = h-product(M[i]); 1674 temp = list(); 1675 //First check for common inhomogeneous factors between hath and h 1676 if (involution(NF(involution(hath,invo), std(involution(ideal(M[i][1]),invo))),invo)==0) 1677 {//hath and h have a common factor on the left 1678 j = 1; 1679 f1 = M[i]; 1680 if (j+1<=size(f1)) 1681 {//Checking for more than one common factor 1682 while(involution(NF(involution(hath,invo),std(involution(ideal(product(f1[1..(j+1)])),invo))),invo)==0) 1683 { 1684 if (j+1<size(f1)) 1685 { 1686 j++; 1687 } 1688 else 1689 { 1690 break; 1691 } 1692 } 1693 }//Checking for more than one common factor 1694 if (deg(product(f1[1..j]),intvec(1,1))!=0) 1695 { 1696 f2 = list(f1[1..j])+list(involution(lift(involution(product(f1[1..j]),invo),involution(hath,invo))[1,1],invo)); 1697 } 1698 else 1699 { 1700 f2 = list(f1[1..j])+list(involution(lift(product(f1[1..j]),involution(hath,invo))[1,1],invo)); 1701 } 1702 temp = temp + merge_cf(f2,f1,limits); 1703 }//hath and h have a common factor on the left 1704 if (reduce(hath, std(ideal(M[i][size(M[i])])))==0) 1705 {//hath and h have a common factor on the right 1706 j = size(M[i]); 1707 f1 = M[i]; 1708 if (j-1>0) 1709 {//Checking for more than one factor 1710 while(reduce(hath,std(ideal(product(f1[(j-1)..size(f1)]))))==0) 1711 { 1712 if (j-1>1) 1713 { 1714 j--; 1715 } 1716 else 1717 { 1718 break; 1719 } 1720 } 1721 }//Checking for more than one factor 1722 f2 = list(lift(product(f1[j..size(f1)]),hath)[1,1])+list(f1[j..size(f1)]); 1723 temp = temp + merge_cf(f2,M[i],limits); 1724 }//hath and h have a common factor on the right 1725 //and now the homogeneous 1726 maxh = jet(hath,deg(hath,intvec(0,1)),intvec(0,1))-jet(hath,deg(hath,intvec(0,1))-1,intvec(0,1)); 1727 minh = jet(hath,deg(hath,intvec(0,-1)),intvec(0,-1))-jet(hath,deg(hath,intvec(0,-1))-1,intvec(0,-1)); 1728 f1 = homogfacFirstShift_all(maxh); 1729 f2 = homogfacFirstShift_all(minh); 1730 for (j = 1; j<=size(f1);j++) 1731 { 1732 for (k=1; k<=size(f2);k++) 1733 { 1734 homogtemp = mergence(f1[j],f2[k],limits); 1735 } 1736 } 1737 for (j = 1; j<= size(homogtemp); j++) 1738 { 1739 temp = temp + mergence(homogtemp[j],M[i],limits); 1740 } 1741 for (j = 1; j<=size(temp); j++) 1742 {//filtering invalid entries in temp 1743 if(product(temp[j])==h) 1744 {//This is already a result 1745 result = result + list(temp[j]); 1746 temp = delete(temp,j); 1747 continue; 1748 }//This is already a result 1749 if (deg(hath,intvec(0,1))<=deg(hath-product(temp[j]),intvec(0,1))) 1750 { 1751 temp = delete(temp,j); 1752 continue; 1753 } 1754 }//filtering invalid entries in temp 1755 hatM = hatM + temp; 1756 }//iterate over all elements of M 1757 M = hatM; 1758 for (i = 1; i<=size(M);i++) 1759 {//checking for complete factorizations 1760 if (h == product(M[i])) 1761 { 1762 result = result + list(M[i]); 1763 M = delete(M,i); 1764 continue; 1765 } 1766 }//checking for complete factorizations 1767 M = delete_dublicates_eval(M); 1768 }//work on the elements of M 1769 //In the case, that there is none, write a constant factor before the factor of interest. 1770 for (i = 1 ; i<=size(result);i++) 1771 {//add a constant factor 1772 if (deg(result[i][1],intvec(1,1))!=0) 1773 { 1774 result[i] = insert(result[i],1); 1775 } 1776 }//add a constant factor 1777 result = delete_dublicates_noteval(result); 1778 return(result); 1779 }//proc facshift 1780 1323 1781 /* 1324 Example polynomials where one can find factorizations 1325 (x^2+y)*(x^2+y); 1326 (x^2+x)*(x^2+y); 1327 (x^3+x+1)*(x^4+y*x+2); 1328 (x^2*y+y)*(y+x*y); 1329 y^3+x*y^3+2*y^2+2*(x+1)*y^2+y+(x+2)*y; //Example 5 Gregoriev Schwarz. 1330 (y+1)*(y+1)*(y+x*y); //Landau Example projected to the first dimension. 1331 */ 1782 Example polynomials where one can find factorizations 1783 (x^2+y)*(x^2+y); 1784 (x^2+x)*(x^2+y); 1785 (x^3+x+1)*(x^4+y*x+2); 1786 (x^2*y+y)*(y+x*y); 1787 y^3+x*y^3+2*y^2+2*(x+1)*y^2+y+(x+2)*y; //Example 5 Grigoriev-Schwarz. 1788 (y+1)*(y+1)*(y+x*y); //Landau Example projected to the first dimension. 1789 */ 1790 1791 static proc bugFromMLee() 1792 { 1793 LIB "ncfactor.lib"; 1794 ring s = 0,(z,x),dp; 1795 def S = nc_algebra(1,1); setring S; 1796 poly f= -60z4x2-54z4-56zx3-59z2x-64; 1797 option(prot); option(mem); 1798 def l= facFirstWeyl (f); 1799 l; // before: empty list; after fix: 1 entry, f is irreducible 1800 poly g = 75z3x2+92z3+24; 1801 def l= facFirstWeyl (g); 1802 l; //before: empty list 1803 } 1804 1805 /* very hard things from Martin Lee: 1806 // ex1, ex2 1807 ring s = 0,(z,x),Ws(-1,1); 1808 def S = nc_algebra(1,1); setring S; 1809 poly a = 10z5x4+26z4x5+47z5x2-97z4x3; //Abgebrochen nach einer Stunde; yes, it takes long 1810 def l= facFirstWeyl (a); l; 1811 kill l; 1812 poly b = -5328z8x5-5328z7x6+720z9x2+720z8x3-16976z7x4-38880z6x5-5184z7x3-5184z6x4-3774z5x5+2080z8x+5760z7x2-6144z6x3-59616z5x4+3108z3x6-4098z6x2-25704z5x3-21186z4x4+8640z6x-17916z4x3+22680z2x5+2040z5x-4848z4x2-9792z3x3+3024z2x4-10704z3x2-3519z2x3+34776zx4+12096zx3+2898x4-5040z2x+8064x3+6048x2; //Abgebrochen nach 1.5 Stunden; seems to be very complicated 1813 def l= facFirstWeyl (b); l; 1814 1815 // ex3: there was difference in answers => fixed 1816 LIB "ncfactor.lib"; 1817 ring r = 0,(x,y,z),dp; 1818 matrix D[3][3]; D[1,3]=-1; 1819 def R = nc_algebra(1,D); 1820 setring R; 1821 poly g= 7*z4*x+62*z3+26*z; 1822 def l1= facSubWeyl (g, x, z); 1823 l1; 1824 //---- other ring 1825 ring s = 0,(x,z),dp; 1826 def S = nc_algebra(1,-1); setring S; 1827 poly g= 7*z4*x+62*z3+26*z; 1828 def l2= facFirstWeyl (g); 1829 l2; 1830 map F = R,x,0,z; 1831 list l1 = F(l1); 1832 l1; 1833 //---- so the answers look different, check them! 1834 testNCfac(l2); // ok 1835 testNCfac(l1); // was not ok, but now it's been fixed!!! 1836 1837 // selbst D und X so vertauschen dass sie erfuellt ist : ist gemacht 1838 1839 */
Note: See TracChangeset
for help on using the changeset viewer.
