/////////////////////////////////////////////////////////// version="version ncfactor.lib 4.0.0.0 Jun_2014 "; // $Id$ category="Noncommutative"; info=" LIBRARY: ncfactor.lib Tools for factorization in some noncommutative algebras AUTHORS: Albert Heinle, aheinle@uwaterloo.ca @* Viktor Levandovskyy, levandov@math.rwth-aachen.de OVERVIEW: In this library, new methods for factorization on polynomials are implemented for three types of algebras, all generated by 2*n, n in NN, generators (n'th Weyl, n'th shift and graded polynomials in n'th q-Weyl algebras) over a field K. @* More detailled description of the algorithms and related publications can be found at @url{https://cs.uwaterloo.ca/\~aheinle/}. PROCEDURES: facWeyl(h); Factorization in the n'th Weyl algebra facFirstWeyl(h); Factorization in the first Weyl algebra testNCfac(l[,h[,1]]); Tests factorizations from a given list for correctness facSubWeyl(h,X,D); Factorization in the first Weyl algebra as a subalgebra facShift(h); Factorization in the n'th shift algebra facFirstShift(h); Factorization in the first shift algebra homogfacNthQWeyl(h); Homogeneous factorization in the n'th Q-Weyl algebra homogfacFirstQWeyl(h); Homogeneous factorization in the first Q-Weyl algebra homogfacNthQWeyl_all(h); Homogeneous factorization (complete) in the n'th Q-Weyl algebra homogfacFirstQWeyl_all(h); Homogeneous factorization (complete) in the first Q-Weyl algebra tst_ncfactor(); Runs the examples of all contained not static functions. Test thing. "; LIB "rootsur.lib"; LIB "general.lib"; LIB "nctools.lib"; LIB "involut.lib"; LIB "freegb.lib"; // for isVar LIB "crypto.lib"; //for introot LIB "matrix.lib"; //for submatrix LIB "solve.lib"; //for solve LIB "poly.lib"; //for content proc tst_ncfactor() " A little test if the library works correct. Runs simply all examples of non-static functions. " { example facWeyl; example facShift; example facSubWeyl; example testNCfac; example homogfacFirstQWeyl; example homogfacFirstQWeyl_all; } example { "EXAMPLE:";echo=2; tst_ncfactor(); } ///////////////////////////////////////////////////// //==================================================* //deletes double-entries in a list of factorization //without evaluating the product. static proc delete_dublicates_noteval(list l) " INPUT: A list of lists; Output same as e.g. FacFirstWeyl. Containing different factorizations of a polynomial OUTPUT: If there are dublicates in this list, this procedure deletes them and returns the list without double entries " {//proc delete_dublicates_noteval list result= l; int j; int k; int i; int deleted = 0; int is_equal; for (i = 1; i<= size(l); i++) {//Iterate over the different factorizations for (j = i+1; j<= size(l); j++) {//Compare the i'th factorization to the j'th if (size(l[i])!= size(l[j])) {//different sizes => not equal j++; continue; }//different sizes => not equal is_equal = 1; for (k = 1; k <= size(l[i]);k++) {//Compare every entry if (l[i][k]!=l[j][k]) { is_equal = 0; break; } }//Compare every entry if (is_equal == 1) {//Delete this entry, because there is another equal one int the list result = delete(result, i-deleted); deleted = deleted+1; break; }//Delete this entry, because there is another equal one int the list }//Compare the i'th factorization to the j'th }//Iterate over the different factorizations return(result); }//proc delete_dublicates_noteval //================================================== static proc ncfactor_isWeyl() " INPUT: None OUTPUT: Returns 1 if the basering is a Weyl algebra, 0 otherwise. " {//ncfactor_isWeyl if (nvars(basering) % 2 != 0) {//Ring cannot be a Weyl algebra, as we need an even number of vars return(0); }//Ring cannot be a Weyl algebra, as we need an even number of vars list tempRingList = ringlist(basering); if (size(tempRingList)<=4) {//Given ring is commutative return(0); }//Given ring is commutative matrix C = tempRingList[5]; matrix D = tempRingList[6]; int i; int j; //Checking if C is okay. for (i = 1; i 1) {return(0);} countOnes = 0; for (j = i-1; j>=1; j--) {//Checking whether there is at most one nontrivial entry in each col if (D[j,i] == 1 or D[j,i]==-1) {countOnes++;} }//Checking whether there is at most one nontrivial entry in each col if (countOnes > 1) {return(0);} }//Iterating the diagonal //Now check if all variables are paired up. list the_vars; for (i = 1; i<=nvars(basering); i++) { the_vars[i] = var(i); } int noPairForFirstOne; while(size(the_vars)>0) {//We did not pair up all variables yet noPairForFirstOne = 1; for (i = 2; i<=size(the_vars); i++) {//Compare the commutation relations of the jth variable to the first one if (the_vars[1]*the_vars[i] - the_vars[i]*the_vars[1] == 1 or the_vars[1]*the_vars[i] - the_vars[i]*the_vars[1] == -1) {//Found a counter matching one noPairForFirstOne = 0; the_vars = delete(the_vars,i); the_vars = delete(the_vars,1); break; }//Found a counter matching one }//Compare the commutation relations of the jth variable to the first one if(noPairForFirstOne) {return(0);} }//We did not pair up all variables yet return(1); }//ncfactor_isWeyl //================================================== static proc ncfactor_isQWeyl() " INPUT: None OUTPUT: Returns 1 if the basering is a Weyl algebra, 0 otherwise. " {//ncfactor_isqWeyl if (nvars(basering) % 2 != 0) {//Ring cannot be a Weyl algebra, as we need an even number of vars return(0); }//Ring cannot be a Weyl algebra, as we need an even number of vars /* if (npars(basering)!=nvars(basering) div 2) */ /* {//Ring must have enough parameters */ /* return(0); */ /* }//Ring must have enough parameters */ list tempRingList = ringlist(basering); if (size(tempRingList)<=4) {//Given ring is commutative return(0); }//Given ring is commutative matrix C = tempRingList[5]; matrix D = tempRingList[6]; int i; int j; int k; int validentry; //Checking if C is okay. for (i = 1; i 1) {return(0);} countOnes = 0; for (j = i-1; j>=1; j--) {//Checking whether there is at most one nontrivial entry in each col if (D[j,i] == 1 or D[j,i]==-1) {countOnes++;} }//Checking whether there is at most one nontrivial entry in each col if (countOnes > 1) {return(0);} }//Iterating the diagonal //Now check if all variables are paired up. list the_vars; for (i = 1; i<=nvars(basering); i++) { the_vars[i] = var(i); } int noPairForFirstOne; while(size(the_vars)>0) {//We did not pair up all variables yet noPairForFirstOne = 1; for (i = 2; i<=size(the_vars); i++) {//Compare the commutation relations of the jth variable to the first one for (j = 1; j<=npars(basering);j++) { for (k = 1; k<= npars(basering); k++) { if (the_vars[1]*the_vars[i] - par(k)*the_vars[i]*the_vars[1] == 1 or par(k)*the_vars[1]*the_vars[i] - the_vars[i]*the_vars[1] == -1) {//Found a counter matching one noPairForFirstOne = 0; break; }//Found a counter matching one } if (noPairForFirstOne == 0) { break; } } if(noPairForFirstOne==0) { the_vars = delete(the_vars,i); the_vars = delete(the_vars,1); break; } }//Compare the commutation relations of the jth variable to the first one if(noPairForFirstOne) {return(0);} }//We did not pair up all variables yet return(1); }//ncfactor_isqWeyl //////////////////////////////////////////////////// //================================================== //////////////////////////////////////////////////// static proc ncfactor_isShift() " INPUT: None OUTPUT: Returns 1 if the basering is a Shift algebra, 0 otherwise. " {//ncfactor_isShift if (nvars(basering) % 2 != 0) {//Ring cannot be a Shift algebra, as we need an even number of vars return(0); }//Ring cannot be a Shift algebra, as we need an even number of vars list tempRingList = ringlist(basering); if (size(tempRingList)<=4) {//Given ring is commutative return(0); }//Given ring is commutative matrix C = tempRingList[5]; matrix D = tempRingList[6]; int i; int j; int k; //Checking if C is okay. for (i = 1; i 1) {return(0);} countOperators = 0; for (j = i-1; j>=1; j--) {//Checking whether there is at most one nontrivial entry in each col if (D[j,i]!= 0) {//Nontrivial relation detected isValidNCRelation = 0; for (k =1; k<= nvars(basering); k++) {//checking validity for each ring element if (D[j,i] == var(k) or D[j,i]==-var(k)) {//valid relation isValidNCRelation = 1; break; }//valid relation }//checking validity for each ring element if (!isValidNCRelation) {//We have not a valid nc-relation at hand return(0); }//We have not a valid nc-relation at hand else {//Increase the numbers of operators occurring in that row countOperators++; }//Increase the numbers of operators occurring in that row }//Nontrivial relation detected }//Checking whether there is at most one nontrivial entry in each col if (countOperators > 1) {return(0);} }//Iterating the diagonal //Now check if all variables are paired up. list the_vars; for (i = 1; i<=nvars(basering); i++) { the_vars[i] = var(i); } int noPairForFirstOne; while(size(the_vars)>0) {//We did not pair up all variables yet noPairForFirstOne = 1; for (i = 2; i<=size(the_vars); i++) {//Compare the commutation relations of the jth variable to the first one if (the_vars[1]*the_vars[i] - the_vars[i]*the_vars[1] == the_vars[1] or the_vars[1]*the_vars[i] - the_vars[i]*the_vars[1] == -the_vars[1] or the_vars[1]*the_vars[i] - the_vars[i]*the_vars[1] == the_vars[i] or the_vars[1]*the_vars[i] - the_vars[i]*the_vars[1] == -the_vars[i]) {//Found a counter matching one noPairForFirstOne = 0; the_vars = delete(the_vars,i); the_vars = delete(the_vars,1); break; }//Found a counter matching one }//Compare the commutation relations of the jth variable to the first one if(noPairForFirstOne) {return(0);} }//We did not pair up all variables yet return(1); }//ncfactor_isShift //================================================== static proc isInCommutativeSubRing(poly h) " INPUT: A polynomial h in an arbitrary polynomial ring. OUTPUT: 1, if all variables that appear at least in 1 monomial in h are all commuting, 0 else. " {//isInCommutativeSubRing list tempRingList = ringlist(basering); if (size(tempRingList)<=4) {//In this case, the given ring was commutative return(1); }//In this case, the given ring was commutative list appearing_variables; int i; int j; intvec degreeIntVec; for (i = 1; i<=nvars(basering);i++) {//checking for variables that appear degreeIntVec = 0:nvars(basering); degreeIntVec[i] = 1; if (deg(h, degreeIntVec) >0) {//Variable does appear in h appearing_variables = appearing_variables + list(var(i)); }//Variable does appear in h }//checking for variables that appear for (i = 1; i=size(result[j])) {//result[i] has more entries result = delete(result,j); continue; }//result[i] has more entries else {//result[j] has more entries result = delete(result,i); i--; break; }//result[j] has more entries }//There are two equal results; throw away that one with the smaller size }//comparing with the other elements }//Iterating over all elements in result return(result); }//proc delete_dublicates_eval //==================================================* static proc combinekfinlf(list g, int nof) //nof stands for "number of factors" " given a list of factors g and a desired size nof, this procedure combines the factors, such that we recieve a list of the length nof. INPUT: A list of containing polynomials or any type where the *-operator is existent OUTPUT: All possibilities (without permutation of the given list) to combine the polynomials into nof polynomials given by the user. " {//Procedure combinekfinlf list result; int i; int j; int k; //iteration variables list fc; //fc stands for "factors combined" list temp; //a temporary store for factors def nofgl = size(g); //nofgl stands for "number of factors of the given list" if (nofgl == 0) {//g was the empty list return(result); }//g was the empty list if (nof <= 0) {//The user wants to recieve a negative number or no element as a result return(result); }//The user wants to recieve a negative number or no element as a result if (nofgl == nof) {//There are no factors to combine result = result + list(g); return(result); }//There are no factors to combine if (nof == 1) {//User wants to get just one factor result = result + list(list(product(g))); return(result); }//User wants to get just one factor for (i = nof; i > 1; i--) {//computing the possibilities that have at least one original factor from g for (j = i; j>=1; j--) {//shifting the window of combinable factors to the left //fc below stands for "factors combined" fc = combinekfinlf(list(g[(j)..(j+nofgl - i)]),nof - i + 1); for (k = 1; k<=size(fc); k++) {//iterating over the different solutions of the smaller problem if (j>1) {//There are g_i before the combination if (j+nofgl -i < nofgl) {//There are g_i after the combination temp = list(g[1..(j-1)]) + fc[k] + list(g[(j+nofgl-i+1)..nofgl]); }//There are g_i after the combination else {//There are no g_i after the combination temp = list(g[1..(j-1)]) + fc[k]; }//There are no g_i after the combination }//There are g_i before the combination if (j==1) {//There are no g_i before the combination if (j+ nofgl -i 0) {//set of equal pairs is not empty temp = M[1]; temppos = 1; for (i = 2; i<=size(M); i++) {//finding the minimal element of M if (M[i][1]<=temp[1]) {//a possible candidate that is smaller than temp could have been found if (M[i][1]==temp[1]) {//In this case we must look at the second number if (M[i][2]< temp[2]) {//the candidate is smaller temp = M[i]; temppos = i; }//the candidate is smaller }//In this case we must look at the second number else {//The candidate is definately smaller temp = M[i]; temppos = i; }//The candidate is definately smaller }//a possible candidate that is smaller than temp could have been found }//finding the minimal element of M M = delete(M, temppos); if(temp[1]>1) {//There are factors to combine before the equal factor if (temp[1]0) {//There are factors to combine for (i = 1; i <= size(pre); i++) {//all possible pre's... for (j = 1; j<= size(post); j++) {//...combined with all possible post's candidate= pre[i]+list(f[temp[1]])+post[j]; if (limitcheck(candidate,limits)) { result = result + list(candidate); } }//...combined with all possible post's }//all possible pre's... }//There are factors to combine }//The most common case else {//the last factor is the common one pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits); for (i = 1; i<= size(pre); i++) {//iterating over the possible pre-factors candidate = pre[i]+list(f[temp[1]]); if (limitcheck(candidate,limits)) { result = result + list(candidate); } }//iterating over the possible pre-factors }//the last factor is the common one }//There are factors to combine before the equal factor else {//There are no factors to combine before the equal factor if (temp[1]0) {//we could find other combinations for (i = 1; i<=size(post); i++) { candidate = list(f[temp[1]])+post[i]; if (limitcheck(candidate,limits)) { result = result + list(candidate); } } }//we could find other combinations }//Just a check for security }//There are no factors to combine before the equal factor }//set of equal pairs is not empty for (i = 1; i <= size(result); i++) {//delete those combinations, who have an entry with degree less or equal 0 for (j = 1; j<=size(result[i]);j++) {//Delete entry if there is a zero or an integer as a factor if (deg(result[i][j]) <= 0) {//found one result = delete(result,i); i--; break; }//found one }//Delete entry if there is a zero as factor }//delete those combinations, who have an entry with degree less or equal 0 return(result); }//proc merge_cf //==================================================* //merges two sets of factors static proc mergence(list l1, list l2, intvec limits) " DEPRECATED " {//Procedure mergence list g; list f; int k; int i; int j; list F = list(); list G = list(); list tempEntry; list comb; if (size(l2)<=size(l1)) {//l1 will be our g, l2 our f g = l1; f = l2; }//l1 will be our g, l2 our f else {//l1 will be our f, l2 our g g = l2; f = l1; }//l1 will be our f, l2 our g if (size(f)==1 or size(g)==1) {//One of them just has one entry if (size(f)== 1) {f = list(1) + f;} if (size(g) == 1) {g = list(1) + g;} }//One of them just has one entry //first, we need to add some latent -1's to the list f and to the list g in order //to get really all possibilities of combinations later for (i=1;i<=size(f)-1;i++) {//first iterator for (j=i+1;j<=size(f);j++) {//second iterator tempEntry = f; tempEntry[i] = (-1)*tempEntry[i]; tempEntry[j] = (-1)*tempEntry[j]; F = F + list(tempEntry); }//secont iterator }//first iterator F = F + list(f); //And now same game with g for (i=1;i<=size(g)-1;i++) {//first iterator for (j=i+1;j<=size(g);j++) {//second iterator tempEntry = g; tempEntry[i] = (-1)*tempEntry[i]; tempEntry[j] = (-1)*tempEntry[j]; G = G + list(tempEntry); }//secont iterator }//first iterator G = G + list(g); //Done with that list result; for (i = 1; i<=size(F); i++) {//Iterate over all entries in F for (j = 1;j<=size(G);j++) {//Same with G comb = combinekfinlf(F[i],2,limits); for (k = 1; k<= size(comb);k++) {//for all possibilities of combinations of the factors of f result = result + merge_cf(comb[k],G[j],limits); result = result + merge_icf(comb[k],G[j],limits); result = delete_dublicates_noteval(result); }//for all possibilities of combinations of the factors of f }//Same with G }//Iterate over all entries in F return(result); }//Procedure mergence //================================================== //Checks, whether a list of factors doesn't exceed the given limits static proc limitcheck(list g, intvec limits) " DEPRECATED " {//proc limitcheck int i; if (size(limits)!=3) {//check the input return(0); }//check the input if(size(g)==0) { return(0); } def prod = product(g); intvec iv11 = intvec(1,1); intvec iv10 = intvec(1,0); intvec iv01 = intvec(0,1); def limg = intvec(deg(prod,iv11) ,deg(prod,iv10),deg(prod,iv01)); for (i = 1; i<=size(limg);i++) {//the final check if(limg[i]>limits[i]) { return(0); } }//the final check return(1); }//proc limitcheck //==================================================* //one factorization of a homogeneous polynomial //in the first Weyl Algebra static proc homogfacFirstWeyl(poly h) "USAGE: homogfacFirstWeyl(h); h is a homogeneous polynomial in the first Weyl algebra with respect to the weight vector [-1,1] RETURN: list PURPOSE: Computes a factorization of a homogeneous polynomial h with respect to the weight vector [-1,1] in the first Weyl algebra THEORY: @code{homogfacFirstWeyl} returns a list with a factorization of the given, [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with k positive, the last k entries in the output list are the second variable. If k is positive, the last k entries will be x. The other entries will be irreducible polynomials of degree zero or 1 resp. -1. SEE ALSO: homogfacFirstWeyl_all "{//proc homogfacFirstWeyl int p = printlevel-voice+2;//for dbprint def r = basering; poly hath; int i; int j; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} intvec ivm11 = intvec(-1,1); if (!homogwithorder(h,ivm11)) {//The given polynomial is not homogeneous ERROR("Given polynomial was not [-1,1]-homogeneous"); return(list()); }//The given polynomial is not homogeneous if (h==0) { return(list(0)); } list result; int m = deg(h,ivm11); dbprint(p,dbprintWhitespace +" Splitting the polynomial in A_0 and A_k-Part"); if (m!=0) {//The degree is not zero if (m <0) {//There are more x than y hath = lift(var(1)^(-m),h)[1,1]; for (i = 1; i<=-m; i++) { result = result + list(var(1)); } }//There are more x than y else {//There are more y than x hath = lift(var(2)^m,h)[1,1]; for (i = 1; i<=m;i++) { result = result + list(var(2)); } }//There are more y than x }//The degree is not zero else {//The degree is zero hath = h; }//The degree is zero dbprint(p,dbprintWhitespace+" Done"); //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1) list mons; dbprint(p,dbprintWhitespace+" Putting the monomials in the A_0-part in a list."); for(i = 1; i<=size(hath);i++) {//Putting the monomials in a list mons = mons+list(hath[i]); }//Putting the monomials in a list dbprint(p,dbprintWhitespace+" Done"); dbprint(p,dbprintWhitespace+" Mapping this monomials to K[theta]"); ring tempRing = 0,(x,y,theta),dp; setring tempRing; map thetamap = r,x,y; list mons = thetamap(mons); poly entry; for (i = 1; i<=size(mons);i++) {//transforming the monomials as monomials in theta entry = leadcoef(mons[i]); for (j = 0; j1;j--) {//drip x resp. y if (leftpart[j-1]==shiftvar) {//commutative j--; continue; }//commutative if (deg(leftpart[j-1],intvec(-1,1,0))!=0) {//stop here break; }//stop here //Here, we can only have a a0- part leftpart[j] = subst(leftpart[j-1],theta, theta + shift_sign); leftpart[j-1] = shiftvar; lparts = lparts + list(leftpart); }//drip x resp. y //and now deal with the right part if (rightpart[1] == x) { shift_sign = 1; shiftvar = x; } else { shift_sign = -1; shiftvar = y; } for (j = 1 ; j < size(rightpart); j++) { if (rightpart[j+1] == shiftvar) { j++; continue; } if (deg(rightpart[j+1],intvec(-1,1,0))!=0) { break; } rightpart[j] = subst(rightpart[j+1], theta, theta - shift_sign); rightpart[j+1] = shiftvar; rparts = rparts + list(rightpart); } //And now, we put all possibilities together tempadd = list(); for (j = 1; j<=size(lparts); j++) { for (k = 1; k<=size(rparts);k++) { tempadd = tempadd + list(lparts[j]+rparts[k]); } } tempadd = delete(tempadd,1); // The first entry is already in the list result = result + tempadd; continue; //We can may be not be done already with the ith entry }//One entry was theta resp. theta +1 }//checking every entry of result for theta or theta +1 dbprint(p,dbprintWhitespace +" Done"); //map back to the basering dbprint(p,dbprintWhitespace +" Mapping back everything to the basering"); setring(r); map finalmap = tempRing, var(1), var(2),var(1)*var(2); list result = finalmap(result); for (i=1; i<=size(result);i++) {//adding the K factor result[i] = k_factor + result[i]; }//adding the k-factor dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace +" Delete double entries in the list."); result = delete_dublicates_noteval(result); dbprint(p,dbprintWhitespace +" Done"); return(result); }//proc HomogfacFirstWeylAll /* example */ /* { */ /* "EXAMPLE:";echo=2; */ /* ring R = 0,(x,y),Ws(-1,1); */ /* def r = nc_algebra(1,1); */ /* setring(r); */ /* poly h = (x^2*y^2+1)*(x^4); */ /* homogfacFirstWeyl_all(h); */ /* } */ static proc homogfacNthWeyl_all(poly h) "USAGE: homogfacNthWeyl_all(h); h is a homogeneous polynomial in the nth Weyl algebra with respect to the ZZ-grading on the nth Weyl algebra. RETURN: list PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect to the ZZ-grading on the nth Weyl algebra THEORY: @code{homogfacFirstWeyl} returns a list with all factorization of the given, homogeneous polynomial. It uses the output of homogfacNthWeyl and permutes its entries with respect to the commutation rule. Furthermore, if a factor of degree zero is irreducible in K[theta_1, ..., theta_n], but reducible in the nth Weyl algebra, the permutations of this element with the other entries will also be computed. SEE ALSO: homogfacFirstWeyl "{//proc HomogfacNthWeylAll int p=printlevel-voice+2;//for dbprint intvec iv11= 1:nvars(basering); if (deg(h,iv11) <= 0 ) {//h is a constant dbprint(p,"Given polynomial was not homogeneous"); return(list(list(h))); }//h is a constant def r = basering; list one_hom_fac; //stands for one homogeneous factorization int i; int j; int k; int l; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} dbprint(p,dbprintWhitespace +" Calculate one homogeneous factorization using homogfacNthWeyl"); //Compute again a homogeneous factorization one_hom_fac = homogfacNthWeyl(h); dbprint(p,dbprintWhitespace +"Successful"); if (size(one_hom_fac) == 0) {//there is no homogeneous factorization or the polynomial was not homogeneous return(list()); }//there is no homogeneous factorization or the polynomial was not homogeneous //divide list in A0-Part and a list of x_i's resp. y_i's list list_not_azero = list(); list list_azero; list k_factor; int is_list_not_azero_empty = 1; int is_list_azero_empty = 1; k_factor = list(number(one_hom_fac[1])); dbprint(p, dbprintWhitespace + "Determine whether there is an A0 part or not."); int absValueOfDegree = 0; intvec degVecH = degreeOfNthWeylPoly(h); intvec lExp; for (i = 1; i<=size(degVecH); i++) {//adding up the absolute values of the degrees of the respective variables absValueOfDegree = absValueOfDegree + absValue(degVecH[i]); }//adding up the absolute values of the degrees of the respective variables if (absValueOfDegree < size(one_hom_fac) - 1) {//There is a nontrivial A0 part list_azero = one_hom_fac[2..(size(one_hom_fac)-absValueOfDegree)]; is_list_azero_empty = 0; }//There is a nontrivial A0 part dbprint(p,dbprintWhitespace +" Combine x_i,d_i to x_id_i in the factorization again."); dbprint(p,dbprintWhitespace + " The corresponding list of A0 factors is: " + string(list_azero)); for (i = 1; i1;j--) {//drip x resp. y if (leftpart[j-1]==shiftvar) {//commutative j--; continue; }//commutative if (leadexp(leftpart[j-1])[thetaIndex + nvars(r) div 2] - leadexp(leftpart[j-1])[thetaIndex]!=0) {//stop here break; }//stop here //Here, we can only have a a0- part leftpart[j] = subst(leftpart[j-1],theta(thetaIndex), theta(thetaIndex) + shift_sign); leftpart[j-1] = shiftvar; lparts = lparts + list(leftpart); }//drip x resp. y //and now deal with the right part if (rightpart[1] == x(thetaIndex)) { shift_sign = 1; shiftvar = x(thetaIndex); } else { shift_sign = -1; shiftvar = y(thetaIndex); } for (j = 1 ; j < size(rightpart); j++) { if (rightpart[j+1] == shiftvar) { j++; continue; } if (leadexp(rightpart[j+1])[thetaIndex + nvars(r) div 2] - leadexp(rightpart[j+1])[thetaIndex]!=0) { break; } rightpart[j] = subst(rightpart[j+1], theta(thetaIndex), theta(thetaIndex) - shift_sign); rightpart[j+1] = shiftvar; rparts = rparts + list(rightpart); } //And now, we put all possibilities together tempadd = list(); for (j = 1; j<=size(lparts); j++) { for (k = 1; k<=size(rparts);k++) { tempadd = tempadd + list(lparts[j]+rparts[k]); } } tempadd = delete(tempadd,1); // The first entry is already in the list result = result + tempadd; continue; //We can may be not be done already with the ith entry }//One entry was theta resp. theta +1 }//checking every entry of result for theta or theta +1 dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace + "The new result list is:"); dbprint(result); setring(r); ideal finalMapList; for(i = 1; i<=nvars(r);i++) { finalMapList[i] = var(i); } for (i = 1; i<=nvars(r) div 2; i++) { finalMapList[i + nvars(r)] = var(i)*var(i + (nvars(r) div 2)); } map finalmap = tempRing,finalMapList; list result = finalmap(result); for (i=1; i<=size(result);i++) {//adding the K factor result[i] = k_factor + result[i]; }//adding the k-factor dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace +" Delete double entries in the list."); result = delete_dublicates_noteval(result); dbprint(p,dbprintWhitespace +" Done"); return(result); }//proc HomogfacNthWeylAll /* Interesting Test-Polys: ring R = 0,(x1,x2,x3,d1,d2,d3),dp; matrix C[6][6] = 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1; matrix D[6][6] = 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1, -1,0,0,0,0,0, 0,-1,0,0,0,0, 0,0,-1,0,0,0; def r = nc_algebra(C,D); setring(r); poly h =x1*x2^2*x3^3*d1*d2^2+x2*x3^3*d2; h = (x1*x2^2*x3 + d1^2*d2*d3^3*x1^3*x2^3*x3^4)*(x1*d1 + x2*d2 + x3*d3); h = x1^2*d1+x1*x2*d2; */ //==================================================* //Computes all permutations of a given list static proc perm(list l) " DEPRECATED " {//proc perm int i; int j; list tempresult; list result; if (size(l)==0) { return(list()); } if (size(l)==1) { return(list(l)); } for (i = 1; i<=size(l); i++ ) { tempresult = perm(delete(l,i)); for (j = 1; j<=size(tempresult);j++) { tempresult[j] = list(l[i])+tempresult[j]; } result = result+tempresult; } return(result); }//proc perm //================================================== //computes all permutations of a given list by //ignoring equal entries (faster than perm) static proc permpp(list l) " INPUT: A list with entries of a type, where the ==-operator is defined OUTPUT: A list with all permutations of this given list. " {//proc permpp int i; int j; list tempresult; list l_without_double; list l_without_double_pos; int double_entry; list result; if (size(l)==0) { return(list()); } if (size(l)==1) { return(list(l)); } for (i = 1; i<=size(l);i++) {//Filling the list with unique entries double_entry = 0; for (j = 1; j<=size(l_without_double);j++) { if (l_without_double[j] == l[i]) { double_entry = 1; break; } } if (!double_entry) { l_without_double = l_without_double + list(l[i]); l_without_double_pos = l_without_double_pos + list(i); } }//Filling the list with unique entries for (i = 1; i<=size(l_without_double); i++ ) { tempresult = permpp(delete(l,l_without_double_pos[i])); for (j = 1; j<=size(tempresult);j++) { tempresult[j] = list(l_without_double[i])+tempresult[j]; } result = result+tempresult; } return(result); }//proc permpp //================================================== static proc checkIfProperNthWeyl() " INPUT: None OUTPUT: Checks whether the given basering is a proper Weyl algebra. Proper means in the sense of our algorithms, i.e. fulfilling the assumption that o ur basering is the Nth Weyl algebra and that the xs are the first n variables, the differential operators are the last n. Returns 1 if proper, 0 otherwise. " {//checkIfProperNthWeyl if (!ncfactor_isWeyl()) {return(0);} int i; for (i = 1; i<=nvars(basering) div 2; i++) { if (var(i + nvars(basering) div 2)*var(i) - var(i)*var(i+nvars(basering) div 2)!=1) { return(0); } } return(1); }//checkIfProperNthWeyl //================================================== static proc checkIfProperNthQWeyl() " INPUT: None OUTPUT: Checks whether the given basering is a proper q-Weyl algebra. Proper means in the sense of our algorithms, i.e. fulfilling the assumption that o ur basering is the Nth Weyl algebra and that the xs are the first n variables, the differential operators are the last n. Returns 1 if proper, 0 otherwise. " {//checkIfProperNthQWeyl if (!ncfactor_isQWeyl()) {return(0);} int i; for (i = 1; i<=nvars(basering) div 2; i++) { if (var(i + nvars(basering) div 2)*var(i) - par(i)*var(i)*var(i+nvars(basering) div 2)!=1) { return(0); } } return(1); }//checkIfProperNthQWeyl //================================================== static proc checkIfProperNthShift() " INPUT: None OUTPUT: Checks whether the given basering is a proper shift algebra. Proper means in the sense of our algorithms, i.e. fulfilling the assumption that our basering is the Nth shift algebra and that the xs are the first n variables, the shift operators are the last n. Returns 1 if proper, 0 otherwise. " {//checkIfProperNthShift if (!ncfactor_isShift()) {return(0);} int i; for (i = 1; i<=nvars(basering) div 2; i++) { if (var(i + nvars(basering) div 2)*var(i) - var(i)*var(i+nvars(basering) div 2)!=var(i+nvars(basering) div 2)) { return(0); } } return(1); }//checkIfProperNthShift //================================================== proc facWeyl(poly h) "USAGE: facWeyl(h); h a polynomial in the nth Weyl algebra RETURN: list PURPOSE: compute all factorizations of a polynomial in the first Weyl algebra THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle ASSUME: basering is the nth Weyl algebra, where n in NN. NOTE: Every entry of the output list is a list with factors for one possible factorization. The first factor is always a constant (1, if no nontrivial constant could be excluded). EXAMPLE: example facFirstWeyl; shows examples SEE ALSO: facSubWeyl, testNCfac, facFirstShift, facFirstWeyl "{//proc facWeyl //Definition of printlevel variable int p = printlevel-voice+2; int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} dbprint(p,dbprintWhitespace +" Checking if the given algebra is a Weyl algebra"); //Redefine the ring in my standard form if (!ncfactor_isWeyl()) {//Our basering is not the Weyl algebra ERROR("Ring was not a Weyl algebra"); return(list()); }//Our basering is not the Weyl algebra dbprint(p,dbprintWhitespace +" Successful"); //A last check before we start the real business: Is h maybe just //dependable on commutative variables? if (isInCommutativeSubRing(h)) {//h is in a commutative subring list hdepvars; intvec tempIntVec; for (i = 1; i<=nvars(basering) ; i++) { tempIntVec = 0:nvars(basering); tempIntVec[i] = 1; if (deg(h,tempIntVec)>0) { hdepvars = hdepvars + list(var(i)); } } if (size(hdepvars) ==0) {//We just have a constant return(list(list(h))); }//We just have a constant dbprint(p,dbprintWhitespace+"Polynomial was given commutative subring. Performing commutative factorization."); def r = basering; def rList = ringlist(basering); rList = delete(rList,5); rList = delete(rList,5); def tempRing = ring(rList); setring(tempRing); poly h = imap(r,h); list tempResult = factorize(h); list result = list(list()); int j; for (i = 1; i<=size(tempResult[1]); i++) { for (j = 1; j<=tempResult[2][i]; j++) { result[1] = result[1] + list(tempResult[1][i]); } } //mapping back setring(r); def result = imap(tempRing,result); dbprint(p,dbprintWhitespace+"result:"); dbprint(p,result); dbprint(p,dbprintWhitespace+"Computing all permutations of this factorization"); poly constantFactor = result[1][1]; result[1] = delete(result[1],1);//Deleting the constant factor result=permpp(result[1]); for (i = 1; i<=size(result);i++) {//Insert constant factor result[i] = insert(result[i],constantFactor); }//Insert constant factor dbprint(p,dbprintWhitespace+"Done."); return(result); }//h is in a commutative subring dbprint(p,dbprintWhitespace +" Successful"); list result = list(); int j; int k; int l; //counter if (!checkIfProperNthWeyl()) {//The given ring was not a proper nth Weyl algebra dbprint(p,dbprintWhitespace +" positions of the variables have to be switched"); dbprint(p,dbprintWhitespace + "Constructing the the proper ring."); def r = basering; list tempRingList = ringlist(r); tempRingList = delete(tempRingList,6); list the_vars; for (i = 1; i<=nvars(r); i++) {the_vars[i] = var(i);} int maybeDInWrongPos; poly tempVariable; for (i = 1; i<=size(the_vars) div 2; i++) {//Swapping the variables as needed maybeDInWrongPos = 1; if (the_vars[i + size(the_vars) div 2]*the_vars[i] -the_vars[i]*the_vars[i + size(the_vars) div 2] == 1) { i++; continue; } //If we enter this line, there is a break with our property //condition for (j = i+1; j<=size(the_vars); j++) { if (the_vars[j]*the_vars[i]-the_vars[i]*the_vars[j]==1) {//In this case, we matched a var x to a repective d tempVariable = the_vars[i + size(the_vars) div 2]; the_vars[i + size(the_vars) div 2] = the_vars[j]; the_vars[j] = tempVariable; maybeDInWrongPos = 0; break; }//In this case, we matched a var x to a repective d } if (maybeDInWrongPos) {//var(i) is actually a d, not an x print("i has to be pushed to the end."); tempVariable = the_vars[i]; the_vars = delete(the_vars, i); the_vars = the_vars + list(tempVariable); continue; }//var(i) is actually a d, not an x }//Swapping the variables as needed for (i = 1; i<=size(the_vars); i++) {tempRingList[2][i] = string(the_vars[i]);} matrix DTemp[nvars(r)][nvars(r)]; for (i = 1; i<=ncols(DTemp) div 2; i++) { DTemp[i,i + nvars(r) div 2] = 1; } tempRingList = tempRingList + list(DTemp); def tempRing = ring(tempRingList); dbprint(p,dbprintWhitespace + "Done. The altered ring is the following:"); dbprint(p,tempRing); setring(tempRing); poly h = imap(r,h); dbprint(p,dbprintWhitespace +" Successful"); list resulttemp = facWeyl(h); setring(r); result = imap(tempRing,resulttemp); return (result); }//The given ring was not a proper nth Weyl algebra dbprint(p, dbprintWhitespace +" factorization of the polynomial with the routine sfacwaNthWeyl"); result = sfacwaNthWeyl(h); dbprint(p,dbprintWhitespace +" Done"); if (homogwithorderNthWeyl(h)) { dbprint(p, dbprintWhitespace + " Polynomial was homogeneous, therefore we have already a complete factorization and do not have to go through the factors recursively."); return(result); } result = normalizeFactors(result); result = delete_dublicates_noteval(result); dbprint(p,dbprintWhitespace + "We have the following intermediate list of inhomogeneous factorizations:"); dbprint(p,result); dbprint(p,dbprintWhitespace +" recursively check factors for irreducibility"); list recursivetemp; int changedSomething; for(i = 1; i<=size(result);i++) {//recursively factorize factors if(size(result[i])>2) {//Nontrivial factorization for (j=2;j<=size(result[i]);j++) {//Factorize every factor recursivetemp = facWeyl(result[i][j]); //if(size(recursivetemp)>1) //{//we have a nontrivial factorization changedSomething = 0; for(k=1; k<=size(recursivetemp);k++) {//insert factorized factors if(size(recursivetemp[k])>2) {//nontrivial changedSomething = 1; result = insert(result,result[i],i); for(l = size(recursivetemp[k]);l>=2;l--) { result[i+1] = insert(result[i+1],recursivetemp[k][l],j); } result[i+1] = delete(result[i+1],j); }//nontrivial }//insert factorized factors if (changedSomething) { result = delete(result,i); } //}//we have a nontrivial factorization }//Factorize every factor }//Nontrivial factorization }//recursively factorize factors dbprint(p,dbprintWhitespace +" Done"); if (size(result)==0) {//only the trivial factorization could be found result = list(list(1,h)); }//only the trivial factorization could be found list resultWithInterchanges; dbprint(p,dbprintWhitespace+ "And the result without interchanges with homogeneous factors is:"); dbprint(p,result); for (i = 1; i <= size(result) ; i++) {//applying the interchanges to result resultWithInterchanges = resultWithInterchanges + checkForHomogInhomogInterchangabilityNthWeyl(result[i], 2, size(result[i])); }//applying the interchanges to result dbprint(p,dbprintWhitespace + "With interchanges, the result is:"); dbprint(p,resultWithInterchanges); //now, refine the possible redundant list return( delete_dublicates_noteval(resultWithInterchanges) ); }//proc facWeyl example { "EXAMPLE:";echo=2; ring R = 0,(x1,x2,d1,d2),dp; matrix C[4][4] = 1,1,1,1, 1,1,1,1, 1,1,1,1, 1,1,1,1; matrix D[4][4] = 0,0,1,0, 0,0,0,1, -1,0,0,0, 0,-1,0,0; def r = nc_algebra(C,D); setring(r); poly h = (d1+1)^2*(d1 + x1*d2); facWeyl(h); } //================================================== static proc normalizeFactors(list factList) "INPUT: A list of factorizations, as outputted e.g. by facWeyl OUTPUT: If any entry in a factorization is not primitive, this function divides the common divisor out and multiplies the first entry with it. " {//normalizeFactors int i; int j; list result = factList; for (i = 1; i<=size(result); i++) {//iterating through every different factorization for (j=2; j<=size(result[i]); j++) {//Iterating through all respective factors if (content(result[i][j])!=number(1)) {//Got one where the content is not equal to 1 result[i][1] = result[i][1] * content(result[i][j]); result[i][j] = result[i][j] / content(result[i][j]); }//Got one where the content is not equal to 1 }//Iterating through all respective factors }//iterating through every different factorization return(result); }//normalizeFactors //================================================== //factorization of the first Weyl Algebra //The following procedure just serves the purpose to //transform the input into an appropriate input for //the procedure sfacwa, where the ring must contain the //variables in a certain order. static proc facFirstWeyl_old(poly h) "USAGE: facFirstWeyl(h); h a polynomial in the first Weyl algebra RETURN: list PURPOSE: compute all factorizations of a polynomial in the first Weyl algebra THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle ASSUME: basering is the first Weyl algebra NOTE: Every entry of the output list is a list with factors for one possible factorization. The first factor is always a constant (1, if no nontrivial constant could be excluded). EXAMPLE: example facFirstWeyl; shows examples SEE ALSO: facSubWeyl, testNCfac, facFirstShift "{//proc facFirstWeyl_old //Definition of printlevel variable int p = printlevel-voice+2; int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} dbprint(p,dbprintWhitespace +" Checking if the given algebra is a Weyl algebra"); //Redefine the ring in my standard form if (!isWeyl()) {//Our basering is not the Weyl algebra ERROR("Ring was not the first Weyl algebra"); return(list()); }//Our basering is not the Weyl algebra dbprint(p,dbprintWhitespace +" Successful"); dbprint(p,dbprintWhitespace +" Checking, if the given ring is the first Weyl algebra"); if(nvars(basering)!=2) {//Our basering is the Weyl algebra, but not the first ERROR("Ring is not the first Weyl algebra"); return(list()); }//Our basering is the Weyl algebra, but not the first //A last check before we start the real business: Is h already given as a polynomial just //in one variable? if (deg(h,intvec(1,0))== 0 or deg(h,intvec(0,1)) == 0) {//h is in K[x] or in K[d] if (deg(h,intvec(1,0))== 0 and deg(h,intvec(0,1)) == 0) {//We just have a constant return(list(list(h))); }//We just have a constant dbprint(p,dbprintWhitespace+"Polynomial was given in one variable. Performing commutative factorization."); int theCommVar; if (deg(h,intvec(1,0)) == 0) {//The second variable is the variable to factorize theCommVar = 2; }//The second variable is the variable to factorize else{theCommVar = 1;} def r = basering; ring tempRing = 0,(var(theCommVar)),dp; if (theCommVar == 1){map mapToCommutative = r,var(1),1;} else {map mapToCommutative = r,1,var(1);} poly h = mapToCommutative(h); list tempResult = factorize(h); list result = list(list()); int j; for (i = 1; i<=size(tempResult[1]); i++) { for (j = 1; j<=tempResult[2][i]; j++) { result[1] = result[1] + list(tempResult[1][i]); } } //mapping back setring(r); map mapBackFromCommutative = tempRing,var(theCommVar); def result = mapBackFromCommutative(result); dbprint(p,dbprintWhitespace+"result:"); dbprint(p,result); dbprint(p,dbprintWhitespace+"Computing all permutations of this factorization"); poly constantFactor = result[1][1]; result[1] = delete(result[1],1);//Deleting the constant factor result=permpp(result[1]); for (i = 1; i<=size(result);i++) {//Insert constant factor result[i] = insert(result[i],constantFactor); }//Insert constant factor dbprint(p,dbprintWhitespace+"Done."); return(result); }//h is in K[x] or in K[d] dbprint(p,dbprintWhitespace +" Successful"); list result = list(); int j; int k; int l; //counter if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why { dbprint(p,dbprintWhitespace +" positions of the variables have to be switched"); def r = basering; ring tempRing = ringlist(r)[1][1],(x,y),Ws(-1,1); // very strange: // setting Wp(-1,1) leads to SegFault; to clarify why!!! def NTR = nc_algebra(1,1); setring NTR ; map transf = r, var(2), var(1); dbprint(p,dbprintWhitespace +" Successful"); list resulttemp = sfacwa(h); setring(r); map transfback = NTR, var(2),var(1); result = transfback(resulttemp); } else { dbprint(p, dbprintWhitespace +" factorization of the polynomial with the routine sfacwa"); result = sfacwa(h); dbprint(p,dbprintWhitespace +" Done"); } if (homogwithorder(h,intvec(-1,1))) { dbprint(p, dbprintWhitespace + " Polynomial was homogeneous, therefore we have already a complete factorization and do not have to go through the factors recursively."); return(result); } result = normalizeFactors(result); result = delete_dublicates_noteval(result); dbprint(p,dbprintWhitespace + "We have the following intermediate list of inhomogeneous factorizations:"); dbprint(p,result); dbprint(p,dbprintWhitespace +" recursively check factors for irreducibility"); list recursivetemp; int changedSomething; for(i = 1; i<=size(result);i++) {//recursively factorize factors if(size(result[i])>2) {//Nontrivial factorization for (j=2;j<=size(result[i]);j++) {//Factorize every factor recursivetemp = facFirstWeyl(result[i][j]); //if(size(recursivetemp)>1) //{//we have a nontrivial factorization changedSomething = 0; for(k=1; k<=size(recursivetemp);k++) {//insert factorized factors if(size(recursivetemp[k])>2) {//nontrivial changedSomething = 1; result = insert(result,result[i],i); for(l = size(recursivetemp[k]);l>=2;l--) { result[i+1] = insert(result[i+1],recursivetemp[k][l],j); } result[i+1] = delete(result[i+1],j); }//nontrivial }//insert factorized factors if (changedSomething) { result = delete(result,i); } //}//we have a nontrivial factorization }//Factorize every factor }//Nontrivial factorization }//recursively factorize factors dbprint(p,dbprintWhitespace +" Done"); if (size(result)==0) {//only the trivial factorization could be found result = list(list(1,h)); }//only the trivial factorization could be found list resultWithInterchanges; dbprint(p,dbprintWhitespace+ "And the result without interchanges with homogeneous factors is:"); dbprint(p,result); for (i = 1; i <= size(result) ; i++) {//applying the interchanges to result resultWithInterchanges = resultWithInterchanges + checkForHomogInhomogInterchangability(result[i],2,size(result[i])); }//applying the interchanges to result dbprint(p,dbprintWhitespace + "With interchanges, the result is:"); dbprint(p,resultWithInterchanges); //now, refine the possible redundant list return( delete_dublicates_noteval(resultWithInterchanges) ); }//proc facFirstWeyl_old proc facFirstWeyl(poly h) "USAGE: facFirstWeyl(h); h a polynomial in the first Weyl algebra RETURN: list PURPOSE: compute all factorizations of a polynomial in the first Weyl algebra THEORY: This function is a wrapper for facWeyl. It exists to make this library downward-compatible with older versions. ASSUME: basering is the first Weyl algebra NOTE: Every entry of the output list is a list with factors for one possible factorization. The first factor is always a constant (1, if no nontrivial constant could be excluded). EXAMPLE: example facFirstWeyl; shows examples SEE ALSO: facSubWeyl, testNCfac, facShift" {//facFirstWeyl return(facWeyl(h)); }//facFirstWeyl example { "EXAMPLE:";echo=2; ring R = 0,(x,y),dp; def r = nc_algebra(1,1); setring(r); poly h = (x^2*y^2+x)*(x+1); facFirstWeyl(h); } static proc checkForHomogInhomogInterchangabilityNthWeyl(list factors, posLeft, posRight) " INPUT: A list consisting of factors of a certain polynomial in the nth Weyl algebra, factors, and a position from the left and the right, where the last swap was done. OUTPUT: A list containing lists consisting of factors of a certain polynomial in the nth Weyl algebra. The purpose of this function is to check whether we can interchange certain inhomogeneous factors with homogeneous ones. If it is possible, this function returns a list of lists of possible interchanges. The idea came because of an example, where we need an extra swap in the end, otherwise we would not capture all factorizations. The example was h = x4d7+11x3d6+x2d7+x2d6+x3d4+29x2d5+xd6+8xd5+d6+5x2d3+14xd4+13d4+5xd2+d3+d; ASSUMPTIONS: - All factors are irreducible - Our basering is the Nth Weyl algebra; the xs are the first n variables, the differential operators are the last n. - No entry in the list factors is 0. " {//checkForHomogInhomogInterchangabilityNthWeyl int p = printlevel-voice+2; string dbprintWhitespace = ""; int i; int j; int k; int l; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (size(factors) <= 2 || posLeft >= posRight - 1) {//easiest case: There is nothing to swap return (list(factors)); }//easiest case: There is nothing to swap list result = list(factors); list tempResultEntries; list tempSwaps; list tempSwapsTempEntry; list attemptToSwap; int posHomogBegin; int posHomogEnd; poly leftHomogFactorProduct; poly rightHomogFactorProduct; dbprint(p, dbprintWhitespace+"We try to swap elements in the following list:"); dbprint(p, factors); dbprint(p, "The left border is at position: "+string(posLeft)); dbprint(p, "The right border is at position: " + string(posRight)); for (i = posLeft; i < posRight; i++) {//checking within the window posLeft <--> posRight, if there are interchanges possible leftHomogFactorProduct = 1; while (homogwithorderNthWeyl(factors[i])) {//We have a homogeneous polynomial somewhere on the left if (leftHomogFactorProduct == 1) {posHomogBegin = i;} leftHomogFactorProduct = leftHomogFactorProduct * factors[i]; i = i+1; if (i>=posRight) {break;} }//We have a homogeneous polynomial somewhere on the left if ((leftHomogFactorProduct !=1) && (!homogwithorderNthWeyl(factors[i]))) {//We have a group of homogeneous polynomials and an inhomogeneous one that we can try to swap. attemptToSwap = extractHomogeneousDivisorsRightNthWeyl(leftHomogFactorProduct*factors[i]); for (l = 1; l<=size(attemptToSwap); l++) { if (size(attemptToSwap[l])>1) {//Bingo, we were able to swap this one element dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" with its predecessors"); dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:"); dbprint(p,attemptToSwap); tempSwapsTempEntry = list(); for (j = size(factors); j >=1; j--) {//creating a new entry for the resulting list, replacing the swap in factors if (j==i) { for (k = size(attemptToSwap[l]); k >=1 ; k--) { tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[l][k]); } j=posHomogBegin;//Because of the change } else { tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]); } }//creating a new entry for the resulting list, replacing the swap in factors if (posRight <= size(tempSwapsTempEntry)) { tempSwaps = insert(tempSwaps,list(list(posHomogBegin+1,posRight),tempSwapsTempEntry)); } else { tempSwaps = insert(tempSwaps,list(list(posHomogBegin+1,size(tempSwapsTempEntry)), tempSwapsTempEntry)); } }//Bingo, we were able to swap this one element } }//We have a group of homogeneous polynomials and an inhomogeneous one that we can try to swap. else { if (i trying to swap j = i+1; // print(j); // print(size(factors)); // print(posRight); // print("==="); while (homogwithorderNthWeyl(factors[j])) { rightHomogFactorProduct = rightHomogFactorProduct * factors[j]; posHomogEnd = j; j = j+1; if (j > posRight) {break;} } attemptToSwap = extractHomogeneousDivisorsLeftNthWeyl(factors[i]*rightHomogFactorProduct); for (l =1; l<=size(attemptToSwap);l++) { if (size(attemptToSwap[l])>1) {//Bingo, we were able to swap this one element dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and its successors"); dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:"); dbprint(p,attemptToSwap); tempSwapsTempEntry = list(); for (j = size(factors); j >=1; j--) {//creating a new entry for the resulting list, replacing the swap in factors if (j==posHomogEnd) { for (k = size(attemptToSwap[l]); k >=1 ; k--) { tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[l][k]); } j = i; //Because we changed entry i+1 and i } else { tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]); } }//creating a new entry for the resulting list, replacing the swap in factors tempSwaps=insert(tempSwaps, list(list(posLeft,i+size(attemptToSwap[l])-1),tempSwapsTempEntry)); }//Bingo, we were able to swap this one element } }//position i+1 is homogeneous, position i is not ==> trying to swap } } }//checking within the window posLeft <--> posRight, if there are - interchanges possible //Now we will recursively call the function for all swapped entries. dbprint(p,dbprintWhitespace+ "Our list of different factorizations is now:"); dbprint(p,tempSwaps); for (i = 1; i<=size(tempSwaps);i++) {//recursive call to all formerly attempted swaps. dbprint(p, "Calling checkForHomogInterchangabilityNthWeyl recursively with values:"); dbprint(p, tempSwaps); tempResultEntries=checkForHomogInhomogInterchangabilityNthWeyl(tempSwaps[i][2], tempSwaps[i][1][1],tempSwaps[i][1][2]); result = result + tempResultEntries; }//recursive call to all formerly attempted swaps. result = delete_dublicates_noteval(result); return(result); }//checkForHomogInhomogInterchangabilityNthWeyl static proc checkForHomogInhomogInterchangability(list factors, posLeft, posRight) " INPUT: A list consisting of factors of a certain polynomial in the first Weyl algebra, factors, and a position from the left and the right, where the last swap was done. OUTPUT: A list containing lists consisting of factors of a certain polynomial in the first Weyl algebra. The purpose of this function is to check whether we can interchange certain inhomogeneous factors with homogeneous ones. If it is possible, this function returns a list of lists of possible interchanges. The idea came because of an example, where we need an extra swap in the end, otherwise we would not capture all factorizations. The example was h = x4d7+11x3d6+x2d7+x2d6+x3d4+29x2d5+xd6+8xd5+d6+5x2d3+14xd4+13d4+5xd2+d3+d; ASSUMPTIONS: - All factors are irreducible " {//checkForHomogInhomogInterchangability int p = printlevel-voice+2; string dbprintWhitespace = ""; int i; int j; int k; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (size(factors) <= 2 || posLeft >= posRight - 1) {//easiest case: There is nothing to swap return (list(factors)); }//easiest case: There is nothing to swap list result = list(factors); list tempResultEntries; list tempSwaps; list tempSwapsTempEntry; list attemptToSwap; intvec ivm11 = intvec(-1,1); dbprint(p, dbprintWhitespace+"We try to swap elements in the following list:"); dbprint(p, factors); for (i = posLeft; i < posRight; i++) {//checking within the window posLeft <--> posRight, if there are interchanges possible if (homogwithorder(factors[i],ivm11) && !homogwithorder(factors[i+1],ivm11)) {//position i is homogeneous, position i+1 is not ==> trying to swap attemptToSwap = extractHomogeneousDivisorsRight(factors[i]*factors[i+1]); if (size(attemptToSwap[1])>1) {//Bingo, we were able to swap this one element dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1)); dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:"); dbprint(p,attemptToSwap); tempSwapsTempEntry = list(); for (j = size(factors); j >=1; j--) {//creating a new entry for the resulting list, replacing the swap in factors if (j==i+1) { for (k = size(attemptToSwap[1]); k >=1 ; k--) { tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][k]); } j--; //Because we changed entry i+1 and i } else { tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]); } }//creating a new entry for the resulting list, replacing the swap in factors tempSwaps = insert(tempSwaps,list(list(i+1,posRight),tempSwapsTempEntry)); }//Bingo, we were able to swap this one element }//position i is homogeneous, position i+1 is not ==> trying to swap else { if(!homogwithorder(factors[i],ivm11) && homogwithorder(factors[i+1],ivm11)) {//position i+1 is homogeneous, position i is not ==> trying to swap attemptToSwap = extractHomogeneousDivisorsLeft(factors[i]*factors[i+1]); if (size(attemptToSwap[1])>1) {//Bingo, we were able to swap this one element dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1)); dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:"); dbprint(p,attemptToSwap); tempSwapsTempEntry = list(); for (j = size(factors); j >=1; j--) {//creating a new entry for the resulting list, replacing the swap in factors if (j==i+1) { for (k = size(attemptToSwap[1]); k >=1 ; k--) { tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][k]); } j--; //Because we changed entry i+1 and i } else { tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]); } }//creating a new entry for the resulting list, replacing the swap in factors tempSwaps = insert(tempSwaps,list(list(posLeft,i),tempSwapsTempEntry)); }//Bingo, we were able to swap this one element }//position i+1 is homogeneous, position i is not ==> trying to swap } }//checking within the window posLeft <--> posRight, if there are interchanges possible //Now we will recursively call the function for all swapped entries. dbprint(p,dbprintWhitespace+ "Our list of different factorizations is now:"); dbprint(p,tempSwaps); for (i = 1; i<=size(tempSwaps);i++) {//recursive call to all formerly attempted swaps. tempResultEntries=checkForHomogInhomogInterchangability(tempSwaps[i][2], tempSwaps[i][1][1],tempSwaps[i][1][2]); result = result + tempResultEntries; }//recursive call to all formerly attempted swaps. result = delete_dublicates_noteval(result); return(result); }//checkForHomogInhomogInterchangability static proc sfacwa(poly h) "INPUT: A polynomial h in the first Weyl algebra OUTPUT: A list of factorizations, where the factors might still be reducible. ASSUMPTIONS: - Our basering is the first Weyl algebra; the x is the first variable, the differential operator the second. " {//proc sfacwa int i; int j; int k; int p = printlevel-voice+2; string dbprintWhitespace = ""; number commonCoefficient = content(h); for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} dbprint(p,dbprintWhitespace + " Extracting homogeneous left and right factors"); if(homogwithorder(h,intvec(-1,1))) {//we are already dealing with a -1,1 homogeneous poly dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous. Start homog. fac. and ret. its result"); return(homogfacFirstWeyl_all(h)); }//we are already dealing with a -1,1 homogeneous poly list resulttemp = extractHomogeneousDivisors(h/commonCoefficient); //resulttemp = resulttemp + list(list(h/commonCoefficient)); list inhomogeneousFactorsToFactorize; int isAlreadyInInhomogList; dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace +" Making Set of inhomogeneous polynomials we have to factorize."); for (i = 1; i<=size(resulttemp); i++) {//Going through all different kinds of factorizations where we extracted homogeneous factors for (j = 1;j<=size(resulttemp[i]);j++) {//searching for the inhomogeneous factor if (!homogwithorder(resulttemp[i][j],intvec(-1,1))) {//We have found our candidate isAlreadyInInhomogList = 0; for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++) {//Checking if our candidate is already in our tofactorize-list if (inhomogeneousFactorsToFactorize[k]==resulttemp[i][j]) {//The candidate was already in the list isAlreadyInInhomogList = 1; break; }//The candidate was already in the list }//Checking if our candidate is already in our tofactorize-list if (!isAlreadyInInhomogList) { inhomogeneousFactorsToFactorize=inhomogeneousFactorsToFactorize + list(resulttemp[i][j]); } }//We have found our candidate }//searching for the inhomogeneous factor }//Going through all different kinds of factorizations where we extracted homogeneous factors dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace + "The set is:"); dbprint(p,inhomogeneousFactorsToFactorize); dbprint(p,dbprintWhitespace+ "Factorizing the different occuring inhomogeneous factors"); for (i = 1; i<= size(inhomogeneousFactorsToFactorize); i++) {//Factorizing all kinds of inhomogeneous factors inhomogeneousFactorsToFactorize[i] = sfacwa2(inhomogeneousFactorsToFactorize[i]); for (j = 1; j<=size(inhomogeneousFactorsToFactorize[i]);j++) {//Deleting the leading coefficient since we don't need him if (deg(inhomogeneousFactorsToFactorize[i][j][1],intvec(1,1))==0) { inhomogeneousFactorsToFactorize[i][j] = delete(inhomogeneousFactorsToFactorize[i][j],1); } }//Deleting the leading coefficient since we don't need him }//Factorizing all kinds of inhomogeneous factors dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace +" Putting the factorizations in the lists"); list result; int posInhomogPoly; int posInhomogFac; for (i = 1; i<=size(resulttemp); i++) {//going through all by now calculated factorizations for (j = 1;j<=size(resulttemp[i]); j++) {//Finding the inhomogeneous factor if (!homogwithorder(resulttemp[i][j],intvec(-1,1))) {//Found it posInhomogPoly = j; break; }//Found it }//Finding the inhomogeneous factor for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++) {//Finding the matching inhomogeneous factorization we already determined if(product(inhomogeneousFactorsToFactorize[k][1]) == resulttemp[i][j]) {//found it posInhomogFac = k; break; }//Found it }//Finding the matching inhomogeneous factorization we already determined for (j = 1; j <= size(inhomogeneousFactorsToFactorize[posInhomogFac]); j++) { result = insert(result, resulttemp[i]); result[1] = delete(result[1],posInhomogPoly); for (k =size(inhomogeneousFactorsToFactorize[posInhomogFac][j]);k>=1; k--) {//Inserting factorizations result[1] = insert(result[1],inhomogeneousFactorsToFactorize[posInhomogFac][j][k], posInhomogPoly-1); }//Inserting factorizations dbprint(p,dbprintWhitespace + "Added a factorization to result, namely:"); dbprint(p, result[1]); } }//going through all by now calculated factorizations dbprint(p,dbprintWhitespace +" Done"); result = delete_dublicates_noteval(result); for (i = 1; i<=size(result);i++) {//Putting the content everywhere result[i] = insert(result[i],commonCoefficient); }//Putting the content everywhere return(result); }//proc sfacwa static proc sfacwaNthWeyl(poly h) "INPUT: A polynomial h in the Nth Weyl algebra OUTPUT: A list of factorizations, where the factors might still be reducible. ASSUMPTIONS: - Our basering is the Nth Weyl algebra; the xs are the first n variables, the differential operators are the last n. " {//proc sfacwaNthWeyl int i; int j; int k; int p = printlevel-voice+2; string dbprintWhitespace = ""; number commonCoefficient = content(h); for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} dbprint(p,dbprintWhitespace + " Extracting homogeneous left and right factors"); if(homogwithorderNthWeyl(h)) {//we are already dealing with a -1,1 homogeneous poly dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous. Start homog. fac. and ret. its result"); return(homogfacNthWeyl_all(h)); }//we are already dealing with a -1,1 homogeneous poly list resulttemp = extractHomogeneousDivisorsNthWeyl(h/commonCoefficient); //resulttemp = resulttemp + list(list(h/commonCoefficient)); list inhomogeneousFactorsToFactorize; int isAlreadyInInhomogList; dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace +" Making Set of inhomogeneous polynomials we have to factorize."); for (i = 1; i<=size(resulttemp); i++) {//Going through all different kinds of factorizations where we extracted homogeneous factors for (j = 1;j<=size(resulttemp[i]);j++) {//searching for the inhomogeneous factor if (!homogwithorderNthWeyl(resulttemp[i][j])) {//We have found our candidate isAlreadyInInhomogList = 0; for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++) {//Checking if our candidate is already in our tofactorize-list if (inhomogeneousFactorsToFactorize[k]==resulttemp[i][j]) {//The candidate was already in the list isAlreadyInInhomogList = 1; break; }//The candidate was already in the list }//Checking if our candidate is already in our tofactorize-list if (!isAlreadyInInhomogList) { inhomogeneousFactorsToFactorize=inhomogeneousFactorsToFactorize + list(resulttemp[i][j]); } }//We have found our candidate }//searching for the inhomogeneous factor }//Going through all different kinds of factorizations where we extracted homogeneous factors dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace + "The set is:"); dbprint(p,inhomogeneousFactorsToFactorize); dbprint(p,dbprintWhitespace+ "Factorizing the different occuring inhomogeneous factors"); for (i = 1; i<= size(inhomogeneousFactorsToFactorize); i++) {//Factorizing all kinds of inhomogeneous factors inhomogeneousFactorsToFactorize[i] = sfacwa2NthWeyl(inhomogeneousFactorsToFactorize[i]); for (j = 1; j<=size(inhomogeneousFactorsToFactorize[i]);j++) {//Deleting the leading coefficient since we don't need him if (deg(inhomogeneousFactorsToFactorize[i][j][1],intvec(1,1))==0) { inhomogeneousFactorsToFactorize[i][j] = delete(inhomogeneousFactorsToFactorize[i][j],1); } }//Deleting the leading coefficient since we don't need him }//Factorizing all kinds of inhomogeneous factors dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace +" Putting the factorizations in the lists"); list result; int posInhomogPoly; int posInhomogFac; for (i = 1; i<=size(resulttemp); i++) {//going through all by now calculated factorizations for (j = 1;j<=size(resulttemp[i]); j++) {//Finding the inhomogeneous factor if (!homogwithorderNthWeyl(resulttemp[i][j])) {//Found it posInhomogPoly = j; break; }//Found it }//Finding the inhomogeneous factor for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++) {//Finding the matching inhomogeneous factorization we already determined if(product(inhomogeneousFactorsToFactorize[k][1]) == resulttemp[i][j]) {//found it posInhomogFac = k; break; }//Found it }//Finding the matching inhomogeneous factorization we already determined for (j = 1; j <= size(inhomogeneousFactorsToFactorize[posInhomogFac]); j++) { result = insert(result, resulttemp[i]); result[1] = delete(result[1],posInhomogPoly); for (k =size(inhomogeneousFactorsToFactorize[posInhomogFac][j]);k>=1; k--) {//Inserting factorizations result[1] = insert(result[1],inhomogeneousFactorsToFactorize[posInhomogFac][j][k], posInhomogPoly-1); }//Inserting factorizations dbprint(p,dbprintWhitespace + "Added a factorization to result, namely:"); dbprint(p, result[1]); } }//going through all by now calculated factorizations dbprint(p,dbprintWhitespace +" Done"); result = delete_dublicates_noteval(result); for (i = 1; i<=size(result);i++) {//Putting the content everywhere result[i] = insert(result[i],commonCoefficient); }//Putting the content everywhere return(result); }//proc sfacwaNthWeyl static proc sfacwa2(poly h) " Subprocedure of sfacwa Assumptions: - h is not in K[x] or in K[d], or even in K. These cases are caught by the input - The coefficients are integer values and the gcd of the coefficients is 1 " {//proc sfacwa2 int p=printlevel-voice+2; // for dbprint int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} intvec ivm11 = intvec(-1,1); intvec iv11 = intvec(1,1); intvec iv10 = intvec(1,0); intvec iv01 = intvec(0,1); intvec iv1m1 = intvec(1,-1); poly p_max; poly p_min; poly q_max; poly q_min; map invo = basering,-var(1),var(2); list calculatedRightFactors; if(homogwithorder(h,ivm11)) {//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous. Start homog. fac. and ret. its result"); return(homogfacFirstWeyl_all(h)); }//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here list result = list(); int j; int k; int l; dbprint(p,dbprintWhitespace+" Computing the degree-limits of the factorization"); //end finding the limits dbprint(p,dbprintWhitespace+" Computing the maximal and the minimal homogeneous part of the given polynomial"); list M = computeCombinationsMinMaxHomog(h); dbprint(p,dbprintWhitespace+" Done."); dbprint(p,dbprintWhitespace+" Filtering invalid combinations in M."); for (i = 1 ; i<= size(M); i++) {//filter valid combinations if (product(M[i]) == h) {//We have one factorization result = result + divides(M[i][1],h,invo,1); dbprint(p,dbprintWhitespace+"Result list updated:"); dbprint(p,dbprintWhitespace+string(result)); M = delete(M,i); continue; }//We have one factorization }//filter valid combinations dbprint(p,dbprintWhitespace+"Done."); dbprint(p,dbprintWhitespace+"The size of M is "+string(size(M))); for (i = 1; i<=size(M); i++) {//Iterate over all first combinations (p_max + p_min)(q_max + q_min) dbprint(p,dbprintWhitespace+" Combination No. "+string(i)+" in M:" ); p_max = jet(M[i][1],deg(M[i][1],ivm11),ivm11)-jet(M[i][1],deg(M[i][1],ivm11)-1,ivm11); p_min = jet(M[i][1],deg(M[i][1],iv1m1),iv1m1)-jet(M[i][1],deg(M[i][1],iv1m1)-1,iv1m1); q_max = jet(M[i][2],deg(M[i][2],ivm11),ivm11)-jet(M[i][2],deg(M[i][2],ivm11)-1,ivm11); q_min = jet(M[i][2],deg(M[i][2],iv1m1),iv1m1)-jet(M[i][2],deg(M[i][2],iv1m1)-1,iv1m1); dbprint(p,dbprintWhitespace+" pmax = "+string(p_max)); dbprint(p,dbprintWhitespace+" pmin = "+string(p_min)); dbprint(p,dbprintWhitespace+" qmax = "+string(q_max)); dbprint(p,dbprintWhitespace+" qmin = "+string(q_min)); //Check, whether p_max + p_min or q_max and q_min are already left or right divisors. if (divides(p_min + p_max,h,invo)) { dbprint(p,dbprintWhitespace+" Got one result."); result = result + divides(p_min + p_max,h,invo,1); } else { if (divides(q_min + q_max,h,invo)) { dbprint(p,dbprintWhitespace+" Got one result."); result = result + divides(q_min + q_max, h , invo, 1); } } //Now the check, if deg(p_max) = deg(p_min)+1 (and the same with q_max and q_min) if (deg(p_max, ivm11) == deg(p_min, ivm11) +1 or deg(q_max, ivm11) == deg(q_min, ivm11) +1 ) {//Therefore, p_max + p_min must be a left factor or we can dismiss the combination dbprint(p,dbprintWhitespace+" There are no homogeneous parts we can put between pmax and pmin resp. qmax and qmin."); //TODO: Prove, that then also a valid right factor is not possible M = delete(M,i); continue; }//Therefore, p_max + p_min must be a left factor or we can dismiss the combination //Done with the Check //If we come here, there are still homogeneous parts to be added to p_max + p_min //AND to q_max and q_min in //order to obtain a real factor //We use the procedure determineRestOfHomogParts to find our q. dbprint(p,dbprintWhitespace+" Solving for the other homogeneous parts in q"); calculatedRightFactors = determineRestOfHomogParts(p_max,p_min,q_max,q_min,h); dbprint(p,dbprintWhitespace+" Done with it. Found "+string(size(calculatedRightFactors)) +" solutions."); for (j = 1; j<=size(calculatedRightFactors);j++) {//Check out whether we really have right factors of h in calculatedRightFactors if (divides(calculatedRightFactors[j],h,invo)) { result = result + divides(calculatedRightFactors[j],h,invo,1); } else { dbprint(p,"Solution for max and min homog found, but not a divisor of h"); //TODO: Proof, why this can happen. } }//Check out whether we really have right factors of h in calculatedRightFactors }//Iterate over all first combinations (p_max + p_min)(q_max + q_min) result = delete_dublicates_noteval(result); //print(M); if (size(result) == 0) {//no factorization found result = list(list(h)); }//no factorization found return(result); }//proc sfacwa2 static proc sfacwa2NthWeyl(poly h) " Subprocedure of sfacwa Assumptions: - h is not part of a commutative subalgebra of the nth Weyl algebra - The coefficients are integer values and the gcd of the coefficients is 1 " {//proc sfacwa2NthWeyl int p=printlevel-voice+2; // for dbprint int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} poly p_max; poly p_min; poly q_max; poly q_min; ideal invoIdeal; intvec maxDegrees; intvec tempIntVec1; for (i = 1; i <= nvars(basering); i++) {//filling maxDegrees tempIntVec1 = 0:nvars(basering); tempIntVec1[i] = 1; maxDegrees[i] = deg(h,tempIntVec1); }//filling maxDegrees for (i = 1; i <= nvars(basering); i++) {//Filling the mapping rules for the involution map if (i <= nvars(basering) div 2) {invoIdeal[i] = -var(i);} else {invoIdeal[i] = var(i);} }//Filling the mapping rules for the involution map map invo = basering,invoIdeal; list calculatedRightFactors; if(homogwithorderNthWeyl(h)) {//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous. Start homog. fac. and ret. its result"); return(homogfacNthWeyl_all(h)); }//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here list result = list(); int j; int k; int l; list tempBetweenDegreesP; list tempBetweenDegreesQ; dbprint(p,dbprintWhitespace+" Computing the degree-limits of the factorization"); //end finding the limits dbprint(p,dbprintWhitespace+" Computing the maximal and the minimal homogeneous part of the given polynomial"); list M = computeCombinationsMinMaxHomogNthWeyl(h); dbprint(p,dbprintWhitespace+" Done."); dbprint(p,dbprintWhitespace+" Filtering invalid combinations in M."); for (i = 1 ; i<= size(M); i++) {//filter valid combinations if (product(M[i]) == h) {//We have one factorization result = result + divides(M[i][1],h,invo,1); dbprint(p,dbprintWhitespace+"Result list updated:"); dbprint(p,dbprintWhitespace+string(result)); M = delete(M,i); continue; }//We have one factorization }//filter valid combinations dbprint(p,dbprintWhitespace+"Done."); dbprint(p,dbprintWhitespace+"The size of M is "+string(size(M))); for (i = 1; i<=size(M); i++) {//Iterate over all first combinations (p_max + p_min)(q_max + q_min) dbprint(p,dbprintWhitespace+" Combination No. "+string(i)+" in M:" ); p_max = homogDistributionNthWeyl(M[i][1])[2][2]; p_min = homogDistributionNthWeyl(M[i][1])[1][2]; q_max = homogDistributionNthWeyl(M[i][2])[2][2]; q_min = homogDistributionNthWeyl(M[i][2])[1][2]; dbprint(p,dbprintWhitespace+" pmax = "+string(p_max)); dbprint(p,dbprintWhitespace+" pmin = "+string(p_min)); dbprint(p,dbprintWhitespace+" qmax = "+string(q_max)); dbprint(p,dbprintWhitespace+" qmin = "+string(q_min)); //Check, whether p_max + p_min or q_max and q_min are already left or right divisors. if (divides(p_min + p_max,h,invo)) { dbprint(p,dbprintWhitespace+" Got one result."); result = result + divides(p_min + p_max,h,invo,1); } else { if (divides(q_min + q_max,h,invo)) { dbprint(p,dbprintWhitespace+" Got one result."); result = result + divides(q_min + q_max, h , invo, 1); } } //Now the check, if deg(p_max) = deg(p_min)+1 (and the same with q_max and q_min) tempBetweenDegreesP = possibleHomogPartsInBetween(degreeOfNthWeylPoly(p_max), degreeOfNthWeylPoly(p_min), maxDegrees); tempBetweenDegreesQ = possibleHomogPartsInBetween(degreeOfNthWeylPoly(q_max), degreeOfNthWeylPoly(q_min), maxDegrees); if (size(tempBetweenDegreesQ)==2 or size(tempBetweenDegreesP)==2 ) {//Therefore, p_max + p_min must be a left factor or we can dismiss the combination dbprint(p,dbprintWhitespace+" There are no homogeneous parts we can put between pmax and pmin resp. qmax and qmin."); //TODO: Prove, that then also a valid right factor is not possible M = delete(M,i); continue; }//Therefore, p_max + p_min must be a left factor or we can dismiss the combination //Done with the Check //If we come here, there are still homogeneous parts to be added to p_max + p_min //AND to q_max and q_min in //order to obtain a real factor //We use the procedure determineRestOfHomogParts to find our q. dbprint(p,dbprintWhitespace+" Solving for the other homogeneous parts in q"); calculatedRightFactors = determineRestOfHomogPartsNthWeyl(p_max,p_min,q_max,q_min,h); dbprint(p,dbprintWhitespace+" Done with it. Found "+string(size(calculatedRightFactors)) +" solutions."); for (j = 1; j<=size(calculatedRightFactors);j++) {//Check out whether we really have right factors of h in calculatedRightFactors if (divides(calculatedRightFactors[j],h,invo)) { result = result + divides(calculatedRightFactors[j],h,invo,1); } else { dbprint(p,"Solution for max and min homog found, but not a divisor of h"); //TODO: Proof, why this can happen. } }//Check out whether we really have right factors of h in calculatedRightFactors }//Iterate over all first combinations (p_max + p_min)(q_max + q_min) result = delete_dublicates_noteval(result); //print(M); if (size(result) == 0) {//no factorization found result = list(list(h)); }//no factorization found return(result); }//proc sfacwa2NthWeyl static proc determineRestOfHomogParts(poly pmax, poly pmin, poly qmax, poly qmin, poly h) "INPUT: Polynomials p_max, p_min, q_max, q_min and h. The maximum homogeneous part h_max of h is given by p_max*pmin, the minimum homogeneous part h_min of h is given by p_min*q_min. OUTPUT: A list of right factors q of h that have q_max and q_min as their maximum respectively minimum homogeneous part. Empty list, if those elements are not existent ASSUMPTIONS: - deg(p_max,intvec(-1,1))>deg(p_min,intvec(-1,1)) +1 - deg(q_max,intvec(-1,1))>deg(q_min,intvec(-1,1)) +1 - p_max*q_max = h_max - p_min*q_min = h_min - The basering is the first Weyl algebra " {//proc determineRestOfHomogParts int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} int kappa = Min(intvec(deg(h,intvec(1,0)), deg(h,intvec(0,1)))); def R = basering; int n1 = deg(pmax,intvec(-1,1)); int nk = -deg(pmin,intvec(1,-1)); int m1 = deg(qmax,intvec(-1,1)); int ml = -deg(qmin,intvec(1,-1)); int j; int k; ideal mons; dbprint(p,dbprintWhitespace+" Extracting zero homog. parts of pmax, qmax, pmin, qmin and h."); //Extracting the zero homogeneous part of the given polynomials ideal pandqZero = pmax,pmin,qmax,qmin; if (n1 > 0){pandqZero[1] = lift(var(2)^n1,pmax)[1,1];} else{if (n1 < 0){pandqZero[1] = lift(var(1)^(-n1),pmax)[1,1];} else{pandqZero[1] = pmax;}} if (nk > 0){pandqZero[2] = lift(var(2)^nk,pmin)[1,1];} else{if (nk < 0){pandqZero[2] = lift(var(1)^(-nk),pmin)[1,1];} else{pandqZero[2] = pmin;}} if (m1 > 0){pandqZero[3] = lift(var(2)^m1,qmax)[1,1];} else{if (m1 < 0){pandqZero[3] = lift(var(1)^(-m1),qmax)[1,1];} else{pandqZero[3] = qmax;}} if (ml > 0){pandqZero[4] = lift(var(2)^ml,qmin)[1,1];} else{if (ml < 0){pandqZero[4] = lift(var(1)^(-ml),qmin)[1,1];} else{pandqZero[4] = qmin;}} list hZeroinR = homogDistribution(h); for (i = 1; i<=size(hZeroinR);i++) {//Extracting the zero homogeneous parts of the homogeneous summands of h if (hZeroinR[i][1] > 0){hZeroinR[i][2] = lift(var(2)^hZeroinR[i][1],hZeroinR[i][2])[1,1];} if (hZeroinR[i][1] < 0){hZeroinR[i][2] = lift(var(1)^(-hZeroinR[i][1]),hZeroinR[i][2])[1,1];} }//Extracting the zero homogeneous parts of the homogeneous summands of h dbprint(p,dbprintWhitespace+" Done!"); //Moving everything into the ring K[theta] dbprint(p,dbprintWhitespace+" Moving everything into the ring K[theta]"); ring KTheta = 0,(x,d,theta),dp; map thetamap = R, x, d; poly entry; ideal mons; ideal pandqZero; list hZeroinKTheta; setring(R); //Starting with p and q for (k=1; k<=4; k++) {//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys mons = ideal(); for(i = 1; i<=size(pandqZero[k]);i++) {//Putting the monomials in a list mons[size(mons)+1] = pandqZero[k][i]; }//Putting the monomials in a list setring(KTheta); mons = thetamap(mons); for (i = 1; i<=size(mons);i++) {//transforming the monomials as monomials in theta entry = leadcoef(mons[i]); for (j = 0; j= 1;k--) {//Transforming the different homogeneous parts of h into polys in K[theta] mons = ideal(); for(i = 1; i<=size(hZeroinR[k][2]);i++) {//Putting the monomials in a list mons[size(mons)+1] = hZeroinR[k][2][i]; }//Putting the monomials in a list setring(KTheta); mons = thetamap(mons); for (i = 1; i<=size(mons);i++) {//transforming the monomials as monomials in theta entry = leadcoef(mons[i]); for (j = 0; j."); setring(R); return(list()); }//No solution in this case. Return the empty list if(vdim(slimgb(solutionSystemforqs+theta))==-1) {//My conjecture is that this would never happen //ERROR("This is an counterexample to your conjecture. We have infinitely many solutions"); //TODO: See, what we would do here dbprint(p,dbprintWhitespace+"There are infinitely many solution to this system. We will return the empty list."); setring(R); return(list()); }//My conjecture is that this would never happen else {//We have finitely many solutions if(vdim(slimgb(solutionSystemforqs+theta))==1) {//exactly one solution for (i = 2; i<= size(qs)-1;i++) { qs[i][2] = NF(qs[i][2],solutionSystemforqs); } setring(R); map backFromSolutionRing = solutionRing,var(1)*var(2); list qs = backFromSolutionRing(qs); list result = list(0); for (i = 1; i<=size(qs); i++) { if (qs[i][1]>0){qs[i][2] = qs[i][2]*var(2)^qs[i][1];} if (qs[i][1]<0){qs[i][2] = qs[i][2]*var(1)^(-qs[i][1]);} result[1] = result[1] + qs[i][2]; } dbprint(p,dbprintWhitespace+"Found one unique solution. Returning the result."); return(result); }//exactly one solution else {//We have more than one solution, but finitely many def ringForSolveLib = solve(solutionSystemforqs+theta, "nodisplay"); setring ringForSolveLib; list valuesForQs = list(); list tempValues; int validSol; for (i = 1; i<=size(SOL); i++) {//filtering integer solutions validSol = 1; for (j =1; j<=size(SOL[i]); j++) {//Checking every entry if it is integer or not if (SOL[i][j] - int(SOL[i][j])!=0) {//No integer solution validSol = 0; break; }//No integer solution }//Checking every entry if it is integer or not if (validSol) { tempValues = list(); for (j = 1; j<=size(SOL[i]); j++) {//filling the valuesforQs tempValues[j]= int(SOL[i][j]); }//filling the valuesforQs valuesForQs = valuesForQs + list(tempValues); } }//filtering integer solutions print(valuesForQs); setring solutionRing; list differentQs = list(); ideal tempSolutionForQ; for (i = 1; i<=size(valuesForQs); i++) { differentQs[i] = qs; tempSolutionForQ = ideal(); for (j=2; j<=nvars(solutionRing);j++) {//filling solution ideal tempSolutionForQ[j-1] = var(j) - valuesForQs[i][j]; }//filling solution ideal for (j=2; j0){differentQs[k][i][2] = differentQs[k][i][2]*var(2)^differentQs[k][i][1];} if (differentQs[k][i][1]<0){differentQs[k][i][2] = differentQs[k][i][2]*var(1)^(-differentQs[k][i][1]);} result[k] = result[k] + differentQs[k][i][2]; } } dbprint(p,dbprintWhitespace+"Found multiple solutions. Returning the result."); return(result); }//We have more than one solution, but finitely many }//We have finitely many solutions }//proc determineRestOfHomogParts static proc determineRestOfHomogPartsNthWeyl(poly pmax, poly pmin, poly qmax, poly qmin, poly h) "INPUT: Polynomials p_max, p_min, q_max, q_min and h. The maximum homogeneous part h_max of h is given by p_max*pmin, the minimum homogeneous part h_min of h is given by p_min*q_min. OUTPUT: A list of right factors q of h that have q_max and q_min as their maximum respectively minimum homogeneous part. Empty list, if those elements are not existent ASSUMPTIONS: - deg(p_max) >_lex deg(p_min) +1 - deg(q_max) >_lex deg(q_min) +1 - p_max*q_max = h_max - p_min*q_min = h_min - The basering is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//proc determineRestOfHomogPartsNthWeyl int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i; int j; int k; int l; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} intvec kappa = 0:(nvars(basering) div 2); intvec maxDegrees; intvec tempIntVec1; intvec tempIntVec2; for (i = 1; i <= nvars(basering) div 2; i++) {//filling kappa tempIntVec1 = 0:nvars(basering); tempIntVec2 = 0:nvars(basering); tempIntVec1[i] = 1; tempIntVec2[i + nvars(basering) div 2] = 1; kappa[i] = Min(intvec(deg(h,tempIntVec1), deg(h,tempIntVec2))); }//filling kappa dbprint(p, dbprintWhitespace + "The Kappas for the respective variables are:"); dbprint(p,kappa); for (i = 1; i <= nvars(basering); i++) {//filling maxDegrees tempIntVec1 = 0:nvars(basering); tempIntVec1[i] = 1; maxDegrees[i] = deg(h,tempIntVec1); }//filling maxDegrees def R = basering; intvec n1 = degreeOfNthWeylPoly(pmax); intvec nk = degreeOfNthWeylPoly(pmin); list pBetweenDegrees = possibleHomogPartsInBetween(degreeOfNthWeylPoly(pmax), degreeOfNthWeylPoly(pmin), maxDegrees); int numberHomogPartsP = size(pBetweenDegrees); dbprint(p,dbprintWhitespace + "Possible degrees between pmax and pmin:"); dbprint(p,pBetweenDegrees); intvec m1 = degreeOfNthWeylPoly(qmax); intvec ml = degreeOfNthWeylPoly(qmin); list qBetweenDegrees = possibleHomogPartsInBetween(degreeOfNthWeylPoly(qmax), degreeOfNthWeylPoly(qmin), maxDegrees); int numberHomogPartsQ = size(qBetweenDegrees); dbprint(p,dbprintWhitespace + "Possible degrees between qmax and qmin:"); dbprint(p,qBetweenDegrees); ideal mons; list hBetweenDegrees = produceHomogListForProduct(pBetweenDegrees,qBetweenDegrees); int numberOfHomogPartsH = size(hBetweenDegrees); dbprint(p,dbprintWhitespace + "Possible degrees between hmax and hmin:"); dbprint(p,hBetweenDegrees); dbprint(p,dbprintWhitespace+" Extracting zero homog. parts of pmax, qmax, pmin, qmin and h."); //Extracting the zero homogeneous part of the given polynomials ideal pandqZero = pmax,pmin,qmax,qmin; for (i = 1; i<=nvars(basering) div 2; i++) {//Extracting zero homogeneous part variable for variable if (n1[i] > 0){pandqZero[1] = lift(var(i + nvars(basering) div 2)^n1[i],pandqZero[1])[1,1];} else{if (n1[i] < 0){pandqZero[1] = lift(var(i)^(-n1[i]),pandqZero[1])[1,1];} /* else{pandqZero[1] = pmax;} */} if (nk[i] > 0){pandqZero[2] = lift(var(i + nvars(basering) div 2)^nk[i],pandqZero[2])[1,1];} else{if (nk[i] < 0){pandqZero[2] = lift(var(i)^(-nk[i]),pandqZero[2])[1,1];} /* else{pandqZero[2] = pmin;} */} if (m1[i] > 0){pandqZero[3] = lift(var(i + nvars(basering) div 2)^m1[i],pandqZero[3])[1,1];} else{if (m1[i] < 0){pandqZero[3] = lift(var(i)^(-m1[i]),pandqZero[3])[1,1];} /* else{pandqZero[3] = qmax;} */} if (ml[i] > 0){pandqZero[4] = lift(var(i + nvars(basering) div 2)^ml[i],pandqZero[4])[1,1];} else{if (ml < 0){pandqZero[4] = lift(var(i)^(-ml[i]),pandqZero[4])[1,1];} /* else{pandqZero[4] = qmin;} */} }//Extracting zero homogeneous part variable for variable list hZeroinR = homogDistributionNthWeyl(h); for (i = 1; i<=size(hZeroinR);i++) {//Extracting the zero homogeneous parts of the homogeneous summands of h for (j = 1; j<=size(hZeroinR[i][1]); j++) {//Iterating through the different variables if (hZeroinR[i][1][j] > 0) { hZeroinR[i][2] = lift(var(j + nvars(basering) div 2)^hZeroinR[i][1][j],hZeroinR[i][2])[1,1]; } if (hZeroinR[i][1][j] < 0) { hZeroinR[i][2] = lift(var(j)^(-hZeroinR[i][1][j]),hZeroinR[i][2])[1,1]; } }//Iterating through the different variables }//Extracting the zero homogeneous parts of the homogeneous summands of h //Now we need to fill up the space between the h-parts: j = size(hZeroinR); hBetweenDegrees = reverse(hBetweenDegrees); list tempHList; /* for (i = size(hBetweenDegrees); i>0 ;i--) */ /* {//Filling up hzero with 0s in between */ /* ~; */ /* while (hBetweenDegrees[i] < hZeroinR[j][1]) */ /* { */ /* tempHList[i] = list(hBetweenDegrees[i],0); */ /* i--; */ /* } */ /* tempHList[i] = hZeroinR[j]; */ /* j--; */ /* }//Filling up hzero with 0s in between */ for (i = size(hBetweenDegrees); i>0 ;i--) { tempHList[i] = list(hBetweenDegrees[i],0); } for (i = 1; i <= size(hZeroinR); i++) { for (j = 1; j<=size(tempHList);j++) { if (hZeroinR[i][1] == tempHList[j][1]) { tempHList[j][2] = hZeroinR[i][2]; } } } hZeroinR = tempHList; //hBetweenDegrees = reverse(hBetweenDegrees); dbprint(p,dbprintWhitespace+" Done!"); dbprint(p,dbprintWhitespace+" Moving everything into the ring K[theta]"); ring KTheta = 0,(x(1..(nvars(basering) div 2)), d(1..(nvars(basering) div 2)), theta(1..(nvars(basering) div 2))),dp; setring(KTheta); ideal mapList; for (i = 1; i<=nvars(R) ; i++) {//filling the list of elements we want to map mapList[i] = var(i); }//filling the list of elements we want to map map thetamap = R, mapList; poly entry; ideal mons; ideal pandqZero; list hZeroinKTheta; intvec lExp; setring(R); //Starting with p and q for (k=1; k<=4; k++) {//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys mons = ideal(); for(i = 1; i<=size(pandqZero[k]);i++) {//Putting the monomials in a list mons[size(mons)+1] = pandqZero[k][i]; }//Putting the monomials in a list setring(KTheta); mons = thetamap(mons); for (i = 1; i<=size(mons);i++) {//transforming the monomials as monomials in the theta_i entry = leadcoef(mons[i]); lExp = leadexp(mons[i]); for (l = 1; l <=nvars(R) div 2; l++) { for (j = 0; j= 1;k--) {//Transforming the different homogeneous parts of h into polys in K[theta] mons = ideal(); for(i = 1; i<=size(hZeroinR[k][2]);i++) {//Putting the monomials in a list mons[size(mons)+1] = hZeroinR[k][2][i]; }//Putting the monomials in a list setring(KTheta); mons = thetamap(mons); for (i = 1; i<=size(mons);i++) {//transforming the monomials as monomials in theta entry = leadcoef(mons[i]); lExp = leadexp(mons[i]); for (l = 1; l <=nvars(R) div 2; l++) { for (j = 0; j=1; i--) { ps[i] = list(pBetweenDegrees[i],0); } ps[1][2] = pandqZero[2]; ps[numberHomogPartsP][2] = pandqZero[1]; list qs; for (i = numberHomogPartsQ; i>=1 ; i--) { qs[i] = list(qBetweenDegrees[i],0); } qs[1][2] = pandqZero[4]; qs[numberHomogPartsQ][2] = pandqZero[3]; qs = reverse(qs); ps = reverse(ps); pBetweenDegrees = reverse(pBetweenDegrees); qBetweenDegrees = reverse(qBetweenDegrees); tempIntVec2 = 0:nvars(R); for (i = 1; i<=size(kappa); i++) {tempIntVec2[i + size(kappa)] = kappa[i];} list coefficientIndices = possibleHomogPartsInBetween(kappa, 0:size(kappa), tempIntVec2); poly tempPoly; for (i = 2; i."); setring(R); return(list()); }//No solution in this case. Return the empty list if(vdim(slimgb(solutionSystemforqs + theThetas))==-1) {//My conjecture is that this would never happen //ERROR("This is an counterexample to your conjecture. We have infinitely many solutions"); //TODO: See, what we would do here dbprint(p,dbprintWhitespace+"There are infinitely many solution to this system. We will return the empty list."); setring(R); return(list()); }//My conjecture is that this would never happen else {//We have finitely many solutions if(vdim(slimgb(solutionSystemforqs+theThetas))==1) {//exactly one solution for (i = 2; i<= size(qs)-1;i++) { qs[i][2] = NF(qs[i][2],slimgb(solutionSystemforqs)); } setring(R); ideal finalMapList; for (i = 1; i<=nvars(R) div 2; i++) { finalMapList[i] = var(i)*var(i + (nvars(R) div 2)); } map backFromSolutionRing = solutionRing,finalMapList; list qs = backFromSolutionRing(qs); list result = list(0); for (i = 1; i<=size(qs); i++) { for (j= 1; j<=size(qs[i][1]); j++) { if (qs[i][1][j]>0){qs[i][2] = qs[i][2]*var(j + nvars(R) div 2)^qs[i][1][j];} if (qs[i][1][j]<0){qs[i][2] = qs[i][2]*var(j)^(-qs[i][1][j]);} } result[1] = result[1] + qs[i][2]; } dbprint(p,dbprintWhitespace+"Found one unique solution. Returning the result."); return(result); }//exactly one solution else {//We have more than one solution, but finitely many def ringForSolveLib = solve(solutionSystemforqs+theThetas, "nodisplay"); setring ringForSolveLib; list valuesForQs = list(); list tempValues; int validSol; for (i = 1; i<=size(SOL); i++) {//filtering integer solutions validSol = 1; for (j =1; j<=size(SOL[i]); j++) {//Checking every entry if it is integer or not if (SOL[i][j] - int(SOL[i][j])!=0) {//No integer solution validSol = 0; break; }//No integer solution }//Checking every entry if it is integer or not if (validSol) { tempValues = list(); for (j = 1; j<=size(SOL[i]); j++) {//filling the valuesforQs tempValues[j]= int(SOL[i][j]); }//filling the valuesforQs valuesForQs = valuesForQs + list(tempValues); } }//filtering integer solutions setring solutionRing; list differentQs = list(); ideal tempSolutionForQ; for (i = 1; i<=size(valuesForQs); i++) { differentQs[i] = qs; tempSolutionForQ = ideal(); for (j=(nvars(R) div 2)+1; j<=nvars(solutionRing);j++) {//filling solution ideal tempSolutionForQ[j-(nvars(R) div 2)] = var(j) - valuesForQs[i][j]; }//filling solution ideal for (j=2; j0) {differentQs[k][i][2] = differentQs[k][i][2]*var(j + nvars(R) div 2)^differentQs[k][i][1][j];} if (differentQs[k][i][1][j]<0) {differentQs[k][i][2] = differentQs[k][i][2]*var(j)^(-differentQs[k][i][1][j]);} } result[k] = result[k] + differentQs[k][i][2]; } } dbprint(p,dbprintWhitespace+"Found multiple solutions. Returning the result."); return(result); }//We have more than one solution, but finitely many }//We have finitely many solutions }//proc determineRestOfHomogPartsNthWeyl static proc produceHomogListForProduct(list homogParts1, list homogParts2) " INPUT: Two lists of integer vectors, homogParts1 and homogParts2, where the contained integer vectors all have the same size n. OUTPUT: One list of integer vectors, which have size n respectively. The entries are the set of possible integer vectors which result of a sum of an entry y in homogParts1 and an entry x in homogParts2. " {//produceHomogListForProduct int i; int j; int k; list p = reverse(sort(homogParts1)[1]); list q = reverse(sort(homogParts2)[1]); list result; intvec tempIntVec; for (i = 1; i<= size(homogParts1); i++) { for (j = 1; j<= size(homogParts2); j++) { tempIntVec = homogParts1[i] + homogParts2[j]; if (binarySearch(result,tempIntVec)==0) {//Element was not yet in result, add it result = result + list(tempIntVec); result = sort(result)[1]; }//Element was not yet in result, add it } } result = reverse(result); //Till here, we have a complete list with entries that are between hmax and hmin //Now, we need to filter this list. int isIn; for(i = 1; i<=size(p); i++) {//Checking, if homogparts[i] has a counterpart in homogparts2 isIn = 0; for (j = 1; j<=size(q); j++) {//iterating through the counterparts if(p[i] + q[j] == result[i]) { isIn = 1; break; } }//iterating through the counterparts if(!isIn) { result = delete(result,i); continue; } }//Checking, if homogparts[i] has a counterpart in homogparts2 for (i = 0; i inputElem) { last = middle -1; } else { return(middle); } } } return(0); }//binarySearch static proc possibleHomogPartsInBetween(intvec pmax, intvec pmin, intvec maxDegrees) " INPUT: Integer vectors pmax, pmin and maxDegrees. All but maxDegrees are of the same size, except from maxDegrees, which is double of the size of the rest. OUTPUT: A list of possible points in Z^n that can lie between pmax and pmin, sorted with respect to the lexicographic order from biggest to smallest (leftmost entry in the intvec is the biggest), represented as list of intvecs. " {//possibleHomogPartsInBetween if (size(pmax) != size(pmin) || 2*size(pmax)!=size(maxDegrees) || pmax<=pmin) { ERROR("pmax and pmin shall have the same size, and maxDegrees shall be double the size of both. Furhtermore, it must hold that pmax>pmin. The Input was: " + string(pmax) + ", " + string(pmin) + ", " + string(maxDegrees)); } list result; intvec tempIntVec = pmax; while(tempIntVec > pmin) { result = result + list(tempIntVec); tempIntVec = nextSmallerEntry(tempIntVec,pmin,maxDegrees); } result = result + list(pmin); return (sort(result)[1]); }//possibleHomogPartsInBetween static proc nextSmallerEntry(intvec pmax, intvec pmin, intvec maxDegrees) "INPUT: Integer vectors pmax, pmin and maxDegrees. maxDegrees has to have size 2*size(pmax), where pmin has to have the same size as pmax. OUTPUT: Counting down by lexicographical ordering, this function returns the next smaller entry after pmax, which is still bigger than pmin. " {//nextSmallerEntry if (pmax == pmin) {return (pmin);} if (size(pmax) == 1) { if (pmax <= pmin) {return(pmin);} else {return(pmax -1);} } int i; intvec recPmax = pmax[1..(size(pmax)-1)]; intvec recPmin = pmin[1..(size(pmin)-1)]; intvec recMaxDegrees = intvec(maxDegrees[1]); for (i = 2; i<=size(maxDegrees); i++) { if (i!=size(pmax) && i!= size(2*size(pmax))) { recMaxDegrees = recMaxDegrees,maxDegrees[i]; } } if (recPmax == recPmin) {//In this case, we can only possibly count down at our current position if (pmax[size(pmax)] > pmin[size(pmin)]) {return (recPmax,(pmax[size(pmax)]-1));} else {return(pmin);} }//In this case, we can only possibly count down at our current position else {//In this case, we can go down to the bounds given by maxDegrees if (pmax[size(pmax)] > -maxDegrees[size(pmax)]) { return (recPmax,(pmax[size(pmax)]-1)); } else { return(nextSmallerEntry(recPmax,recPmin,recMaxDegrees),maxDegrees[2*size(pmax)]); } }//In this case, we can go down to the bounds given by maxDegrees }//nextSmallerEntry static proc possibleHomogPartsInBetweenNonRecursive(intvec pmax, intvec pmin, intvec maxDegrees) " INPUT: Integer vectors pmax and maxDegrees. All but maxDegrees are of the same size, except from maxDegrees, which is double of the size of the rest. OUTPUT: A list of possible points in Z^n that can lie between pmax and pmin, sorted with respect to the lexicographic order from biggest to smallest (leftmost entry in the intvec is the biggest), represented as list of intvecs. " {//possibleHomogPartsInBetween if (size(pmax) != size(pmin) || 2*size(pmax)!=size(maxDegrees) || pmax<=pmin) { ERROR("pmax and pmin shall have the same size, and maxDegrees shall be double the size of both. Furhtermore, it must hold that pmax>pmin. The Input was: " + string(pmax) + ", " + string(pmin) + ", " + string(maxDegrees)); } list result = list(pmax); int pos; int i; int j; int k; list leftPart; int leftPartIsMaxAndMin; list possibleMiddles; list possibleRightParts = list(list()); list tempRightParts; intvec tempentry; for (pos = size(pmax); pos >=1 ; pos--) { leftPart = list(); leftPartIsMaxAndMin = 1; possibleMiddles = list(); for (i = 1; i < pos; i++) {//filling the left part leftPart[i] = pmax[i]; if(pmax[i]!=pmin[i]) {leftPartIsMaxAndMin = 0;} }//filling the left part if (leftPartIsMaxAndMin) { for (i = pmax[pos]-1; i>pmin[pos]; i--) {//possible entries for position pos possibleMiddles = possibleMiddles + list(i); }//possible entries for position pos } else { for (i = pmax[pos]-1; i>=-maxDegrees[pos]; i--) {//possible entries for position pos possibleMiddles = possibleMiddles + list(i); }//possible entries for position pos } for (i = 1; i<=size(possibleMiddles); i++) {//Adding possibilities to result for (j = 1; j<=size(possibleRightParts); j++) {//going through the right parts tempentry = intvec(0); for (k = 1; k<=size(leftPart);k++) { tempentry[k] = leftPart[k]; } if (size(leftPart)==0) {tempentry = possibleMiddles[i];} else {tempentry = tempentry,possibleMiddles[i];} for (k = 1; k<=size(possibleRightParts[j]) ; k++) { tempentry = tempentry, possibleRightParts[j][k]; } result = result + list(tempentry); }//going through the right parts }//Adding possibilities to result tempRightParts = list(); for (i = 1; i<=size(possibleRightParts);i++) { for (j = -maxDegrees[pos]; j <= maxDegrees[pos + size(pmax)]; j++) { tempRightParts = tempRightParts + list(list(j) + possibleRightParts[i]); } } possibleRightParts = tempRightParts; } //Now the last possible guys result = result + list(pmin); leftPart = list(); int positionWhereDifference = 1; for (i = 1; i<= size(pmin); i++) { if(pmax[i] != pmin[i]) { positionWhereDifference = i; break; } } possibleRightParts = list(list()); for (pos = size(pmin); pos > positionWhereDifference; pos--) {//Dealing with the minimal parts that we left out in the loop above leftPart = list(); possibleMiddles = list(); for (i = 1; i < pos; i++) {//filling up the left part leftPart[i] = pmin[i]; }//filling up the left part for (i = pmin[pos] +1 ; i<=maxDegrees[pos + size(pmin)] ; i++) { possibleMiddles = possibleMiddles + list(i); } for (i = 1; i<=size(possibleMiddles); i++) {//Adding possibilities to result for (j = 1; j<=size(possibleRightParts); j++) {//going through the right parts tempentry = intvec(0); for (k = 1; k<=size(leftPart);k++) { tempentry[k] = leftPart[k]; } if (size(leftPart)==0) {tempentry = possibleMiddles[i];} else {tempentry = tempentry,possibleMiddles[i];} for (k = 1; k<=size(possibleRightParts[j]) ; k++) { tempentry = tempentry, possibleRightParts[j][k]; } result = result + list(tempentry); }//going through the right parts }//Adding possibilities to result tempRightParts = list(); for (i = 1; i<=size(possibleRightParts);i++) { for (j = -maxDegrees[pos]; j <= maxDegrees[pos + size(pmax)]; j++) { tempRightParts = tempRightParts + list(list(j) + possibleRightParts[i]); } } possibleRightParts = tempRightParts; }//Dealing with the minimal parts that we left out in the loop above result = sort(result)[1]; if (result[1] != pmin) {result = insert(result,pmin);} return(result); }//possibleHomogPartsInBetween static proc gammaForTheta(int j1,int j2) " INPUT: Two integers j1 and j2 OUTPUT: A polynomial in the first variable of the given ring. It calculates the following function: / 1, if j1,j2>0 or j1,j2 <= 0 | prod_{kappa = 0}^{|j1|-1}(var(1)-kappa), if j1<0, j2>0, |j1|<=|j2| gamma_{j1,j2}:= < prod_{kappa = 0}^{j2-1}(var(1)-kappa-|j1|+|j2|), if j1<0, j2>0, |j1|>|j2| | prod_{kappa = 1}^{j1}(var(1)+kappa), if j1>0, j2<0, |j1|<=|j2| \ prod_{kappa = 1}^{|j2|}(\var(1)+kappa+|j1|-|j2|), if j1>0, j2<0, |j1|>|j2| ASSUMPTION: - Ring has at least one variable " {//gammaForTheta if (j1<=0 && j2 <=0) {return(1);} if (j1>=0 && j2 >=0){return(1);} poly result; int i; if (j1<0 && j2>0) {//case 2 or 3 from description above if (absValue(j1)<=absValue(j2)) {//Case 2 holds here result = 1; for (i = 0;i nvars(basering)) {ERROR("The size of the input vectors is bigger than the number of variables.");} int i; int j; ideal separateGammas; poly polyForPositionJ; for (j = 1; j <= size(j1); j++) {//Iterate through the \theta_j if (j1[j]<=0 && j2[j] <=0) { separateGammas[j] = 1; j++; continue; } if (j1[j]>=0 && j2[j] >=0) { separateGammas[j] = 1; j++; continue; } if (j1[j]<0 && j2[j]>0) {//case 2 or 3 from description above if (absValue(j1[j])<=absValue(j2[j])) {//Case 2 holds here polyForPositionJ = 1; for (i = 0;i=1;j--) { hath = result[i][j]*hath; tempList = extractHomogeneousDivisorsRight(hath); if(size(tempList[1])==1) {//We could not swap this element to the right break; }//We could not swap this element to the right dbprint(p,dbprintWhitespace+"A swapping (left) of an element was possible"); for(k = 1; k<=size(tempList);k++) { tempResult = insert(tempResult,result[i]); for (l = j;l<=posInhomog;l++) { tempResult[1] = delete(tempResult[1],j); } for (l = size(tempList[k]);l>=1;l--) { tempResult[1] = insert(tempResult[1],tempList[k][l],j-1); } } } hath = result[i][posInhomog]; for(j=posInhomog+1;j<=size(result[i]);j++) { hath = hath*result[i][j]; tempList = extractHomogeneousDivisorsLeft(hath); if(size(tempList[1])==1) {//We could not swap this element to the right break; }//We could not swap this element to the right dbprint(p,dbprintWhitespace+"A swapping (right) of an element was possible"); for(k = 1; k<=size(tempList);k++) { tempResult = insert(tempResult,result[i]); for (l=posInhomog; l<=j;l++) { tempResult[1] = delete(tempResult[1],posInhomog); } for (l = size(tempList[k]);l>=1;l--) { tempResult[1] = insert(tempResult[1],tempList[k][l],posInhomog-1); } } } }//Checking if we can swap left resp. right divisors result = result + tempResult; result = delete_dublicates_noteval(result); return(result); }//extractHomogeneousDivisors static proc extractHomogeneousDivisorsNthWeyl(poly h) "INPUT: A polynomial h in the nth Weyl algebra OUTPUT: If h is homogeneous with respect to the ZZ-grading on the nth Weyl algebra, then all factorizations of h are returned. If h is inhomogeneous, then a list l is returned whose entries are again lists k = [k_1,...,k_n], where k_1*...*k_n = h and there exists an i in {1,...,n}, such that k_i is inhomogeneous and there is no homogeneous polynomial that divides this k_i neither from the left nor from the right. All the other entries in k are homogeneous polynomials. GENERAL ASSUMPTIONS: - The basering is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//extractHomogeneousDivisorsNthWeyl int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i; int j; int k; int l; list result; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (homogwithorderNthWeyl(h)) {//given polynomial was homogeneous already dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations."); result = homogfacNthWeyl_all(h); for (i = 1; i<=size(result);i++) {//removing the first entry (coefficient) from the list result result[i] = delete(result[i],1); }//removing the first entry (coefficient) from the list result return(result); }//given polynomial was homogeneous already dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous left divisors extracted"); list leftDivisionPossibilities = extractHomogeneousDivisorsLeftNthWeyl(h); dbprint(p,dbprintWhitespace+"Done. The result is:"); dbprint(p,leftDivisionPossibilities); dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous Right divisors extracted"); list rightDivisionPossibilities = extractHomogeneousDivisorsRightNthWeyl(h); dbprint(p,dbprintWhitespace+"Done. The result is:"); dbprint(p,rightDivisionPossibilities); list tempList; dbprint(p,dbprintWhitespace+"Calculating remaining right and left homogeneous divisors"); for (i = 1; i<=size(leftDivisionPossibilities); i++) {//iterating through the list with extracted left divisors tempList = extractHomogeneousDivisorsRightNthWeyl (leftDivisionPossibilities[i][size(leftDivisionPossibilities[i])]); leftDivisionPossibilities[i] = delete(leftDivisionPossibilities[i], size(leftDivisionPossibilities[i])); for (j=1;j<=size(tempList);j++) {//Updating the list for Result tempList[j] = leftDivisionPossibilities[i] + tempList[j]; }//Updating the list for Result result = result + tempList; }//iterating through the list with extracted left divisors for (i = 1; i<=size(rightDivisionPossibilities); i++) {//iterating through the list with extracted right divisors tempList = extractHomogeneousDivisorsLeftNthWeyl(rightDivisionPossibilities[i][1]); rightDivisionPossibilities[i] = delete(rightDivisionPossibilities[i],1); for (j=1;j<=size(tempList);j++) {//Updating the list for Result tempList[j] = tempList[j]+rightDivisionPossibilities[i]; }//Updating the list for Result result = result + tempList; }//iterating through the list with extracted right divisors dbprint(p,dbprintWhitespace+"Done"); int posInhomog; poly hath = 1; list tempResult; dbprint(p,dbprintWhitespace+"Checking if we can swap left resp. right divisors and updating result."); for (i = 1; i<= size(result); i++) {//Checking if we can swap left resp. right divisors for (j = 1; j<=size(result[i]);j++) {//finding the position of the inhomogeneous element in the list if(!homogwithorderNthWeyl(result[i][j])) { posInhomog = j; break; } }//finding the position of the inhomogeneous element in the list hath = result[i][posInhomog]; for(j=posInhomog-1;j>=1;j--) { hath = result[i][j]*hath; tempList = extractHomogeneousDivisorsRightNthWeyl(hath); if(size(tempList[1])==1) {//We could not swap this element to the right break; }//We could not swap this element to the right dbprint(p,dbprintWhitespace+"A swapping (left) of an element was possible"); for(k = 1; k<=size(tempList);k++) { tempResult = insert(tempResult,result[i]); for (l = j;l<=posInhomog;l++) { tempResult[1] = delete(tempResult[1],j); } for (l = size(tempList[k]);l>=1;l--) { tempResult[1] = insert(tempResult[1],tempList[k][l],j-1); } } } hath = result[i][posInhomog]; for(j=posInhomog+1;j<=size(result[i]);j++) { hath = hath*result[i][j]; tempList = extractHomogeneousDivisorsLeftNthWeyl(hath); if(size(tempList[1])==1) {//We could not swap this element to the right break; }//We could not swap this element to the right dbprint(p,dbprintWhitespace+"A swapping (right) of an element was possible"); for(k = 1; k<=size(tempList);k++) { tempResult = insert(tempResult,result[i]); for (l=posInhomog; l<=j;l++) { tempResult[1] = delete(tempResult[1],posInhomog); } for (l = size(tempList[k]);l>=1;l--) { tempResult[1] = insert(tempResult[1],tempList[k][l],posInhomog-1); } } } }//Checking if we can swap left resp. right divisors result = result + tempResult; result = delete_dublicates_noteval(result); return(result); }//extractHomogeneousDivisorsNthWeyl static proc extractHomogeneousDivisorsLeft(poly h) "INPUT: A polynomial h in the first Weyl algebra OUTPUT: If h is homogeneous, then all factorizations of h are returned. If h is inhomogeneous, then a list l is returned whose entries are again lists k = [k_1,...,k_n], where k_1*...*k_n = h. The entry k_n is inhomogeneous and has no other homogeneous left divisors any more. All the other entries in k are homogeneous polynomials. " {//extractHomogeneousDivisorsLeft int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i;int j; int k; list result; poly hath; list recResult; map invo = basering,-var(1),var(2); for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (homogwithorder(h,intvec(-1,1))) {//given polynomial was homogeneous already dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations."); result = homogfacFirstWeyl_all(h); for (i = 1; i<=size(result);i++) {//removing the first entry (coefficient) from the list result result[i] = delete(result[i],1); }//removing the first entry (coefficient) from the list result return(result); }//given polynomial was homogeneous already list hlist = homogDistribution(h); dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands."); for (i = 1; i<= size(hlist); i++) { hlist[i] = homogfacFirstWeyl_all(hlist[i][2]); if (size(hlist[i][1])==1) {//One homogeneous part just has a trivial factorization if(hlist[i][1][1] == 0) { hlist = delete(hlist,i); continue; } else { return(list(list(h))); } }//One homogeneous part just has a trivial factorization } dbprint(p,dbprintWhitespace+ " Done."); dbprint(p,dbprintWhitespace+ " Trying to find Left divisors"); list alreadyConsideredCandidates; poly candidate; int isCandidate; for (i = 1; i<=size(hlist[1]);i++) {//Finding candidates for homogeneous left divisors of h candidate = hlist[1][i][2]; isCandidate = 0; for (j=1;j<=size(alreadyConsideredCandidates);j++) { if(alreadyConsideredCandidates[j] == candidate) { isCandidate =1; break; } } if(isCandidate) { i++; continue; } else { alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate); } dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous left divisor"); for (j = 2; j<=size(hlist);j++) {//Iterating through the other homogeneous parts isCandidate = 0; for(k=1; k<=size(hlist[j]);k++) { if(hlist[j][k][2]==candidate) { isCandidate = 1; break; } } if(!isCandidate) { break; } }//Iterating through the other homogeneous parts if(isCandidate) {//candidate was really a left divisor dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous left divisor"); hath = involution(lift(involution(candidate,invo),involution(h,invo))[1,1],invo); recResult = extractHomogeneousDivisorsLeft(hath); for (j = 1; j<=size(recResult); j++) { recResult[j] = insert(recResult[j],candidate); } result = result + recResult; }//Candidate was really a left divisor }//Finding candidates for homogeneous left divisors of h if (size(result)==0) { return(list(list(h))); } return(result); }//extractHomogeneousDivisorsLeft static proc extractHomogeneousDivisorsLeftNthWeyl(poly h) "INPUT: A polynomial h in the nth Weyl algebra OUTPUT: If h is homogeneous with respect to the ZZ grading on the nth Weyl algebra, then all factorizations of h are returned. If h is inhomogeneous, then a list l is returned whose entries are again lists k = [k_1,...,k_n], where k_1*...*k_n = h. The entry k_n is inhomogeneous and has no other homogeneous left divisors any more. All the other entries in k are homogeneous polynomials. GENERAL ASSUMPTIONS: - The basering is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//extractHomogeneousDivisorsLeftNthWeyl int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i;int j; int k; list result; poly hath; list recResult; ideal invoIdeal; for (i = 1; i <= nvars(basering); i++) {//Filling the mapping rules for the involution map if (i <= nvars(basering) div 2) {invoIdeal[i] = -var(i);} else {invoIdeal[i] = var(i);} }//Filling the mapping rules for the involution map map invo = basering,invoIdeal; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (homogwithorderNthWeyl(h)) {//given polynomial was homogeneous already dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations."); result = homogfacNthWeyl_all(h); for (i = 1; i<=size(result);i++) {//removing the first entry (coefficient) from the list result result[i] = delete(result[i],1); }//removing the first entry (coefficient) from the list result return(result); }//given polynomial was homogeneous already list hlist = homogDistributionNthWeyl(h); dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands."); for (i = 1; i<= size(hlist); i++) { hlist[i] = homogfacNthWeyl_all(hlist[i][2]); if (size(hlist[i][1])==1) {//One homogeneous part just has a trivial factorization if(hlist[i][1][1] == 0) { hlist = delete(hlist,i); continue; } else { return(list(list(h))); } }//One homogeneous part just has a trivial factorization } dbprint(p,dbprintWhitespace+ " Done."); dbprint(p,dbprintWhitespace+ " Trying to find Left divisors"); list alreadyConsideredCandidates; poly candidate; int isCandidate; for (i = 1; i<=size(hlist[1]);i++) {//Finding candidates for homogeneous left divisors of h candidate = hlist[1][i][2]; isCandidate = 0; for (j=1;j<=size(alreadyConsideredCandidates);j++) { if(alreadyConsideredCandidates[j] == candidate) { isCandidate =1; break; } } if(isCandidate) { i++; continue; } else { alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate); } dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous left divisor"); for (j = 2; j<=size(hlist);j++) {//Iterating through the other homogeneous parts isCandidate = 0; for(k=1; k<=size(hlist[j]);k++) { if(hlist[j][k][2]==candidate) { isCandidate = 1; break; } } if(!isCandidate) { break; } }//Iterating through the other homogeneous parts if(isCandidate) {//candidate was really a left divisor dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous left divisor"); hath = involution(lift(involution(candidate,invo),involution(h,invo))[1,1],invo); recResult = extractHomogeneousDivisorsLeftNthWeyl(hath); for (j = 1; j<=size(recResult); j++) { recResult[j] = insert(recResult[j],candidate); } result = result + recResult; }//Candidate was really a left divisor }//Finding candidates for homogeneous left divisors of h if (size(result)==0) { return(list(list(h))); } return(result); }//extractHomogeneousDivisorsLeftNthWeyl /* Interesting test cases: h=x1^2*x2*x3*d1*d2*d3+x1*x2^2*x3*d1*d2*d3+x1*x2*x3*d1^2*d2*d3+x1*x2*x3*d1*d3+x1*x2*x3*d2*d3 h = x1*(d1+x1); h = x1*d2 + x3*d3; h = x1*x2*d2^2+x2^2*d2*d3+x1*d2+2*x2*d3 */ static proc extractHomogeneousDivisorsRight(poly h) "INPUT: A polynomial h in the first Weyl algebra OUTPUT: If h is homogeneous, then all factorizations of h are returned. If h is inhomogeneous, then a list l is returned whose entries are again lists k = [k_1,...,k_n], where k_1*...*k_n = h. The entry k_1 is inhomogeneous and has no other homogeneous right divisors any more. All the other entries in k are homogeneous polynomials. " {//extractHomogeneousDivisorsRight int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i;int j; int k; list result; poly hath; list recResult; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (homogwithorder(h,intvec(-1,1))) {//given polynomial was homogeneous already dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations."); result = homogfacFirstWeyl_all(h); for (i = 1; i<=size(result);i++) {//removing the first entry (coefficient) from the list result result[i] = delete(result[i],1); }//removing the first entry (coefficient) from the list result return(result); }//given polynomial was homogeneous already list hlist = homogDistribution(h); dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands."); for (i = 1; i<= size(hlist); i++) { hlist[i] = homogfacFirstWeyl_all(hlist[i][2]); if (size(hlist[i][1])==1) {//One homogeneous part just has a trivial factorization if(hlist[i][1][1] == 0) { hlist = delete(hlist,i); continue; } else { return(list(list(h))); } }//One homogeneous part just has a trivial factorization } dbprint(p,dbprintWhitespace+ " Done."); dbprint(p,dbprintWhitespace+ " Trying to find right divisors"); list alreadyConsideredCandidates; poly candidate; int isCandidate; for (i = 1; i<=size(hlist[1]);i++) {//Finding candidates for homogeneous left divisors of h candidate = hlist[1][i][size(hlist[1][i])]; isCandidate = 0; for (j=1;j<=size(alreadyConsideredCandidates);j++) { if(alreadyConsideredCandidates[j] == candidate) { isCandidate =1; break; } } if(isCandidate) { i++; continue; } else { alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate); } dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous r-divisor"); for (j = 2; j<=size(hlist);j++) {//Iterating through the other homogeneous parts isCandidate = 0; for(k=1; k<=size(hlist[j]);k++) { if(hlist[j][k][size(hlist[j][k])]==candidate) { isCandidate = 1; break; } } if(!isCandidate) { break; } }//Iterating through the other homogeneous parts if(isCandidate) {//candidate was really a left divisor dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous right divisor"); hath = lift(candidate,h)[1,1]; recResult = extractHomogeneousDivisorsRight(hath); for (j = 1; j<=size(recResult); j++) { recResult[j] = insert(recResult[j],candidate,size(recResult[j])); } result = result + recResult; }//Candidate was really a left divisor }//Finding candidates for homogeneous left divisors of h if (size(result)==0) { result = list(list(h)); } return(result); }//extractHomogeneousDivisorsRight static proc extractHomogeneousDivisorsRightNthWeyl(poly h) "INPUT: A polynomial h in the nth Weyl algebra OUTPUT: If h is homogeneous with respect to the ZZ grading on the nth Weyl algebra, then all factorizations of h are returned. If h is inhomogeneous, then a list l is returned whose entries are again lists k = [k_1,...,k_n], where k_1*...*k_n = h. The entry k_1 is inhomogeneous and has no other homogeneous right divisors any more. All the other entries in k are homogeneous polynomials. GENERAL ASSUMPTIONS: - The basering is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//extractHomogeneousDivisorsRightNthWeyl int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i;int j; int k; list result; poly hath; list recResult; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (homogwithorderNthWeyl(h)) {//given polynomial was homogeneous already dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations."); result = homogfacNthWeyl_all(h); for (i = 1; i<=size(result);i++) {//removing the first entry (coefficient) from the list result result[i] = delete(result[i],1); }//removing the first entry (coefficient) from the list result return(result); }//given polynomial was homogeneous already list hlist = homogDistributionNthWeyl(h); dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands."); for (i = 1; i<= size(hlist); i++) { hlist[i] = homogfacNthWeyl_all(hlist[i][2]); if (size(hlist[i][1])==1) {//One homogeneous part just has a trivial factorization if(hlist[i][1][1] == 0) { hlist = delete(hlist,i); continue; } else { return(list(list(h))); } }//One homogeneous part just has a trivial factorization } dbprint(p,dbprintWhitespace+ " Done."); dbprint(p,dbprintWhitespace+ " Trying to find right divisors"); list alreadyConsideredCandidates; poly candidate; int isCandidate; for (i = 1; i<=size(hlist[1]);i++) {//Finding candidates for homogeneous left divisors of h candidate = hlist[1][i][size(hlist[1][i])]; isCandidate = 0; for (j=1;j<=size(alreadyConsideredCandidates);j++) { if(alreadyConsideredCandidates[j] == candidate) { isCandidate =1; break; } } if(isCandidate) { i++; continue; } else { alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate); } dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous r-divisor"); for (j = 2; j<=size(hlist);j++) {//Iterating through the other homogeneous parts isCandidate = 0; for(k=1; k<=size(hlist[j]);k++) { if(hlist[j][k][size(hlist[j][k])]==candidate) { isCandidate = 1; break; } } if(!isCandidate) { break; } }//Iterating through the other homogeneous parts if(isCandidate) {//candidate was really a left divisor dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous right divisor"); hath = lift(candidate,h)[1,1]; recResult = extractHomogeneousDivisorsRightNthWeyl(hath); for (j = 1; j<=size(recResult); j++) { recResult[j] = insert(recResult[j],candidate,size(recResult[j])); } result = result + recResult; }//Candidate was really a left divisor }//Finding candidates for homogeneous left divisors of h if (size(result)==0) { result = list(list(h)); } return(result); }//extractHomogeneousDivisorsRightNthWeyl static proc fromZeroHomogToThetaPoly(poly h) " //DEPRECATED INPUT: A polynomial h in the first Weyl algebra, homogeneous of degree 0 OUTPUT: The ring Ktheta with a polynomial result representing h as polynomial in theta ASSUMPTIONS: - h is homogeneous of degree 0 with respect to the [-1,1] weight vector - The basering is the first Weyl algebra " {//proc fromZeroHomogToThetaPoly int i; int j; list mons; if(!homogwithorder(h,intvec(-1,1))) {//Input not a homogeneous polynomial causes an error ERROR("The input was not a homogeneous polynomial"); }//Input not a homogeneous polynomial causes an error if(deg(h,intvec(-1,1))!=0) {//Input does not have degree 0 ERROR("The input did not have degree 0"); }//Input does not have degree 0 for(i = 1; i<=size(h);i++) {//Putting the monomials in a list mons = mons+list(h[i]); }//Putting the monomials in a list ring KTheta = 0,(x,y,theta),dp; setring KTheta; map thetamap = r,x,y; list mons = thetamap(mons); poly entry; for (i = 1; i<=size(mons);i++) {//transforming the monomials as monomials in theta entry = leadcoef(mons[i]); for (j = 0; j0) { tempList = getAllCombOfHomogFact(f1); for (i = 1; i<=size(tempList); i++) {//Every combination combined with the coefficient possibilities for (j = 1; j<=size(coeffTuplesMax); j++) {//iterating through the possible coefficient choices pqmax = pqmax + list(list(coeffTuplesMax[j][1]*tempList[i][1], coeffTuplesMax[j][2]*tempList[i][2])); }//iterating through the possible coefficient choices }//Every combination combined with the coefficient possibilities for (i = 1; i<=size(coeffTuplesMax); i++) { pqmax = pqmax + list(list(coeffTuplesMax[i][1],maxh/coeffTuplesMax[i][1])); pqmax = pqmax + list(list(maxh/coeffTuplesMax[i][2],coeffTuplesMax[i][2])); } } }//the maximal homogeneous factor is not a constant //Now we go to pqmin if (size(f2[1]) == 1) {//the minimal homogeneous factor is a constant pqmin = coeffTuplesMin; }//the minimal homogeneous factor is a constant else {//the minimal homogeneous factor is not a constant for (i = 1; i<=size(f2); i++) {//We can forget about the first coefficient now. Therefore we will delete him from the list. f2[i] = delete(f2[i],1); if(size(f2[i])==1) {//trivial thing for (j = 1; j<=size(coeffTuplesMin); j++) { pqmin = pqmin + list(list(coeffTuplesMin[j][1],coeffTuplesMin[j][2]*f2[i][1])); pqmin = pqmin + list(list(coeffTuplesMin[j][1]*f2[i][1],coeffTuplesMin[j][2])); } f2 = delete(f2,i); continue; } }//We can forget about the first coefficient now. Therefore we will delete him from the list. if(size(f2)>0) { tempList = getAllCombOfHomogFact(f2); for (i = 1; i<=size(tempList); i++) {//Every combination combined with the coefficient possibilities for (j = 1; j<=size(coeffTuplesMin); j++) {//iterating through the possible coefficient choices pqmin = pqmin + list(list(coeffTuplesMin[j][1]*tempList[i][1], coeffTuplesMin[j][2]*tempList[i][2])); }//iterating through the possible coefficient choices }//Every combination combined with the coefficient possibilities for (i = 1; i<=size(coeffTuplesMin); i++) { pqmin = pqmin + list(list(coeffTuplesMin[i][1],minh/coeffTuplesMin[i][1])); pqmin = pqmin + list(list(minh/coeffTuplesMin[i][2],coeffTuplesMin[i][2])); } } }//the minimal homogeneous factor is not a constant //and now we combine them together to obtain all possibilities. for (i = 1; i<=size(pqmax); i++) {//iterate over the maximal homogeneous combination possibilities for (j = 1; j<=size(pqmin); j++) {//iterate over the minimal homogeneous combiniation possibilities if (deg(pqmax[i][1], ivm11)>=deg(pqmin[j][1],ivm11) and deg(pqmax[i][2], ivm11)>=deg(pqmin[j][2],ivm11)) { if (pqmax[i][1]+pqmin[j][1]!=0 and pqmax[i][2]+pqmin[j][2]!=0) { if (deg(h,ivm11)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),ivm11)) { j++; continue; } if (deg(h,iv1m1)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),iv1m1)) { j++; continue; } result = result +list(list(pqmax[i][1]+pqmin[j][1],pqmax[i][2]+pqmin[j][2])); } } }//iterate over the minimal homogeneous combiniation possibilities }//iterate over the maximal homogeneous combination possibilities //Now deleting double entries result = delete_dublicates_noteval(result); return(result); }//proc computeCombinationsMinMaxHomog static proc computeCombinationsMinMaxHomogNthWeyl(poly h) "Input: A polynomial h in the nth Weyl Algebra Output: Combinations of the form (p_max + p_min)(q_max + q_min), such that p_max, p_min, q_max and q_min are homogeneous and p_max*q_max equals the maximal homogeneous part in h, and p_max * q_max equals the minimal homogeneous part in h. GENERAL ASSUMPTIONS: - h is not homogeneous. "{//proc computeCombinationsMinMaxHomogNthWeyl int p=printlevel-voice+2; // for dbprint string dbprintWhitespace = ""; int i; int j; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} list hList = homogDistributionNthWeyl(h); dbprint(p,dbprintWhitespace + "Computing maximal and minimal homogeneous part of h"); poly maxh = hList[size(hList)][2]; poly minh = hList[1][2]; dbprint(p,dbprintWhitespace + "Done. They are:"); dbprint(p,string(maxh) + ", " + string(minh)); dbprint(p, dbprintWhitespace + "Computing their respective homogeneous factorizations"); list f1 = homogfacNthWeyl_all(maxh); list f2 = homogfacNthWeyl_all(minh); dbprint(p,dbprintWhitespace + "Done."); list result = list(); list pqmax = list(); list pqmin = list(); list tempList = list(); //First, we are going to deal with our most hated guys: The Coefficients. // dbprint(p,dbprintWhitespace + "We get all combinations for the coefficient of the maximal homogeneous part"); list coeffTuplesMax = getAllCoeffTuplesComb(factorizeInt(number(f1[1][1]))); //We can assume without loss of generality, that p_max has a //nonnegative leading coefficient for (i = 1; i<=size(coeffTuplesMax);i++) {//Deleting all tuples with negative entries for p_max if (coeffTuplesMax[i][1]<0) { coeffTuplesMax = delete(coeffTuplesMax,i); continue; } }//Deleting all tuples with negative entries for p_max dbprint(p,dbprintWhitespace + "Done. The Combinations are:"); dbprint(p,coeffTuplesMax); dbprint(p,dbprintWhitespace + "We get all combinations for the coefficient of the minimal homogeneous part"); list coeffTuplesMin = getAllCoeffTuplesComb(factorizeInt(number(f2[1][1]))); dbprint(p,dbprintWhitespace + "Done. The Combinations are:"); dbprint(p,coeffTuplesMin); //Now, we will be actally dealing with the Combinations. //Let's start with the pqmax if (size(f1[1]) == 1) {//the maximal homogeneous factor is a constant dbprint(p,dbprintWhitespace + "the maximal homogeneous factor is a constant."); pqmax = coeffTuplesMax; }//the maximal homogeneous factor is a constant else {//the maximal homogeneous factor is not a constant dbprint(p,dbprintWhitespace + "Deleting the first list entry in each list (constant factor)."); for (i = 1; i<=size(f1); i++) {//We can forget about the first coefficient now. Therefore we will delete him from the list. f1[i] = delete(f1[i],1); if(size(f1[i])==1) {//trivial thing for (j = 1; j<=size(coeffTuplesMax); j++) { pqmax = pqmax + list(list(coeffTuplesMax[j][1],coeffTuplesMax[j][2]*f1[i][1])); pqmax = pqmax + list(list(coeffTuplesMax[j][1]*f1[i][1],coeffTuplesMax[j][2])); } f1 = delete(f1,i); continue; }//trivial thing }//We can forget about the first coefficient now. Therefore we will delete him from the list. dbprint(p,dbprintWhitespace + "Done."); dbprint(p,dbprintWhitespace + "Putting all possible pre-coefficients besides the entries in pqmax."); if (size(f1)>0) { tempList = getAllCombOfHomogFact(f1); for (i = 1; i<=size(tempList); i++) {//Every combination combined with the coefficient possibilities for (j = 1; j<=size(coeffTuplesMax); j++) {//iterating through the possible coefficient choices pqmax = pqmax + list(list(coeffTuplesMax[j][1]*tempList[i][1], coeffTuplesMax[j][2]*tempList[i][2])); }//iterating through the possible coefficient choices }//Every combination combined with the coefficient possibilities for (i = 1; i<=size(coeffTuplesMax); i++) { pqmax = pqmax + list(list(coeffTuplesMax[i][1],maxh/coeffTuplesMax[i][1])); pqmax = pqmax + list(list(maxh/coeffTuplesMax[i][2],coeffTuplesMax[i][2])); } } dbprint(p,dbprintWhitespace + "Done."); }//the maximal homogeneous factor is not a constant dbprint(p, dbprintWhitespace + "Doing the same for f2"); //Now we go to pqmin if (size(f2[1]) == 1) {//the minimal homogeneous factor is a constant pqmin = coeffTuplesMin; }//the minimal homogeneous factor is a constant else {//the minimal homogeneous factor is not a constant for (i = 1; i<=size(f2); i++) {//We can forget about the first coefficient now. Therefore we will delete him from the list. f2[i] = delete(f2[i],1); if(size(f2[i])==1) {//trivial thing for (j = 1; j<=size(coeffTuplesMin); j++) { pqmin = pqmin + list(list(coeffTuplesMin[j][1],coeffTuplesMin[j][2]*f2[i][1])); pqmin = pqmin + list(list(coeffTuplesMin[j][1]*f2[i][1],coeffTuplesMin[j][2])); } f2 = delete(f2,i); continue; } }//We can forget about the first coefficient now. Therefore we will delete him from the list. if(size(f2)>0) { tempList = getAllCombOfHomogFact(f2); for (i = 1; i<=size(tempList); i++) {//Every combination combined with the coefficient possibilities for (j = 1; j<=size(coeffTuplesMin); j++) {//iterating through the possible coefficient choices pqmin = pqmin + list(list(coeffTuplesMin[j][1]*tempList[i][1], coeffTuplesMin[j][2]*tempList[i][2])); }//iterating through the possible coefficient choices }//Every combination combined with the coefficient possibilities for (i = 1; i<=size(coeffTuplesMin); i++) { pqmin = pqmin + list(list(coeffTuplesMin[i][1],minh/coeffTuplesMin[i][1])); pqmin = pqmin + list(list(minh/coeffTuplesMin[i][2],coeffTuplesMin[i][2])); } } }//the minimal homogeneous factor is not a constant dbprint(p,dbprintWhitespace + "Done."); //and now we combine them together to obtain all possibilities. for (i = 1; i<=size(pqmax); i++) {//iterate over the maximal homogeneous combination possibilities for (j = 1; j<=size(pqmin); j++) {//iterate over the minimal homogeneous combiniation possibilities if (degreeOfNthWeylPoly(pqmax[i][1])>=degreeOfNthWeylPoly(pqmin[j][1]) and degreeOfNthWeylPoly(pqmax[i][2])>=degreeOfNthWeylPoly(pqmin[j][2])) { if (pqmax[i][1]+pqmin[j][1]!=0 and pqmax[i][2]+pqmin[j][2]!=0) { if (h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2])==0) { /*REMARK: This part here is different from factoring the first Weyl algebra. As the output of degreeOfNthWeylPoly is in general an intvec instead of an int, we cannot say deg(0)=-1; Therefore, we need to catch the 0-case at this point, i.e. min and max homog are combined equal to h. */ result = result +list(list(pqmax[i][1]+pqmin[j][1],pqmax[i][2]+pqmin[j][2])); j++;continue; } if (degreeOfNthWeylPoly(h)<= degreeOfNthWeylPoly(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]))) { j++; continue; } if (degreeOfNthWeylPolyInverted(h)<= degreeOfNthWeylPolyInverted(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]))) { j++; continue; } result = result +list(list(pqmax[i][1]+pqmin[j][1],pqmax[i][2]+pqmin[j][2])); } } }//iterate over the minimal homogeneous combiniation possibilities }//iterate over the maximal homogeneous combination possibilities //Now deleting double entries result = delete_dublicates_noteval(result); return(result); }//proc computeCombinationsMinMaxHomogNthWeyl static proc getAllCombOfHomogFact(list l) "Gets called in computeCombinationsMinMaxHomog. It gets a list of different homogeneous factorizations of one homogeneous polynomial and returns the possibilities to combine them into two factors. Assumptions: - The list does not contain the first coefficient. - The list contains at least one list with two elements." {//proc getAllCombOfHomogFact list result; list leftAndRightHandSides; int i; int j; list tempset; if (size(l)==1 and size(l[1])==2) { result = result + list(list(l[1][1],l[1][2])); return(result); } leftAndRightHandSides = getPossibilitiesForRightSides(l); for (i = 1; i<=size(leftAndRightHandSides); i++) { result =result+list(list(leftAndRightHandSides[i][1],product(leftAndRightHandSides[i][2][1]))); //tidy up the right hand sides, because, if it is just one irreducible factor, we are done for (j = 1; j<=size(leftAndRightHandSides[i][2]);j++) {//Tidy up right hand sides if (size(leftAndRightHandSides[i][2][j])<2) {//Element can be dismissed leftAndRightHandSides[i][2] = delete(leftAndRightHandSides[i][2],j); continue; }//Element can be dismissed }//Tidy up right hand sides if (size(leftAndRightHandSides[i][2])>0) { tempset = getAllCombOfHomogFact(leftAndRightHandSides[i][2]); for (j = 1; j<=size(tempset);j++) {//multiplying the first factor with the left hand side result = result + list(list(leftAndRightHandSides[i][1]*tempset[j][1],tempset[j][2])); }//multiplying the first factor with the left hand side } } return(result); }//proc getAllCombOfHomogFact static proc getPossibilitiesForRightSides(list l) "Given a list of different factorizations l, this function returns a list of the form (a,{(a_2,...,a_n)| (a,a_2,...,a_n) in A})" {//getPossibilitiesForRightSide list templ = l; list result; poly firstElement; list rightSides; list tempRightSide; int i; int j; while (size(templ)>0) { firstElement = templ[1][1]; rightSides = list(); for (i = 1; i<= size(templ); i++) { if (templ[i][1] == firstElement) {//save the right sides tempRightSide = list(); for (j = 2; j<=size(templ[i]);j++) { tempRightSide = tempRightSide + list(templ[i][j]); } if (size(tempRightSide)!=0) { rightSides = rightSides + list(tempRightSide); } templ = delete(templ,i); continue; }//save the right sides } result = result + list(list(firstElement,rightSides)); } return(result); }//getPossibilitiesForRightSide static proc getAllCoeffTuplesComb(list l)" Given the output of factorizeInt ((a_1,...,a_n),(i_1,...,i_n)) , it returns all possible tuples of the set {(a,b) | There exists an real N!=emptyset subset of {1,...,n}, such that a = prod_{i \in N}a_i, b=prod_{i \not\in N} a_i} Assumption: The list is sorted from smallest integer to highest. - it is not the factorization of 0. " {//proc getAllCoeffTuplesComb list result; if (l[1][1] == 0) { ERROR("getAllCoeffTuplesComb: Zero Coefficients as leading and Tail Coeffs? That is not possible. Something went wrong."); } if (size(l[1]) == 1) {//Trivial Factorization, just 1 if (l[1][1] == 1) { return(list(list(1,1),list(-1,-1))); } else { return(list(list(-1,1),list(1,-1))); } }//Trivial Factorization, just 1 if (size(l[1]) == 2 and l[2][2]==1) {//Just a prime number if (l[1][1] == 1) { result = list(list(l[1][2],1),list(1,l[1][2])); result = result + list(list(-l[1][2],-1),list(-1,-l[1][2])); return(result); } else { result = list(list(l[1][2],-1),list(1,-l[1][2])); result = result + list(list(-l[1][2],1),list(-1,l[1][2])); return(result); } }//Just a prime number //Now comes the interesting case: a product of primes list tempPrimeFactors; list tempPowersOfThem; int i; for (i = 2; i<=size(l[1]);i++) {//Removing the starting 1 or -1 to get the N's tempPrimeFactors[i-1] = l[1][i]; tempPowersOfThem[i-1] = l[2][i]; }//Removing the starting 1 or -1 to get the N's list Ns = getAllSubsetsN(list(tempPrimeFactors,tempPowersOfThem)); list tempTuples; number productOfl = multiplyFactIntOutput(l); if (productOfl<0){productOfl = -productOfl;} tempTuples = tempTuples + list(list(1,productOfl),list(productOfl,1)); for (i = 1; i<=size(Ns); i++) { if (productOfl/Ns[i]>Ns[i]) {//TODO: BEWEISEN, dass das die einzigen Combos sind tempTuples = tempTuples + list(list(Ns[i],productOfl/Ns[i]),list(productOfl/Ns[i],Ns[i])); }//TODO: BEWEISEN, dass das die einzigen Combos sind if (productOfl/Ns[i]==Ns[i]) { tempTuples = tempTuples + list(list(Ns[i],Ns[i])); } } //And now, it just remains to get the -1s and 1-s correctly to the tuples list tempEntry; if (l[1][1] == 1) { for (i = 1; i<=size(tempTuples);i++) {//Adding everything to result tempEntry = tempTuples[i]; result = result + list(tempEntry); result = result + list(list(-tempEntry[1], -tempEntry[2])); }//Adding everyThing to Result } else { for (i = 1; i<=size(tempTuples);i++) {//Adding everything to result tempEntry = tempTuples[i]; result = result + list(list(tempEntry[1],-tempEntry[2])); result = result + list(list(-tempEntry[1], tempEntry[2])); }//Adding everyThing to Result } return(result); }//proc getAllCoeffTuplesComb static proc contains(list l, int elem) "Assumption: l is sorted" {//Binary Search in list if (size(l)<=1) { if(size(l) == 0){return(0);} if (l[1]!=elem){return(0);} else{return(1);} } int imax = size(l); int imin = 1; int imid; while(imax >= imin) { imid = (imin + imax)/2; if (l[imid] == elem){return(1);} if (l[imid] =1; i--) { result = insert(result, l[i]); } return(result); }//proc fromIntvecToList static proc factorizeInt(number n) "Given an integer n, factorizeInt computes its factorization. The output is a list containing two intvecs. The first contains the prime factors, the second its powers. ASSUMPTIONS: - n is given as integer number "{ if (n==0) {return(list(list(0),list(1)));} int i; list temp = primefactors(n); if (n<0) {list result = list(list(-1),list(1));} else {list result = list(list(1),list(1));} result[1] = result[1] + temp[1]; result[2] = result[2] + temp[2]; return(result); } static proc homogDistribution(poly h) "Input: A polynomial in the first Weyl Algebra. Output: A two-dimensional list of the following form. Every sublist contains exactly two entries. One for the Z-degree of the corresponding homogeneous part (integer), and the homogeneous polynomial itself, and those sublists are oredered by ascending degree. For example a call of homogDistribution(x+d+1) would have the output [1]: [1]: -1 [2]: x [2]: [1]: 0 [2]: 1 [3]: [1]: 1 [2]: d "{//homogDistribution if (h == 0) {//trivial case where input is 0 return(list(list(0,0))); }//trivial case where input is 0 if (!isWeyl()) {//Our basering is not the Weyl algebra ERROR("Ring was not the first Weyl algebra"); return(list()); }//Our basering is not the Weyl algebra if(nvars(basering)!=2) {//Our basering is the Weyl algebra, but not the first ERROR("Ring is not the first Weyl algebra"); return(list()); }//Our basering is the Weyl algebra, but not the first intvec ivm11 = intvec(-1,1); intvec iv1m1 = intvec(1,-1); poly tempH = h; poly minh; list result = list(); int nextExpectedDegree = -deg(tempH,iv1m1); while (tempH != 0) { minh = jet(tempH,deg(tempH,iv1m1),iv1m1)-jet(tempH,deg(tempH,iv1m1)-1,iv1m1); while (deg(minh,ivm11)>nextExpectedDegree) {//filling empty homogeneous spaces with 0 result = result + list(list(nextExpectedDegree,0)); nextExpectedDegree = nextExpectedDegree +1; }//filling empty homogeneous spaces with 0 result = result + list(list(deg(minh,ivm11),minh)); tempH = tempH - minh; nextExpectedDegree = nextExpectedDegree +1; } return(result); }//homogDistribution static proc homogDistributionNthWeyl(poly h) " INPUT: A polynomial in the n-th Weyl algebra OUTPUT: A two-dimensional list of the following form. Every sublist contains exactly two entries. One for the Z^n-degree of the corresponding homogeneous part (intvec), and the homogeneous polynomial itself, and those sublists are oredered by ascending degree using lexicographical ordering on Z^n. Different from homogDistribution, the 0-summands between the maximum and minimum homogeneous degree are not displayed. For example a call of homogDistribution(x1+d2+1) would have the output (ring is the second weyl algebra with variables x1,x2,d1,d2). [1]: [1]: [0,1] [2]: d2 [2]: [1]: [0,0] [2]: 1 [3]: [1]: [-1,0] [2]: x1 GENERAL ASSUMPTIONS: - The basering is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//proc homogDistributionNthWeyl if (h == 0) {//trivial case where input is 0 return(list(list(0,0))); }//trivial case where input is 0 //TODO: BUG in nctools? /* if (!isWeyl()) */ /* {//Our basering is not the Weyl algebra */ /* ERROR("Ring was not a Weyl algebra"); */ /* return(list()); */ /* }//Our basering is not the Weyl algebra */ list result; poly tempH = h; intvec degVec; poly leadPoly; while(tempH != 0) {//tempH is not equal to zero ==> We have still unconsidered homogeneous summands leadPoly = extractLeadingTermOfNthWeylPoly(tempH); degVec = degreeOfNthWeylPoly(tempH); result = insert(result,list(degVec,leadPoly)); tempH = tempH - leadPoly; }//tempH is not equal to zero ==> We have still unconsidered homogeneous summands return(result); }//proc homogDistributionNthWeyl static proc extractLeadingTermOfNthWeylPoly(poly h) " INPUT: A polynomial h in the nth Weyl algebra. OUTPUT: A polynomial p representing the homogeneous leading polynomial of h with respect to the -1,1 grading on the polynomial nth weyl algebra. GENERAL ASSUMPTIONS: - The ring given is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//extractLeadingTermOfNthWeylPoly poly result = 0; intvec leadDeg = degreeOfNthWeylPoly(h); int i; int j; int isPart = 0; intvec lExp; for (i = 1; i <=size(h); i++) {//iterating through the terms of h isPart = 1; lExp = leadexp(h[i]); for (j = 1; j<=nvars(basering) div 2; j++) {//checking if the term is part of the leading polynomial if(lExp[j + nvars(basering) div 2] - lExp[j] != leadDeg[j]) {//Summand was not part of leading polynomial isPart = 0; break; }//Summand was not part of leading polynomial }//checking if the term is part of the leading degree if (!isPart) {i++; continue;} else {//In this case, h[i] was part of the leading polynomial result = result + h[i]; }//In this case, h[i] was part of the leading polynomial }//iterating through the terms of h return(result); }//extractLeadingTermOfNthWeylPoly /* Test cases for this: ring R = 0,(x1,x2,x3,d1,d2,d3),dp; matrix C[6][6] = 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1; matrix D[6][6] = 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1, -1,0,0,0,0,0, 0,-1,0,0,0,0, 0,0,-1,0,0,0; def r = nc_algebra(C,D); setring(r); poly h = 0; extractLeadingTermOfNthWeylPoly(h); ==> 0 h = x3; extractLeadingTermOfNthWeylPoly(h); ==> x3 h = x1*x2*d1*d2 + x1^4*x2^2*x3^3*d1^4*d2^2*d3^3 + 1 + x1 + x2 + x3; extractLeadingTermOfNthWeylPoly(h); ==> x1^4*x2^2*x3^3*d1^4*d2^2*d3^3+x1*x2*d1*d2+1 h = (x1*d1)^5; extractLeadingTermOfNthWeylPoly(h); ==>x1^5*d1^5+10*x1^4*d1^4+25*x1^3*d1^3+15*x1^2*d1^2+x1*d1 */ static proc degreeOfNthWeylPoly(poly h) " INPUT: A polynomial h in the nth Weyl algebra OUTPUT: An intvector of size n, representing the degree of h with respect to the lexicographical ordering on ZZ^n by considering the -1,1 grading on the nth Weyl algebra. GENERAL ASSUMPTIONS: - The ring given is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//degreeOfNthWeylPoly intvec result = 0:(nvars(basering) div 2); int i; int j; int k; int alreadySmaller; intvec expVec; intvec posUpdated= 0:(nvars(basering) div 2); for(i = 1; i<=nvars(basering) div 2; i++) {//Going through every variable to determine the max. term for (j = 1; j <= size(h); j++) {//Iterating over every term in h alreadySmaller = 0; for (k = 1; k 0,0,0 h = x1 + d1 + x1*d2 + d1*d2; degreeOfNthWeylPoly(h); ==> 1,1,0 h = x1^5 + d1 + x1*d2 + d1*d2; degreeOfNthWeylPoly(h); ==> 1,1,0 h = 1; degreeOfNthWeylPoly(h); ==> 0,0,0 h = x1*d1; degreeOfNthWeylPoly(h); ==> 0,0,0 h = x2*d2; degreeOfNthWeylPoly(h); ==> 0,0,0 h = x3*d3; degreeOfNthWeylPoly(h); ==> 0,0,0 */ static proc degreeOfNthWeylPolyInverted(poly h) " INPUT: A polynomial h in the nth Weyl algebra OUTPUT: An intvector of size n, representing the degree of h with respect to the lexicographical ordering on ZZ^n by considering the 1,-1 grading on the nth Weyl algebra. GENERAL ASSUMPTIONS: - The ring given is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. " {//degreeOfNthWeylPoly intvec result = 0:(nvars(basering) div 2); int i; int j; int k; int alreadySmaller; intvec expVec; intvec posUpdated= 0:(nvars(basering) div 2); for(i = 1; i<=nvars(basering) div 2; i++) {//Going through every variable to determine the max. term for (j = 1; j <= size(h); j++) {//Iterating over every term in h alreadySmaller = 0; for (k = 1; k2) {//We want max. two optional arguments dbprint(p,dbprintWhitespace + " More than two optional arguments"); return(list()); }//We want max. two optional arguments dbprint(p,dbprintWhitespace + " Done"); list result; int j; if (size(#)==0) {//No optional argument is given dbprint(p,dbprintWhitespace + " No optional arguments"); int valid = 1; for (i = size(l);i>=1;i--) {//iterate over the elements of the given list if (size(result)>0) { if (product(l[i])!=result[size(l)-i]) { valid = 0; break; } } result = insert(result, product(l[i])); }//iterate over the elements of the given list return(valid); }//No optional argument is given else { dbprint(p,dbprintWhitespace + " Optional arguments are given."); int valid = 1; for (i = size(l);i>=1;i--) {//iterate over the elements of the given list if (product(l[i])!=#[1]) { valid = 0; } result = insert(result, product(l[i])-#[1]); }//iterate over the elements of the given list if(size(#)==2) { dbprint(p,dbprintWhitespace + " A third argument is given. Output is a list now."); return(result); } return(valid); } }//proc testfac example { "EXAMPLE:";echo=2; ring r = 0,(x,y),dp; def R = nc_algebra(1,1); setring R; poly h = (x^2*y^2+1)*(x^2); def t1 = facFirstWeyl(h); //fist a correct list testNCfac(t1); //now a correct list with the factorized polynomial testNCfac(t1,h); //now we put in an incorrect list without a polynomial t1[3][3] = y; testNCfac(t1); // take h as additional input testNCfac(t1,h); // take h as additional input and output list of differences testNCfac(t1,h,1); } //================================================== //Procedure facSubWeyl: //This procedure serves the purpose to compute a //factorization of a given polynomial in a ring, whose subring //is the first Weyl algebra. The polynomial must only contain //the two arguments, which are also given by the user. proc facSubWeyl(poly h, X, D) "USAGE: facSubWeyl(h,x,y); h, X, D polynomials RETURN: list ASSUME: X and D are variables of a basering, which satisfy DX = XD +1. @* That is, they generate the copy of the first Weyl algebra in a basering. @* Moreover, h is a polynomial in X and D only. PURPOSE: compute factorizations of the polynomial, which depends on X and D. EXAMPLE: example facSubWeyl; shows examples SEE ALSO: facFirstWeyl, testNCfac, facFirstShift "{ int p = printlevel - voice + 2; dbprint(p," Start initial Checks of the input."); // basering can be anything having a Weyl algebra as subalgebra def @r = basering; //We begin to check the input for assumptions // which are: X,D are vars of the basering, if ( (isVar(X)!=1) || (isVar(D)!=1) || (size(X)>1) || (size(D)>1) || (leadcoef(X) != number(1)) || (leadcoef(D) != number(1)) ) { ERROR("expected pure variables as generators of a subalgebra"); } // Weyl algebra: poly w = D*X-X*D-1; // [D,X]=1 poly u = D*X-X*D+1; // [X,D]=1 if (u*w!=0) { // that is no combination gives Weyl ERROR("2nd and 3rd argument do not generate a Weyl algebra"); } // one of two is correct int isReverted = 0; // Reverted Weyl if dx=xd-1 holds if (u==0) { isReverted = 1; } // else: do nothing // DONE with assumptions, Input successfully checked dbprint(p," Successful"); intvec lexpofX = leadexp(X); intvec lexpofD = leadexp(D); int varnumX=1; int varnumD=1; while(lexpofX[varnumX] != 1) { varnumX++; } while(lexpofD[varnumD] != 1) { varnumD++; } /* VL : to add printlevel stuff */ dbprint(p," Change positions of the two variables in the list, if needed"); if (isReverted) { ring firstweyl = 0,(var(varnumD),var(varnumX)),dp; def Firstweyl = nc_algebra(1,1); setring Firstweyl; ideal M = 0:nvars(@r); M[varnumX]=var(2); M[varnumD]=var(1); map Q = @r,M; poly h= Q(h); } else { // that is unReverted ring firstweyl = 0,(var(varnumX),var(varnumD)),dp; def Firstweyl = nc_algebra(1,1); setring Firstweyl; poly h= imap(@r,h); } dbprint(p," Done!"); list result = facFirstWeyl(h); setring @r; list result; if (isReverted) { // map swap back ideal M; M[1] = var(varnumD); M[2] = var(varnumX); map S = Firstweyl, M; result = S(result); } else { // that is unReverted result = imap(Firstweyl,result); } return(result); }//proc facSubWeyl example { "EXAMPLE:";echo=2; ring r = 0,(x,y,z),dp; matrix D[3][3]; D[1,3]=-1; def R = nc_algebra(1,D); // x,z generate Weyl subalgebra setring R; poly h = (x^2*z^2+x)*x; list fact1 = facSubWeyl(h,x,z); // compare with facFirstWeyl: ring s = 0,(z,x),dp; def S = nc_algebra(1,1); setring S; poly h = (x^2*z^2+x)*x; list fact2 = facFirstWeyl(h); map F = R,x,0,z; list fact1 = F(fact1); // it is identical to list fact2 testNCfac(fact1); // check the correctness again } //================================================== //================================================== //************From here: Shift-Algebra************** //================================================== //==================================================* //one factorization of a homogeneous polynomial //in the first Shift Algebra static proc homogfacFirstShift(poly h) {//proc homogfacFirstShift int p=printlevel-voice+2; //for dbprint def r = basering; poly hath; intvec iv01 = intvec(0,1); int i; int j; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (!homogwithorder(h,iv01)) {//The given polynomial is not homogeneous ERROR("The given polynomial is not homogeneous."); return(list()); }//The given polynomial is not homogeneous if (h==0) { return(list(0)); } list result; int m = deg(h,iv01); dbprint(p,dbprintWhitespace+" exclude the homogeneous part of deg. 0"); if (m>0) {//The degree is not zero hath = lift(var(2)^m,h)[1,1]; for (i = 1; i<=m;i++) { result = result + list(var(2)); } }//The degree is not zero else {//The degree is zero hath = h; }//The degree is zero ring tempRing = 0,(x),dp; setring tempRing; map thetamap = r,x,1; poly hath = thetamap(hath); dbprint(p,dbprintWhitespace+" Factorize it using commutative factorization."); list azeroresult = factorize(hath); list azeroresult_return_form; for (i = 1; i<=size(azeroresult[1]);i++) {//rewrite the result of the commutative factorization for (j = 1; j <= azeroresult[2][i];j++) { azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]); } }//rewrite the result of the commutative factorization setring(r); map finalmap = tempRing,var(1); list tempresult = finalmap(azeroresult_return_form); result = tempresult+result; return(result); }//proc homogfacFirstShift //================================================== //Computes all possible homogeneous factorizations static proc homogfacFirstShift_all(poly h) {//proc HomogfacFirstShiftAll int p=printlevel-voice+2; //for dbprint intvec iv11 = intvec(1,1); if (deg(h,iv11) <= 0 ) {//h is a constant return(list(list(h))); }//h is a constant def r = basering; list one_hom_fac; //stands for one homogeneous factorization int i; int j; int k; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} int shiftcounter; //Compute again a homogeneous factorization dbprint(p,dbprintWhitespace+" Computing one homog. factorization of the polynomial"); one_hom_fac = homogfacFirstShift(h); one_hom_fac = delete(one_hom_fac,1); if (size(one_hom_fac) == 0) {//there is no homogeneous factorization or the polynomial was not homogeneous return(list()); }//there is no homogeneous factorization or the polynomial was not homogeneous dbprint(p,dbprintWhitespace+" Permuting the 0-homogeneous part with the s"); list result = permpp(one_hom_fac); for (i = 1; i<=size(result);i++) { shiftcounter = 0; for (j = 1; j<=size(result[i]); j++) { if (result[i][j]==var(2)) { shiftcounter++; } else { result[i][j] = subst(result[i][j], var(1), var(1)-shiftcounter); } } result[i] = insert(result[i],1); } dbprint(p,dbprintWhitespace+" Deleting double entries in the resulting list"); result = delete_dublicates_noteval(result); return(result); }//proc HomogfacFirstShiftAll //================================================== //factorization of the first Shift Algebra static proc facFirstShift_old(poly h) "USAGE: facFirstShift(h); h a polynomial in the first shift algebra RETURN: list PURPOSE: compute all factorizations of a polynomial in the first shift algebra THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle ASSUME: basering is the first shift algebra NOTE: Every entry of the output list is a list with factors for one possible factorization. EXAMPLE: example facFirstShift; shows examples SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl "{//facFirstShift_old int p = printlevel - voice + 2; int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} dbprint(p,dbprintWhitespace +" Checking the input."); if(nvars(basering)!=2) {//Our basering is the Shift algebra, but not the first ERROR("Basering is not the first shift algebra"); return(list()); }//Our basering is the Shift algebra, but not the first def r = basering; setring r; list LR = ringlist(r); number @n = leadcoef(LR[5][1,2]); poly @p = LR[6][1,2]; if ( @n!=number(1) ) { ERROR("Basering is not the first shift algebra"); return(list()); } dbprint(p,dbprintWhitespace +" Done"); list result = list(); int j; int k; int l; //counter // create a ring with the ordering which makes shift algebra // graded // def r = basering; // done before ring tempRing = LR[1][1],(x,s),(a(0,1),Dp); def tempRingnc = nc_algebra(1,s); setring r; // information on relations if (@p == -var(1)) // reverted shift algebra { dbprint(p,dbprintWhitespace +" Reverted shift algebra. Swaping variables in Ringlist"); setring(tempRingnc); map transf = r, var(2), var(1); setring(r); map transfback = tempRingnc, var(2),var(1); // result = transfback(resulttemp); } else { if ( @p == var(2)) // usual shift algebra { setring(tempRingnc); map transf = r, var(1), var(2); // result = facshift(h); setring(r); map transfback = tempRingnc, var(1),var(2); } else { ERROR("Basering is not the first shift algebra"); return(list()); } } // main calls setring(tempRingnc); dbprint(p,dbprintWhitespace +" Factorize the given polynomial with the subroutine sFacShift"); list resulttemp = sFacShift(transf(h)); dbprint(p,dbprintWhitespace +" Successful"); setring(r); result = transfback(resulttemp); return( delete_dublicates_noteval(result) ); }//facFirstShift_old proc facFirstShift(poly h) "USAGE: facFirstShift(h); h a polynomial in the first shift algebra RETURN: list PURPOSE: compute all factorizations of a polynomial in the first shift algebra THEORY: This function is a wrapper for facShift. It exists to make this library downward-compatible with older versions. ASSUME: basering is the first shift algebra NOTE: Every entry of the output list is a list with factors for one possible factorization. EXAMPLE: example facFirstShift; shows examples SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl "{//facFirstShift return(facShift(h)); }//facFirstShift example { "EXAMPLE:";echo=2; ring R = 0,(x,s),dp; def r = nc_algebra(1,s); setring(r); poly h = (s^2*x+x)*s; facFirstShift(h); } proc facShift(poly h) "USAGE: facShift(h); h a polynomial in the n'th shift algebra RETURN: list PURPOSE: compute all factorizations of a polynomial in the nth shift algebra THEORY: h is mapped to the $n$th Weyl algebra and then factorized there. The factorizations are mapped back (S_n in subalgebra of Weyl algebra). ASSUME: basering is the nth shift algebra NOTE: Every entry of the output list is a list with factors for one possible factorization. EXAMPLE: example facFirstShift; shows examples SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl "{//facShift //Definition of printlevel variable int p = printlevel-voice+2; int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} dbprint(p,dbprintWhitespace +" Checking if the given algebra is a Shift algebra"); //Redefine the ring in my standard form if (!ncfactor_isShift()) {//Our basering is not the shift algebra ERROR("Ring was not a Shift algebra"); return(list()); }//Our basering is not the Shift algebra dbprint(p,dbprintWhitespace +" Successful"); //A last check before we start the real business: Is h maybe just //dependable on commutative variables? if (isInCommutativeSubRing(h)) {//h is in a commutative subring list hdepvars; intvec tempIntVec; for (i = 1; i<=nvars(basering) ; i++) { tempIntVec = 0:nvars(basering); tempIntVec[i] = 1; if (deg(h,tempIntVec)>0) { hdepvars = hdepvars + list(var(i)); } } if (size(hdepvars) ==0) {//We just have a constant return(list(list(h))); }//We just have a constant dbprint(p,dbprintWhitespace+"Polynomial was given commutative subring. Performing commutative factorization."); def r = basering; def rList = ringlist(basering); rList = delete(rList,5); rList = delete(rList,5); def tempRing = ring(rList); setring(tempRing); poly h = imap(r,h); list tempResult = factorize(h); list result = list(list()); int j; for (i = 1; i<=size(tempResult[1]); i++) { for (j = 1; j<=tempResult[2][i]; j++) { result[1] = result[1] + list(tempResult[1][i]); } } //mapping back setring(r); def result = imap(tempRing,result); dbprint(p,dbprintWhitespace+"result:"); dbprint(p,result); dbprint(p,dbprintWhitespace+"Computing all permutations of this factorization"); poly constantFactor = result[1][1]; result[1] = delete(result[1],1);//Deleting the constant factor result=permpp(result[1]); for (i = 1; i<=size(result);i++) {//Insert constant factor result[i] = insert(result[i],constantFactor); }//Insert constant factor dbprint(p,dbprintWhitespace+"Done."); return(result); }//h is in a commutative subring dbprint(p,dbprintWhitespace +" Successful"); list result = list(); int j; int k; int l; //counter if (!checkIfProperNthShift()) {//The given ring was not a proper nth shift algebra dbprint(p,dbprintWhitespace +" positions of the variables have to be switched"); dbprint(p,dbprintWhitespace + "Constructing the proper ring."); def r = basering; list tempRingList = ringlist(r); tempRingList = delete(tempRingList,6); list the_vars; for (i = 1; i<=nvars(r); i++) {the_vars[i] = var(i);} int maybeDInWrongPos; poly tempVariable; for (i = 1; i<=size(the_vars) div 2; i++) {//Swapping the variables as needed maybeDInWrongPos = 1; if (the_vars[i + size(the_vars) div 2]*the_vars[i] -the_vars[i]*the_vars[i + size(the_vars) div 2] == the_vars[i + size(the_vars) div 2]) { i++; continue; } //If we enter this line, there is a break with our property //condition for (j = i+1; j<=size(the_vars); j++) { if (the_vars[j]*the_vars[i]-the_vars[i]*the_vars[j]==the_vars[j]) {//In this case, we matched a var x to a repective s tempVariable = the_vars[i + size(the_vars) div 2]; the_vars[i + size(the_vars) div 2] = the_vars[j]; the_vars[j] = tempVariable; maybeDInWrongPos = 0; break; }//In this case, we matched a var x to a repective s } if (maybeDInWrongPos) {//var(i) is actually a s, not an x print("i has to be pushed to the end."); tempVariable = the_vars[i]; the_vars = delete(the_vars, i); the_vars = the_vars + list(tempVariable); continue; }//var(i) is actually a s, not an x }//Swapping the variables as needed for (i = 1; i<=size(the_vars); i++) {tempRingList[2][i] = string(the_vars[i]);} matrix DTemp[nvars(r)][nvars(r)]; for (i = 1; i<=ncols(DTemp) div 2; i++) { DTemp[i,i + nvars(r) div 2] = the_vars[i + nvars(r) div 2]; } tempRingList = tempRingList + list(DTemp); def tempRing = ring(tempRingList); dbprint(p,dbprintWhitespace + "Done. The altered ring is the following:"); dbprint(p,tempRing); setring(tempRing); //We have to go through an intermediate ring, as there is some strange //behaviour in Singular concerning the correctness of the matrix D. // See http://www.singular.uni-kl.de:8002/trac/ticket/542 for details. matrix DTemp = imap(r, DTemp); list tempRingList2 = ringlist(tempRing); tempRingList2[6]= DTemp; def tempRing2 = ring(tempRingList2); setring(tempRing2); poly h = imap(r,h); dbprint(p,dbprintWhitespace +" Successful"); list resulttemp = facShift(h); setring(r); result = imap(tempRing2,resulttemp); return (result); }//The given ring was not a proper nth Shift algebra dbprint(p, dbprintWhitespace +" factorization of the polynomial with the routine sfacNthShift"); result = sFacNthShift(h); dbprint(p,dbprintWhitespace +" Done"); return(result); }//facShift example { "EXAMPLE:";echo=2; ring R = 0,(x1,x2,s1,s2),dp; matrix C[4][4] = 1,1,1,1, 1,1,1,1, 1,1,1,1, 1,1,1,1; matrix D[4][4] = 0,0,s1,0, 0,0,0,s2, -s1,0,0,0, 0,-s2,0,0; def r = nc_algebra(C,D); setring(r); poly h = x1*(x1+1)*s1^2-2*x1*(x1+100)*s1+(x1+99)*(x1+100); facShift(h); } static proc sFacShift(poly h) " USAGE: A static procedure to factorize a polynomial in the first Shift algebra, where all the validity checks were made in advance. INPUT: A polynomial h in the first Shift Algebra. OUTPUT: A list of different factorizations of h, where the factors are irreducible ASSUMPTIONS: - The basering is the first Shift algebra and has n as first, and s as second variable, i.e. we have var(2)*var(1) = var(1)*var(2)+1 THEORY: If the given polynomial h is [0,1]-homogeneous, the routines for homogeneous factorizations are called. Otherwise we map the polynomial into the first Weyl algebra (the first shift algebra is a subring of the first Weyl algebra), and use facFirstWeyl to factorize it. Later we map the factors back, if possible. " {//proc sFacShift int p = printlevel - voice + 2; int i; int j ; string dbprintWhitespace = ""; number commonCoefficient = content(h); for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} //Checking if given polynomial is homogeneous if(homogwithorder(h,intvec(0,1))) {//The given polynomial is [0,1]-homogeneous dbprint(p,dbprintWhitespace+"The polynomial is [0,1]-homogeneous. Returning the homogeneous factorization"); return(homogfacFirstShift_all(h)); }//The given polynomial is [0,1]-homogeneous //---------- Start of interesting part ---------- dbprint(p,dbprintWhitespace+"Mapping the polynomial h into the first Weyl algebra."); poly temph = h/commonCoefficient; def ourBaseRing = basering; ring tempWeylAlgebraComm = 0,(x,d),dp; def tempWeylAlgebra = nc_algebra(1,1); setring(tempWeylAlgebra); map shiftMap = ourBaseRing, x*d, d; poly h = shiftMap(temph); dbprint(p,dbprintWhitespace+"Successful! The polynomial in the Weyl algebra is "+string(h)); dbprint(p,dbprintWhitespace+"Factorizing the polynomial in the first Weyl algebra"); list factorizationInWeyl = facFirstWeyl(h); dbprint(p,dbprintWhitespace+"Successful! The factorization is given by:"); dbprint(p,factorizationInWeyl); list validCombinations; dbprint(p,dbprintWhitespace+"Now we will map this back to the shift algebra and filter valid results"); //-Now we map the results back to the shift algebra. But first, we need to combine them properly. for (i = 1; i<=size(factorizationInWeyl); i++) {//Deleting the first Coefficient factor factorizationInWeyl[i] = delete(factorizationInWeyl[i],1); validCombinations = validCombinations + combineNonnegative(factorizationInWeyl[i]); }//Deleting the first Coefficient factor if (size(validCombinations) == 0) {//There are no valid combinations, therefore we can directly say, that h is irreducible setring(ourBaseRing); return(list(list(commonCoefficient, h/commonCoefficient))); }//There are no valid combinations, therefore we can directly say, that h is irreducible validCombinations = delete_dublicates_noteval(validCombinations); setring(ourBaseRing); map backFromWeyl = tempWeylAlgebra, var(1),var(2); list validCombinations = backFromWeyl(validCombinations); for (i = 1; i<=size(validCombinations); i++) { for (j = 1; j<=size(validCombinations[i]);j++) { setring(tempWeylAlgebra); fromWeylToShiftPoly(validCombinations[i][j],ourBaseRing); validCombinations[i][j] = result; kill result; kill tempResult; kill zeroPoly; kill fromWeyl; } } for (i = 1; i<=size(validCombinations); i++) {//Adding the common factor in the first position of the list validCombinations[i] = insert(validCombinations[i],commonCoefficient); }//Adding the common factor in the first position of the list dbprint(dbprintWhitespace+"Done."); //mapping return(validCombinations); }//proc sFacShift static proc sFacNthShift(poly h) " USAGE: A static procedure to factorize a polynomial in the nth Shift algebra, where all the validity checks were made in advance. INPUT: A polynomial h in the nth Shift Algebra. OUTPUT: A list of different factorizations of h, where the factors are irreducible ASSUMPTIONS: - The basering is the nth Shift algebra and the variables are given in the order (x_1, ... , x_n, s_1, ..., s_n), where s_i*x_i = (x_i+1)*s_i THEORY: We map the polynomial into the nth Weyl algebra (the nth shift algebra is a subring of the nth Weyl algebra), and use facWeyl to factorize it. Later we map the factors back, if possible. " {//proc sFacNthShift int p = printlevel - voice + 2; int i; int j ; string dbprintWhitespace = ""; number commonCoefficient = content(h); for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} //---------- Start of interesting part ---------- dbprint(p,dbprintWhitespace+"Mapping the polynomial h into the nth Weyl algebra."); poly temph = h/commonCoefficient; def ourBaseRing = basering; ring tempWeylAlgebraComm = 0,(x(1..(nvars(ourBaseRing) div 2)), d(1..(nvars(ourBaseRing) div 2))),dp; matrix C[nvars(ourBaseRing)][nvars(ourBaseRing)]; matrix D[nvars(ourBaseRing)][nvars(ourBaseRing)]; for (i = 1; i<= nvars(ourBaseRing); i++) {//Filling the matrices D and C for (j = 1; j<= nvars(ourBaseRing); j++) { C[i,j] = 1; if (j == (nvars(ourBaseRing) div 2) + i) { D[i,j] = 1; } else { D[i,j] = 0; } } }//Filling the matrices D and C def tempWeylAlgebra = nc_algebra(C,D); setring(tempWeylAlgebra); ideal shiftMapIdeal; for (i = 1; i<= nvars(tempWeylAlgebra) div 2; i++) { shiftMapIdeal[i] = x(i)*d(i); shiftMapIdeal[i + (nvars(tempWeylAlgebra) div 2)] = d(i); } map shiftMap = ourBaseRing, shiftMapIdeal; poly h = shiftMap(temph); dbprint(p,dbprintWhitespace+"Successful! The polynomial in the Weyl algebra is "+string(h)); dbprint(p,dbprintWhitespace+"Factorizing the polynomial in the nth Weyl algebra"); list factorizationInWeyl = facWeyl(h); dbprint(p,dbprintWhitespace+"Successful! The factorization is given by:"); dbprint(p,factorizationInWeyl); list validCombinations; dbprint(p,dbprintWhitespace+"Now we will map this back to the shift algebra and filter valid results"); //-Now we map the results back to the shift algebra. But first, we need to combine them properly. for (i = 1; i<=size(factorizationInWeyl); i++) {//Deleting the first Coefficient factor factorizationInWeyl[i] = delete(factorizationInWeyl[i],1); validCombinations = validCombinations + combineNonnegativeNthShift(factorizationInWeyl[i]); }//Deleting the first Coefficient factor if (size(validCombinations) == 0) {//There are no valid combinations, therefore we can directly say, that h is irreducible setring(ourBaseRing); return(list(list(commonCoefficient, h/commonCoefficient))); }//There are no valid combinations, therefore we can directly say, that h is irreducible validCombinations = delete_dublicates_noteval(validCombinations); setring(ourBaseRing); ideal backFromWeylIdeal; for (i = 1; i<=nvars(ourBaseRing); i++) { backFromWeylIdeal[i] = var(i); } map backFromWeyl = tempWeylAlgebra, backFromWeylIdeal; list validCombinations = backFromWeyl(validCombinations); for (i = 1; i<=size(validCombinations); i++) { for (j = 1; j<=size(validCombinations[i]);j++) { setring(tempWeylAlgebra); fromNthWeylToNthShiftPoly(validCombinations[i][j],ourBaseRing); validCombinations[i][j] = result; kill result; kill tempResult; kill zeroPoly; kill fromWeyl; } } for (i = 1; i<=size(validCombinations); i++) {//Adding the common factor in the first position of the list validCombinations[i] = insert(validCombinations[i],commonCoefficient); }//Adding the common factor in the first position of the list dbprint(dbprintWhitespace+"Done."); //mapping return(validCombinations); }//proc sFacNthShift //Tests: //(x1^2*s1 + s2)*(x2*s1+1); //(x1^2 + x2)*s1*s2; //(x2*x1 + s1)*(x2^2 + s2*s1) static proc combineNonnegative(list l) " USAGE: In sFacShift, when we want to map back the results of the factorization of the polynomial in the first Weyl algebra to the shift algebra. We need to recombine the factors such that we can map it back to the shift algebra without any problems. INPUT: A list l containing one factorization of a polynomial in the first Weyl algebra. For example for the polynomial (1+x)*(1+x+d) we would have the list [1,x+1,x+d+1]. OUTPUT:If we can map every factor without a problem back to the shift algebra (i.e. if the smallest homogeneous summand of every factor is of nonnegative degree), a list containing the same list as given in the input is returned. If otherwise some factors cause problems, we consider every possible combination (i.e. products of the factors) and extract those where all factors have a smallest homogeneous summand of nonnegative degree. ASSUMPTIONS: - Weyl algebra is given, and we have var(2)*var(1)=var(1)*var(2) +1 " {//combineNonnegative int p = printlevel - voice + 2; int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} //First the easy case: all of the factors fulfill the condition of mapping to shift: dbprint(p,dbprintWhitespace+"Checking, if the given factors can already be mapped without a problem."); int isValid = 1; for (i = 1; i<=size(l);i++) {//Checking for every entry if the condition is fulfilled. if (deg(l[i],intvec(1,-1))>0) {//Found one, where it is not fulfilled isValid = 0; break; }//Found one, where it is not fulfilled }//Checking for every entry if the condition is fulfilled. dbprint(p,dbprintWhitespace+"Done."); if (isValid) {//We can map every factor to the shift algebra and do not need to combine anything dbprint(p,dbprintWhitespace+"They can be mapped. Therefore we return them directly."); return(list(l)); }//We can map every factor to the shift algebra and do not need to combine anything dbprint(p,dbprintWhitespace+"They cannot be mapped. Looking for valid combinations."); //Starting with the case, where l only consists of 1 or two elements. if(size(l)<=2) {//The case where we won't call the function a second time if (deg(product(l),intvec(1,-1))>0) {//No way of a valid combination return(list()); }//No way of a valid combination else {//The product is the only possible and valid combination return(list(list(product(l)))); }//The product is the only possible and valid combination }//The case where we won't call the function a second time //---------- Easy pre-stuff done. now we combine the factors.---------- int pos; int j; int k; dbprint(p,dbprintWhitespace+"Making combinations of two."); list combinationsOfTwo = combinekfinlf(l,2); dbprint(p,dbprintWhitespace+"Done. Now checking, if there are valid ones in between."); list result; list validLHS; list validRHS; for (i = 1; i<=size(combinationsOfTwo); i++) {//go through all combinations and detect the valid ones if(deg(combinationsOfTwo[i][1],intvec(1,-1))>0 or deg(combinationsOfTwo[i][2],intvec(1,-1))>0) {//No chance, so no further treatment needed i++; continue; }//No chance, so no further treatment needed for (pos = 1; pos<=size(l);pos++) {//find the position where the combination splits if (product(l[1..pos]) == combinationsOfTwo[i][1]) {//Found the position break; }//Found the position }//find the position where the combination splits dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " + string(list(l[1..pos]))); validLHS = combineNonnegative(list(l[1..pos])); dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " + string(list(l[pos+1..size(l)]))); validRHS = combineNonnegative(list(l[pos+1..size(l)])); for (j = 1; j<=size(validLHS); j++) {//Combining the left hand side valid combnations... for (k = 1; k<=size(validRHS); k++) {//... with the right hand side valid combinations result = insert(result, validLHS[j]+validRHS[k]); }//... with the right hand side valid combinations }//Combining the left hand side valid combnations... }//go through all combinations and detect the valid ones result = delete_dublicates_noteval(result); dbprint(p,dbprintWhitespace+"Done."); return(result); }//combineNonnegative static proc combineNonnegativeNthShift(list l) " USAGE: In sFacNthShift, when we want to map back the results of the factorization of the polynomial in the nth Weyl algebra to the shift algebra. We need to recombine the factors such that we can map it back to the shift algebra without any problems. INPUT: A list l containing one factorization of a polynomial in the nth Weyl algebra. For example for the polynomial (1+x)*(1+x+d) we would have the list [1,x+1,x+d+1]. OUTPUT:If we can map every factor without a problem back to the shift algebra (i.e. if the smallest homogeneous summand of every factor is of nonnegative degree), a list containing the same list as given in the input is returned. If otherwise some factors cause problems, we consider every possible combination (i.e. products of the factors) and extract those where all factors have a smallest homogeneous summand of nonnegative degree. ASSUMPTIONS: - The nth Weyl algebra is given, and the variables are ordered in the form (x_1, ... , x_n, d_1, ..., d_n), with d_i*x_i = x_i * d_i + 1; " {//combineNonnegative int p = printlevel - voice + 2; int i; int j; intvec degreeOfInputPoly; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} //First the easy case: all of the factors fulfill the condition of mapping to shift: dbprint(p,dbprintWhitespace+"Checking, if the given factors can already be mapped without a problem."); int isValid = 1; for (i = 1; i<=size(l);i++) {//Checking for every entry if the condition is fulfilled. degreeOfInputPoly = degreeOfNthWeylPolyInverted(l[i]); for (j = 1; j<=size(degreeOfInputPoly); j++) {//checking if each entry has positive degree if (degreeOfInputPoly[j] > 0) { isValid = 0; break; } }//Found one, where it is not fulfilled }//Checking for every entry if the condition is fulfilled. dbprint(p,dbprintWhitespace+"Done."); if (isValid) {//We can map every factor to the shift algebra and do not need to combine anything dbprint(p,dbprintWhitespace+"They can be mapped. Therefore we return them directly."); return(list(l)); }//We can map every factor to the shift algebra and do not need to combine anything dbprint(p,dbprintWhitespace+"They cannot be mapped. Looking for valid combinations."); //Starting with the case, where l only consists of 1 or two elements. if(size(l)<=2) {//The case where we won't call the function a second time degreeOfInputPoly = degreeOfNthWeylPolyInverted(product(l)); for (j =1; j <= size(degreeOfInputPoly); j++) {//Checking if each entry has degree greater than zero or not if (degreeOfInputPoly[j] > 0) {//in this case, we can return false return(list()); }//in this case, we can return false }//Checking if each entry has degree greater than zero or not //If we encounter this line of code, we have two or less factors, //and they can be combined such that we can map it back to the nth //shift algebra return(list(list(product(l)))); }//The case where we won't call the function a second time //---------- Easy pre-stuff done. now we combine the factors.---------- int pos; int k; dbprint(p,dbprintWhitespace+"Making combinations of two."); list combinationsOfTwo = combinekfinlf(l,2); dbprint(p,dbprintWhitespace+"Done. Now checking, if there are valid ones in between."); list result; list validLHS; list validRHS; intvec degOfInpPoly1; intvec degOfInpPoly2; int noChance; for (i = 1; i<=size(combinationsOfTwo); i++) {//go through all combinations and detect the valid ones degOfInpPoly1 = degreeOfNthWeylPolyInverted(combinationsOfTwo[i][1]); degOfInpPoly2 = degreeOfNthWeylPolyInverted(combinationsOfTwo[i][2]); noChance = 0; for (j = 1; j<=size(degOfInpPoly1); j++) { if (degOfInpPoly1[j] > 0 or degOfInpPoly2[j] >0) {//No chance, so no further treatment needed noChance = 1; break; }//No chance, so no further treatment needed } if (noChance) { i++; continue; } for (pos = 1; pos<=size(l);pos++) {//find the position where the combination splits if (product(l[1..pos]) == combinationsOfTwo[i][1]) {//Found the position break; }//Found the position }//find the position where the combination splits dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " + string(list(l[1..pos]))); validLHS = combineNonnegativeNthShift(list(l[1..pos])); dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " + string(list(l[pos+1..size(l)]))); validRHS = combineNonnegativeNthShift(list(l[pos+1..size(l)])); for (j = 1; j<=size(validLHS); j++) {//Combining the left hand side valid combnations... for (k = 1; k<=size(validRHS); k++) {//... with the right hand side valid combinations result = insert(result, validLHS[j]+validRHS[k]); }//... with the right hand side valid combinations }//Combining the left hand side valid combnations... }//go through all combinations and detect the valid ones result = delete_dublicates_noteval(result); dbprint(p,dbprintWhitespace+"Done."); return(result); }//combineNonnegativeNthShift static proc fromWeylToShiftPoly(poly h, sAlgebra) " USAGE: Given a polynomial in the first Weyl algebra, this method returns it -- if possible -- as an element in the first shift algebra, which is given in the method header. INPUT: A polynomial h, and the first shift algebra as a ring OUTPUT: The correct mapping in the shift Algebra ASSUMPTIONS: - The lowest [-1,1]-homogeneous summand of h is of nonnegative degree - The shift algebra is given in the way that var(2)*var(1) = (var(1)+1)*var(2) " {//fromWeylToShiftPoly int p = printlevel - voice + 2; int i; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} if (deg(h,intvec(1,-1))>0) {//Wrong input polynomial ERROR("The lowest [-1,1] homogeneous summand of "+string(h)+" is of negative degree."); }//Wrong input polynomial def ourHomeBase = basering; list hDist = homogDistribution(h); setring(sAlgebra); poly result = 0; poly tempResult; poly zeroPoly; map fromWeyl = ourHomeBase, var(1), var(2); setring(ourHomeBase); poly zeroPoly; poly tempZeroPoly; int j; int k; int derDeg; for (i = 1; i<=size(hDist);i++) { derDeg = hDist[i][1]; setring(sAlgebra); tempResult = 1; setring(ourHomeBase); zeroPoly = lift(d^derDeg, hDist[i][2])[1,1]; for (j = 1; j<=size(zeroPoly); j++) { tempZeroPoly = zeroPoly[j]; setring(sAlgebra); zeroPoly = fromWeyl(tempZeroPoly); tempResult = tempResult * leadcoef(zeroPoly); setring(ourHomeBase); for (k = 1; k<=deg(zeroPoly[j],intvec(0,1));k++) { setring(sAlgebra); tempResult = tempResult*(var(1)-(k-1)); setring(ourHomeBase); } setring(sAlgebra); result = result + tempResult*var(2)^derDeg; tempResult = 1; setring(ourHomeBase); } } setring(sAlgebra); keepring(sAlgebra); }//fromWeylToShiftPoly static proc fromNthWeylToNthShiftPoly(poly h, sAlgebra) " USAGE: Given a polynomial in the nth Weyl algebra, this method returns it -- if possible -- as an element in the nth shift algebra, which is given in the method header. INPUT: A polynomial h, and the nth shift algebra as a ring OUTPUT: The correct mapping in the nth shift Algebra ASSUMPTIONS: - The lowest [-1,1]-homogeneous summand of h is of nonnegative degree - The shift algebra is given in the way, that the variables are sorted by (x_1, ... , x_n, s_1, ..., s_n) with s_i x_i = (x_i + 1) s_i " {//fromNthWeylToNthShiftPoly int p = printlevel - voice + 2; int i; string dbprintWhitespace = ""; intvec degOfh; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} degOfh = degreeOfNthWeylPolyInverted(h); for (i = 1; i<=size(degOfh); i++) { if (degOfh[i]>0) {//Wrong input polynomial ERROR("The lowest [-1,1] homogeneous summand of "+string(h)+" is of negative degree."); }//Wrong input polynomial } def ourHomeBase = basering; list hDist = homogDistributionNthWeyl(h); setring(sAlgebra); poly result = 0; poly tempResult; poly zeroPoly; ideal fromWeylMapideal; for (i = 1; i<=nvars(sAlgebra); i++) { fromWeylMapideal[i] = var(i); } map fromWeyl = ourHomeBase, fromWeylMapideal; setring(ourHomeBase); poly zeroPoly; poly tempZeroPoly; int j; int k; int l; intvec derDeg; intvec tempIntVec; for (i = 1; i<=size(hDist);i++) { derDeg = hDist[i][1]; setring(sAlgebra); tempResult = 1; setring(ourHomeBase); zeroPoly = hDist[i][2]; for (j = 1; j<=nvars(ourHomeBase) div 2; j++) {//lifting all the powers of the different ds away zeroPoly = lift(var((nvars(ourHomeBase) div 2) + j)^derDeg[j],zeroPoly)[1,1]; }//lifting all the powers of the different ds away for (j = 1; j<=size(zeroPoly); j++) { tempZeroPoly = zeroPoly[j]; setring(sAlgebra); zeroPoly = fromWeyl(tempZeroPoly); tempResult = tempResult * leadcoef(zeroPoly); setring(ourHomeBase); for (l = 1; l <= nvars(ourHomeBase) div 2; l++) {//iterating through all the d's tempIntVec = 0:nvars(ourHomeBase); tempIntVec[(nvars(ourHomeBase) div 2) + l] = 1; for (k = 1; k<=deg(zeroPoly[j],tempIntVec);k++) { setring(sAlgebra); tempResult = tempResult*(var(l)-(k-1)); setring(ourHomeBase); } }//iterating through all the d's setring(sAlgebra); for (l = 1; l<=nvars(sAlgebra) div 2; l++) { tempResult = tempResult*var((nvars(sAlgebra) div 2) + l)^derDeg[l]; } result = result + tempResult; tempResult = 1; setring(ourHomeBase); } } setring(sAlgebra); kill fromWeylMapideal; keepring(sAlgebra); }//fromNthWeylToNthShiftPoly static proc refineFactList(list L) { // assume: list L is an output of factorization proc // doing: remove doubled entries int s = size(L); int sm; int i,j,k,cnt; list M, U, A, B; A = L; k = 0; cnt = 1; for (i=1; i<=s; i++) { if (size(A[i]) != 0) { M = A[i]; // "probing with"; M; i; B[cnt] = M; cnt++; for (j=i+1; j<=s; j++) { if ( isEqualList(M,A[j]) ) { k++; // U consists of intvecs with equal pairs U[k] = intvec(i,j); A[j] = 0; } } } } kill A,U,M; return(B); } example { "EXAMPLE:";echo=2; ring R = 0,(x,s),dp; def r = nc_algebra(1,1); setring(r); list l,m; l = list(1,s2+1,x,s,x+s); m = l,list(1,s,x,s,x),l; refineFactList(m); } static proc isEqualList(list L, list M) { // int boolean: 1=yes, 0 =no : test whether two lists are identical int s = size(L); if (size(M)!=s) { return(0); } int j=1; while ( (L[j]==M[j]) && (j 0) {//We have two distict formulas for x and y. In this case use formula for x if (shift == 1) { result[i][j] = subst(result[i][j],theta,par(1)*theta + 1); } else { result[i][j] = subst(result[i][j], theta,par(1)^shift*theta+(par(1)^shift-1)/(par(1)-1)); } }//We have two distict formulas for x and y. In this case use formula for x } } }//adjust the a_0-parts }//Compute all possibilities to permute the x's resp. the y's in the list else {//The result is just all the permutations of the a_0-part result = permpp(list_azero); }//The result is just all the permutations of the a_0 part if (size(result)==0) { return(result); } dbprint(p," Done"); dbprint(p," Searching for theta resp. theta + 1 in the list and factorize them"); //Now we are going deeper and search for theta resp. theta + 1, substitute //them by xy resp. yx and go on permuting int found_theta; int thetapos; list leftpart; list rightpart; list lparts; list rparts; list tempadd; for (i = 1; i<=size(result) ; i++) {//checking every entry of result for theta or theta +1 found_theta = 0; for(j=1;j<=size(result[i]);j++) { if (result[i][j]==theta) {//the jth entry is theta and can be written as x*y thetapos = j; result[i]= insert(result[i],var(1),j-1); j++; result[i][j] = var(2); found_theta = 1; break; }//the jth entry is theta and can be written as x*y if(result[i][j] == par(1)*theta +1) { thetapos = j; result[i] = insert(result[i],var(2),j-1); j++; result[i][j] = var(1); found_theta = 1; break; } } if (found_theta) {//One entry was theta resp. theta +1 leftpart = result[i]; leftpart = leftpart[1..thetapos]; rightpart = result[i]; rightpart = rightpart[(thetapos+1)..size(rightpart)]; lparts = list(leftpart); rparts = list(rightpart); //first deal with the left part if (leftpart[thetapos] == var(1)) { shift_sign = 1; shiftvar = var(1); } else { shift_sign = -1; shiftvar = var(2); } for (j = size(leftpart); j>1;j--) {//drip x resp. y if (leftpart[j-1]==shiftvar) {//commutative j--; continue; }//commutative if (deg(leftpart[j-1],intvec(-1,1,0))!=0) {//stop here break; }//stop here //Here, we can only have a a0- part if (shift_sign<0) { leftpart[j] = subst(leftpart[j-1],theta, 1/par(1)*(theta +shift_sign)); } if (shift_sign>0) { leftpart[j] = subst(leftpart[j-1],theta, par(1)*theta + shift_sign); } leftpart[j-1] = shiftvar; lparts = lparts + list(leftpart); }//drip x resp. y //and now deal with the right part if (rightpart[1] == var(1)) { shift_sign = 1; shiftvar = var(1); } else { shift_sign = -1; shiftvar = var(2); } for (j = 1 ; j < size(rightpart); j++) { if (rightpart[j+1] == shiftvar) { j++; continue; } if (deg(rightpart[j+1],intvec(-1,1,0))!=0) { break; } if (shift_sign<0) { rightpart[j] = subst(rightpart[j+1], theta, par(1)*theta - shift_sign); } if (shift_sign>0) { rightpart[j] = subst(rightpart[j+1], theta, 1/par(1)*(theta - shift_sign)); } rightpart[j+1] = shiftvar; rparts = rparts + list(rightpart); } //And now, we put all possibilities together tempadd = list(); for (j = 1; j<=size(lparts); j++) { for (k = 1; k<=size(rparts);k++) { tempadd = tempadd + list(lparts[j]+rparts[k]); } } tempadd = delete(tempadd,1); // The first entry is already in the list result = result + tempadd; continue; //We can may be not be done already with the ith entry }//One entry was theta resp. theta +1 }//checking every entry of result for theta or theta +1 dbprint(p," Done"); //map back to the basering dbprint(p," Mapping back everything to the basering"); setring(r); map finalmap = tempRing, var(1), var(2),var(1)*var(2); list result = finalmap(result); for (i=1; i<=size(result);i++) {//adding the K factor result[i] = k_factor + result[i]; }//adding the k-factor dbprint(p," Done"); dbprint(p," Delete double entries in the list."); result = delete_dublicates_noteval(result); dbprint(p," Done"); return(result); }//proc HomogfacFirstQWeylAll_old proc homogfacFirstQWeyl_all(poly h) "USAGE: homogfacFirstQWeyl_all(h); h is a homogeneous polynomial in the first q-Weyl algebra with respect to the weight vector [-1,1] RETURN: list PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect to the weight vector [-1,1] in the first q-Weyl algebra THEORY: This function is a wrapper for homogFacNthQWeyl_all. It exists to make this library downward-compatible with older versions. SEE ALSO: homogfacFirstQWeyl "{//proc HomogfacFirstQWeylAll_old return(homogfacNthQWeyl_all(h)); }//proc HomogfacFirstQWeylAll_old example { "EXAMPLE:";echo=2; ring R = (0,q),(x,d),dp; def r = nc_algebra (q,1); setring(r); poly h = q^25*x^10*d^10+q^16*(q^4+q^3+q^2+q+1)^2*x^9*d^9+ q^9*(q^13+3*q^12+7*q^11+13*q^10+20*q^9+26*q^8+30*q^7+ 31*q^6+26*q^5+20*q^4+13*q^3+7*q^2+3*q+1)*x^8*d^8+ q^4*(q^9+2*q^8+4*q^7+6*q^6+7*q^5+8*q^4+6*q^3+ 4*q^2+2q+1)*(q^4+q^3+q^2+q+1)*(q^2+q+1)*x^7*d^7+ q*(q^2+q+1)*(q^5+2*q^4+2*q^3+3*q^2+2*q+1)*(q^4+q^3+q^2+q+1)*(q^2+1)*(q+1)*x^6*d^6+ (q^10+5*q^9+12*q^8+21*q^7+29*q^6+33*q^5+31*q^4+24*q^3+15*q^2+7*q+12)*x^5*d^5+ 6*x^3*d^3+24; homogfacFirstQWeyl_all(h); } //================================================== // Homogeneous factorization of the nth q-Weyl algebra //================================================== proc homogfacNthQWeyl(poly h) "USAGE: homogfacNthQWeyl(h); h is a homogeneous polynomial in the n'th q-Weyl algebra with respect to the weight vector [-1,...,-1,1,...,1]. \__ __/ \__ __/ \/ \/ n/2 n/2 RETURN: list PURPOSE: Computes a factorization of a homogeneous polynomial h in the n'th q-Weyl algebra THEORY:@code{homogfacNthQWeyl} returns a list with a factorization of the given, [-1,1]-homogeneous polynomial. For every i in 1..n: If the degree of the polynomial in [d_i,x_i] is k with k positive, the last entries in the output list are the second variable. If k is positive, the last k entries will be x_i. The other entries will be irreducible polynomials of degree zero or 1 resp. -1. resp. other variables GENERAL ASSUMPTIONS: - The basering is the nth Weyl algebra and has the form, that the first n variables represent x1, ..., xn, and the second n variables do represent the d1, ..., dn. - We have n parameters q_1,..., q_n given. SEE ALSO: homogfacFirstQWeyl, homogfacFirstQWeyl_all, homogfacNthQWeyl_all " {//proc homogfacNthQWeyl int p = printlevel-voice+2;//for dbprint poly hath = h; def r = basering; int i; int j; int k; string dbprintWhitespace = ""; for (i = 1; i<=voice;i++) {dbprintWhitespace = dbprintWhitespace + " ";} intvec ivm11 = intvec(-1,1); if(!checkIfProperNthQWeyl()) {//checking whether the given ring is proper q-Weyl ERROR("Assumptions on the ring structure not met."); return(list()); }//checking whether the given ring is proper q-Weyl if (!homogwithorderNthWeyl(h)) {//The given polynomial is not homogeneous ERROR("Given polynomial was not [-1,...,-1,1,...,1]-homogeneous"); return(list()); }//The given polynomial is not homogeneous if (h==0) { return(list(0)); } list result; intvec m = degreeOfNthWeylPoly(h); dbprint(p,dbprintWhitespace +" Splitting the polynomial in Q_0 and Q_k-Part"); dbprint(p,dbprintWhitespace + "Its [-1,...,-1,1,...,1] degree is "+string(m)); for (j = 1; j<=nvars(basering) div 2; j++) {//extracting the respective variable for every position dbprint(p,dbprintWhitespace + "Considering variables x_"+string(j)+" and d" + string(j)); if (m[j]!=0) {//The degree is not zero if (m[j] <0) {//There are more x than d hath = lift(var(j)^(-m[j]),hath)[1,1]; for (i = 1; i<=-m[j]; i++) { result = result + list(var(j)); } }//There are more x than d else {//There are more d than x hath = lift(var(nvars(basering) div 2 + j)^m[j],hath)[1,1]; for (i = 1; i<=m[j];i++) { result = result + list(var(nvars(basering) div 2 + j)); } }//There are more d than x }//The degree is not zero }//extracting the respective variable for every position dbprint(p,dbprintWhitespace+" Done"); //beginning to factor the zero-homogeneous part list mons; dbprint(p,dbprintWhitespace+" Putting the monomials in the Q_0-part in a list."); for(i = 1; i<=size(hath);i++) {//Putting the monomials in a list mons = mons+list(hath[i]); }//Putting the monomials in a list dbprint(p,dbprintWhitespace+" Done"); dbprint(p,dbprintWhitespace+" Mapping these monomials to K[theta_1,... , theta_n]"); list parameterNames; for (i = 1; i<=(nvars(basering) div 2); i++) {//saving the names of the parameters parameterNames[i] = ringlist(basering)[1][2][i]; }//saving the names of the parameters ring tempRingPre = (0,q(1..(nvars(basering) div 2))), (x(1..(nvars(basering) div 2)), d(1..(nvars(basering) div 2)), theta(1..(nvars(basering) div 2))),dp; list ringListTempRingPre = ringlist(tempRingPre); for (i = 1; i<=size(ringListTempRingPre[1][2]); i++) {//setting the same name for the parameters ringListTempRingPre[1][2][i] = parameterNames[i]; }//setting the same name for the parameters def tempRing = ring(ringListTempRingPre); setring(tempRing); ideal mapList; for (i = 1; i<=nvars(r) ; i++) {//filling the list of elements we want to map mapList[i] = var(i); }//filling the list of elements we want to map map thetamap = r,mapList; list mons = thetamap(mons); poly entry; intvec lExp; poly tempSummand; for (i = 1; i<=size(mons);i++) {//transforming the monomials as monomials in theta entry = leadcoef(mons[i]); lExp = leadexp(mons[i]); for (k = 1; k<=nvars(r) div 2; k++) {//iterating over the pairs x_kd_k for (j = 0; j0) { result[i][j] = subst(result[i][j],theta(k),par(k)^(shift[k])*theta(k) +(1-par(k)^(shift[k]))/(1-par(k))); } } } }//factor was a theta poly }//iterating through each factor }//adjust the a_0-parts }//Compute all possibilities to permute the x's resp. the y's in the list else {//The result is just all the permutations of the a_0-part result = permpp(list_azero); }//The result is just all the permutations of the a_0 part if (size(result)==0) { return(normalizeFactors(result)); } dbprint(p,dbprintWhitespace +" Done"); dbprint(p, dbprintWhitespace + "The factorization list is now:"); dbprint(p,result); dbprint(p, dbprintWhitespace + "Checking whether the intermediate result is correct or not"); dbprint(p,testNCfac(result)); dbprint(p,dbprintWhitespace +" Searching for theta resp. theta+1 in the list and fact. them"); //Now we are going deeper and search for theta resp. theta + 1, substitute //them by xy resp. yx and go on permuting int found_theta; int thetapos; int shift_sign; int thetaIndex; poly shiftvar; list leftpart; list rightpart; list lparts; list rparts; list tempadd; for (i = 1; i<=size(result) ; i++) {//checking every entry of result for theta or theta +1 found_theta = 0; for(j=1;j<=size(result[i]);j++) {//iterating through all factors for (k = 1; k<=nvars(r) div 2; k++) {//iterating through the variables if (result[i][j]==theta(k)) {//the jth entry is theta and can be written as x*y thetapos = j; thetaIndex = k; result[i]= insert(result[i],x(k),j-1); j++; result[i][j] = d(k); found_theta = 1; break; }//the jth entry is theta and can be written as x*y if(result[i][j] == theta(k) +1/par(k)) { thetapos = j; thetaIndex = k; result[i] = insert(result[i],d(k),j-1); j++; result[i][j] = x(k); found_theta = 1; break; } }//iterating through the variables if(found_theta) {break;} }//iterating through all factors if (found_theta) {//One entry was theta resp. theta +1 leftpart = result[i]; leftpart = leftpart[1..thetapos]; rightpart = result[i]; rightpart = rightpart[(thetapos+1)..size(rightpart)]; lparts = list(leftpart); rparts = list(rightpart); //first deal with the left part if (leftpart[thetapos] == x(thetaIndex)) { shift_sign = 1; shiftvar = x(thetaIndex); } else { shift_sign = -1; shiftvar = d(thetaIndex); } for (j = size(leftpart); j>1;j--) {//drip x resp. y if (leftpart[j-1]==shiftvar) {//commutative j--; continue; }//commutative if (leadexp(leftpart[j-1])[thetaIndex + nvars(r) div 2] - leadexp(leftpart[j-1])[thetaIndex]!=0) {//stop here break; }//stop here //Here, we can only have a a0- part if (shift_sign<0) { leftpart[j] = subst(leftpart[j-1],theta(thetaIndex), 1/par(thetaIndex)*(theta(thetaIndex)+shift_sign)); } else { if (shift_sign>0) { leftpart[j] = subst(leftpart[j-1],theta(thetaIndex), par(thetaIndex)*theta(thetaIndex)+shift_sign); } } leftpart[j-1] = shiftvar; lparts = lparts + list(leftpart); }//drip x resp. y //and now deal with the right part if (rightpart[1] == x(thetaIndex)) { shift_sign = 1; shiftvar = x(thetaIndex); } else { shift_sign = -1; shiftvar = d(thetaIndex); } for (j = 1 ; j < size(rightpart); j++) { if (rightpart[j+1] == shiftvar) { j++; continue; } if (leadexp(rightpart[j+1])[thetaIndex + nvars(r) div 2] - leadexp(rightpart[j+1])[thetaIndex]!=0) { break; } if (shift_sign<0) { rightpart[j] = subst(rightpart[j+1], theta(thetaIndex), par(thetaIndex)*theta(thetaIndex)+1); } else { if (shift_sign > 0) { rightpart[j] = subst(rightpart[j+1], theta(thetaIndex), 1/par(thetaIndex)*(theta(thetaIndex)-1)); } } rightpart[j+1] = shiftvar; rparts = rparts + list(rightpart); } //And now, we put all possibilities together tempadd = list(); for (j = 1; j<=size(lparts); j++) { for (k = 1; k<=size(rparts);k++) { tempadd = tempadd + list(lparts[j]+rparts[k]); } } tempadd = delete(tempadd,1); // The first entry is already in the list result = result + tempadd; continue; //We can may be not be done already with the ith entry }//One entry was theta resp. theta +1 }//checking every entry of result for theta or theta +1 dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace + "The new result list is:"); dbprint(result); setring(r); ideal finalMapList; for(i = 1; i<=nvars(r);i++) { finalMapList[i] = var(i); } for (i = 1; i<=nvars(r) div 2; i++) { finalMapList[i + nvars(r)] = var(i)*var(i + (nvars(r) div 2)); } map finalmap = tempRing,finalMapList; list result = finalmap(result); for (i=1; i<=size(result);i++) {//adding the K factor result[i] = k_factor + result[i]; }//adding the k-factor dbprint(p,dbprintWhitespace +" Done"); dbprint(p,dbprintWhitespace +" Delete double entries in the list."); result = delete_dublicates_noteval(result); dbprint(p,dbprintWhitespace +" Done"); return(normalizeFactors(result)); }//proc homogfacNthQWeyl_all example { "EXAMPLE:";echo=2; ring R = (0,q1,q2,q3),(x1,x2,x3,d1,d2,d3),dp; matrix C[6][6] = 1,1,1,q1,1,1, 1,1,1,1,q2,1, 1,1,1,1,1,q3, 1,1,1,1,1,1, 1,1,1,1,1,1, 1,1,1,1,1,1; matrix D[6][6] = 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1, -1,0,0,0,0,0, 0,-1,0,0,0,0, 0,0,-1,0,0,0; def r = nc_algebra(C,D); setring(r); poly h =x1*x2^2*x3^3*d1*d2^2+x2*x3^3*d2; homogfacNthQWeyl_all(h); } //================================================== // EASY EXAMPLES FOR WEYL ALGEBRA //================================================== /* Easy and fast example polynomials where one can find factorizations: K (x^2+d)*(x^2+d); (x^2+x)*(x^2+d); (x^3+x+1)*(x^4+d*x+2); (x^2*d+d)*(d+x*d); d^3+x*d^3+2*d^2+2*(x+1)*d^2+d+(x+2)*d; //Example 5 Grigoriev-Schwarz. (d+1)*(d+1)*(d+x*d); //Landau Example projected to the first dimension. */ //================================================== //Some Bugs(fixed)/hard examples from Martin Lee: //================================================== // ex1, ex2 /* ring s = 0,(x,d),Ws(-1,1); def S = nc_algebra(1,1); setring S; poly a = 10x5d4+26x4d5+47x5d2-97x4d3; //Not so hard any more... Done in around 4 minutes def l= facFirstWeyl (a); l; kill l; poly b = -5328x8d5-5328x7d6+720x9d2+720x8d3-16976x7d4-38880x6d5 -5184x7d3-5184x6d4-3774x5d5+2080x8d+5760x7d2-6144x6d3-59616x5d4 +3108x3d6-4098x6d2-25704x5d3-21186x4d4+8640x6d-17916x4d3+22680x2d5 +2040x5d-4848x4d2-9792x3d3+3024x2d4-10704x3d2-3519x2d3+34776xd4 +12096xd3+2898d4-5040x2d+8064d3+6048d2; //Still very hard... But it seems to be only because of the //combinatorial explosion def l= facFirstWeyl (b); l; // ex3: there was difference in answers => fixed LIB "ncfactor.lib"; ring r = 0,(x,y,z),dp; matrix D[3][3]; D[1,3]=-1; def R = nc_algebra(1,D); setring R; poly g= 7*z4*x+62*z3+26*z; def l1= facSubWeyl (g, x, z); l1; //---- other ring ring s = 0,(x,z),dp; def S = nc_algebra(1,-1); setring S; poly g= 7*z4*x+62*z3+26*z; def l2= facFirstWeyl (g); l2; map F = R,x,0,z; list l1 = F(l1); l1; //---- so the answers look different, check them! testNCfac(l2); // ok testNCfac(l1); // was not ok, but now it's been fixed!!! // selbst D und X so vertauschen dass sie erfuellt ist : ist gemacht */ /* // bug from M Lee LIB "ncfactor.lib"; ring s = 0,(z,x),dp; def S = nc_algebra(1,1); setring S; poly f= -60z4x2-54z4-56zx3-59z2x-64; def l= facFirstWeyl (f); l; // before: empty list; after fix: 1 entry, f is irreducible poly g = 75z3x2+92z3+24; def l= facFirstWeyl (g); l; //before: empty list, now: correct */ /* more things from Martin Lee; fixed ring R = 0,(x,s),dp; def r = nc_algebra(1,s); setring(r); poly h = (s2*x+x)*s; h= h* (x+s); def l= facFirstShift(h); l; // contained doubled entries: not anymore, fixed! ring R = 0,(x,s),dp; def r = nc_algebra(1,-1); setring(r); poly h = (s2*x+x)*s; h= h* (x+s); def l= facFirstWeyl(h); l; // contained doubled entries: not anymore, fixed! */ //====================================================================== //Examples from TestSuite that are terminating in a reasonable time. //====================================================================== //Counter example for old Algorithm, but now working: /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (1+x^2*d)^4; list lsng = facFirstWeyl(h); print(lsng); */ //Example 2.7. from Master thesis /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (xdd + xd+1+ (xd+5)*x)*(((x*d)^2+1)*d + xd+3+ (xd+7)*x); list lsng = facFirstWeyl(h); print(lsng); */ //Example with high combinatorial income /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (xdddd + (xd+1)*d*d+ (xd+5)*x*d*d)*(((x*d)^2+1)*d*x*x + (xd+3)*x*x+ (xd+7)*x*x*x); list lsng = facFirstWeyl(h); print(lsng); */ //Once a bug, now working /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (x^2*d^2+x)*(x+1); list lsng = facFirstWeyl(h); print(lsng); */ //Another one of that kind /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (x*d*d + (x*d)^5 +x)*((x*d+1)*d-(x*d-1)^5+x); list lsng = facFirstWeyl(h); print(lsng); */ //Example of Victor for Shift Algebra /* ring s = 0,(n,Sn),dp; def S = nc_algebra(1,Sn); setring S; LIB "ncfactor.lib"; list lsng = facFirstShift(n^2*Sn^2+3*n*Sn^2-n^2+2*Sn^2-3*n-2); print(lsng); */ //Interesting example, as there are actually also some complex solutions to it: /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "/Users/albertheinle/Studium/forschung/ncfactor/versionen/ncfactor.lib"; poly h =(x^3+x+1)*(x^4+d*x+2);//Example for finitely many, but more than one solution in between. list lsng = facFirstWeyl(h); print(lsng); */ //Another one of that kind: /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h =(x^2+d)*(x^2+d);//Example for finitely many, but more than one solution in between. list lsng = facFirstWeyl(h); print(lsng); */ //Example by W. Koepf: /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (x^4-1)*x*d^2+(1+7*x^4)*d+8*x^3; list lsng = facFirstWeyl(h); print(lsng); */ //Shift Example from W. Koepf /* ring R = 0,(n,s),dp; def r = nc_algebra(1,s); setring(r); LIB "ncfactor.lib"; poly h = n*(n+1)*s^2-2*n*(n+100)*s+(n+99)*(n+100); list lsng = facFirstShift(h); print(lsng); */ //Tsai Example... Once hard, now easy... /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (x^6+2*x^4-3*x^2)*d^2-(4*x^5-4*x^4-12*x^2-12*x)*d + (6*x^4-12*x^3-6*x^2-24*x-12); list lsng =facFirstWeyl(h); print(lsng); */ //====================================================================== // Hard examples not yet calculatable in feasible amount of time //====================================================================== //Also a counterexample for REDUCE. Very long Groebner basis computation in between. /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); LIB "ncfactor.lib"; poly h = (d^4+x^2+dx+x)*(d^2+x^4+xd+d); list lsng = facFirstWeyl(h); print(lsng); */ //Example from the Mainz-Group /* ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); poly dop6 = 1/35*x^4*(27-70*x+35*x^2)+ 1/35*x*(32+152*x+100*x^2-59*x^3+210*x^4+105*x^5)*d+ (-10368/35-67056/35*x-35512/7*x^2-50328/7*x^3-40240/7*x^4-2400*x^5-400*x^6)*d^2+ (-144/35*(x+1)*(1225*x^5+11025*x^4+37485*x^3+61335*x^2+50138*x+16584)-6912/35*(x+2)* (x+1)*(105*x^4+1155*x^3+4456*x^2+7150*x+4212) -27648/35*(x+3)*(x+1)*(35*x^2+350*x+867)* (x+2)^2)*d^3; LIB "ncfactor.lib"; printlevel = 5; facFirstWeyl(dop6); $;*/ //Another Mainz Example: /* LIB "ncfactor.lib"; ring R = 0,(x,d),dp; def r = nc_algebra(1,1); setring(r); poly dopp = 82547*x^4*d^4+60237*x^3*d^3+26772*x^5*d^5+2231*x^6*d^6+x*(1140138* x^2*d^2-55872*x*d-3959658*x^3*d^3-8381805*x^4*d^4-3089576*x^5*d^5-274786* x^6*d^6)+x^2*(-16658622*x*d-83427714*x^2*d^2-19715033*x^3*d^3+78915395*x^4 *d^4+35337930*x^5*d^5+3354194*x^6*d^6)+x^3*(-99752472-1164881352*x*d+ 4408536996*x^2*d^2+11774185985*x^3*d^3+5262196786*x^4*d^4+1046030561/2*x^5* d^5-10564451/2*x^6*d^6)+x^4*(-1925782272+21995375398*x*d+123415803356*x^2* d^2+302465300831/2*x^3*d^3+34140803907/2*x^4*d^4-15535653409*x^5*d^5-\ 2277687768*x^6*d^6)+x^5*(71273525520+691398212366*x*d+901772633569*x^2*d^2+ 2281275427069*x^3*d^3+2944352819911/2*x^4*d^4+836872370039/4*x^5*d^5+ 9066399237/4*x^6*d^6)+x^6*(2365174430376+9596715855542*x*d+29459572469704*x^ 2*d^2+92502197003786*x^3*d^3+65712473180525*x^4*d^4+13829360193674*x^5*d^5 +3231449477251/4*x^6*d^6)+x^7*(26771079436836+117709870166226*x*d+ 821686455179082*x^2*d^2+1803972139232179*x^3*d^3+1083654460691481*x^4*d^4+ 858903621851785/4*x^5*d^5+50096565802957/4*x^6*d^6)+x^8*(179341727601960+ 2144653944040630*x*d+13123246960284302*x^2*d^2+41138357917778169/2*x^3*d^3+ 20605819587976401/2*x^4*d^4+3677396642905423/2*x^5*d^5+402688260229369/4*x^6 *d^6)+x^9*(2579190935961288+43587063726809764*x*d+157045086382352387*x^2*d^ 2+172175668477370223*x^3*d^3+138636285385875407/2*x^4*d^4+10707836398626232* x^5*d^5+529435530567584*x^6*d^6)+x^10*(41501953525903392+558336731465626084* x*d+1407267553543222268*x^2*d^2+1153046693323226808*x^3*d^3+ 372331468563656085*x^4*d^4+48654019090240214*x^5*d^5+2114661191282167*x^6*d ^6)+x^11*(364526077273381884+4158060401095928464*x*d+8646807662899324262*x^2* d^2+5914675753405705400*x^3*d^3+1631934058875116005*x^4*d^4+ 187371894330537204*x^5*d^5+7366806367019734*x^6*d^6)+x^12*( 1759850321214603648+18265471270535733520*x*d+34201910114871110912*x^2*d^2+ 21265221434709398152*x^3*d^3+5437363546219595036*x^4*d^4+594029113431041060* x^5*d^5+22881659624561644*x^6*d^6)+x^13*(4648382639403200688+ 45699084277107816096*x*d+81049061578449009384*x^2*d^2+48858488665016574368*x ^3*d^3+12515362110098721444*x^4*d^4+1412152747420021048*x^5*d^5+ 57196947123984972*x^6*d^6)+x^14*(5459369397960020544+55837825300341621824*x* d+105671876924055409696*x^2*d^2+71551727420848766624*x^3*d^3+ 21094786205096577808*x^4*d^4+2695663190297032192*x^5*d^5+118791751565613264* x^6*d^6)+x^15*(1023333653580043776+47171127937488813824*x*d+ 157258351906685700352*x^2*d^2+145765192195300531840*x^3*d^3+ 49876215785510342176*x^4*d^4+6647374188802036864*x^5*d^5+287310278455067312* x^6*d^6)+x^16*(11960091747366236160+250326608568269289472*x*d+ 677587171115580981248*x^2*d^2+538246374825683603456*x^3*d^3+ 161380433451548754048*x^4*d^4+19149099315354950144*x^5*d^5+ 746433247985092544*x^6*d^6)+x^17*(42246252365448668160+657220532737851248640* x*d+1531751689216283911680*x^2*d^2+1090829514212206064640*x^3*d^3+ 299280728709430851840*x^4*d^4+32932767387222323200*x^5*d^5+ 1202281367574179840*x^6*d^6)+x^18*(6239106101942784000+320638742839606579200* x*d+873857213570556364800*x^2*d^2+645649080101933721600*x^3*d^3+ 177008238160627276800*x^4*d^4+19165088507111475200*x^5*d^5+ 683600826675660800*x^6*d^6)+x^19*(-60440251454613504000-476055211197689856000 *x*d-733497382597635072000*x^2*d^2-386038662982742016000*x^3*d^3-\ 83361486778142976000*x^4*d^4-7524999543181824000*x^5*d^5-232189492987008000* x^6*d^6)+x^20*(1578562930483200000+12628503443865600000*x*d+ 19732036631040000000*x^2*d^2+10523752869888000000*x^3*d^3+ 2302070940288000000*x^4*d^4+210475057397760000*x^5*d^5+6577345543680000*x^6* d^6); printlevel = 3; facFirstWeyl(dopp); */ //Hard Example by Viktor: /* ring r = 0,(x,d), (dp); def R = nc_algebra(1,1); setring R; LIB "ncfactor.lib"; poly t = x; poly D =d; poly p = 2*t^2*D^8-6*t*D^8+2*t^2*D^7+8*t*D^7+12*D^7-2*t^4*D^6+6*t^3*D^6+12*t*D^6-20*D^6 -2*t^4*D^5-8*t^3*D^5-4*t^2*D^5+12*t*D^5-28*D^5-12*t^3*D^4-4*t^2*D^4-4*t*D^4-24*D^4+4*t^4*D^3 -12*t^3*D^3+2*t^2*D^3-18*t*D^3+16*D^3+6*t^4*D^2-2*t^3*D^2+2*t^2*D^2-2*t*D^2+44*D^2+2*t^4*D +12*t^3*D+2*t*D+4*t^3-8; list lsng = facFirstWeyl(p); print(lsng); */ /*later hard example From Beals-Kartashova paper ring R = 0, (x1,x2,d1,d2),dp;def r = Weyl();setring(r); poly h = d1^2 -d2^2 + x1*d2 + x2*d1 + 1/4*(x2^2 - x1^2) +1; */