source: git/Singular/LIB/ncfactor.lib @ 32f8ae

spielwiese
Last change on this file since 32f8ae was f75297, checked in by Viktor Levandovskyy <levandov@…>, 13 years ago
*levandov: checkin from Aachen, bugfixes and improvements git-svn-id: file:///usr/local/Singular/svn/trunk@13924 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 82.2 KB
Line 
1///////////////////////////////////////////////////////////
2version = "$Id$";
3category="Noncommutative";
4info="
5LIBRARY: ncfactor.lib  Tools for factorization in some noncommutative algebras
6AUTHORS: Albert Heinle,     albert.heinle@rwth-aachen.de
7@*       Viktor Levandovskyy,     levandov@math.rwth-aachen.de
8
9OVERVIEW: In this library, new methods for factorization on polynomials
10  are implemented for two algebras, both generated by two generators (Weyl and
11  shift algebras) over a field K. Recall, that   the first Weyl algebra over K
12  is generated by x,d obeying the relation d*x=x*d+1.
13@* The first shift algebra over K is generated by x,s obeying the relation s*x=x*s+s.
14@* More detailled description of the algorithms can be found at
15@url{http://www.math.rwth-aachen.de/\~Albert.Heinle}.
16
17Guide: We are interested in computing a tree of factorizations, that is at the moment
18a list of all found factorizations is returned. It may contain factorizations, which
19are further reducible.
20
21PROCEDURES:
22  facFirstWeyl(h);    factorization in the first Weyl algebra
23  testNCfac(l[,h]);   tests factorizations from a given list for correctness
24  facSubWeyl(h,X,D);  factorization in the first Weyl algebra as a subalgebra
25  facFirstShift(h);   factorization in the first shift algebra
26";
27
28LIB "general.lib";
29LIB "nctools.lib";
30LIB "involut.lib";
31LIB "freegb.lib"; // for isVar
32
33proc tst_ncfactor()
34{
35  example facFirstWeyl;
36  example facFirstShift;
37  example facSubWeyl;
38  example testNCfac;
39}
40
41/////////////////////////////////////////////////////
42//==================================================*
43//deletes double-entries in a list of factorization
44//without evaluating the product.
45static proc delete_dublicates_noteval(list l)
46{//proc delete_dublicates_noteval
47  list result= l;
48  int j; int k; int i;
49  int deleted = 0;
50  int is_equal;
51  for (i = 1; i<= size(l); i++)
52  {//Iterate over the different factorizations
53    for (j = i+1; j<= size(l); j++)
54    {//Compare the i'th factorization to the j'th
55      if (size(l[i])!= size(l[j]))
56      {//different sizes => not equal
57        j++;
58        continue;
59      }//different sizes => not equal
60      is_equal = 1;
61      for (k = 1; k <= size(l[i]);k++)
62      {//Compare every entry
63        if (l[i][k]!=l[j][k])
64        {
65          is_equal = 0;
66          break;
67        }
68      }//Compare every entry
69      if (is_equal == 1)
70      {//Delete this entry, because there is another equal one int the list
71        result = delete(result, i-deleted);
72        deleted = deleted+1;
73        break;
74      }//Delete this entry, because there is another equal one int the list
75    }//Compare the i'th factorization to the j'th
76  }//Iterate over the different factorizations
77  return(result);
78}//proc delete_dublicates_noteval
79
80//==================================================
81//deletes the double-entries in a list with
82//evaluating the products
83static proc delete_dublicates_eval(list l)
84{//proc delete_dublicates_eval
85  list result=l;
86  int j; int k; int i;
87  int deleted = 0;
88  int is_equal;
89  for (i = 1; i<= size(result); i++)
90  {//Iterating over all elements in result
91    for (j = i+1; j<= size(result); j++)
92    {//comparing with the other elements
93      if (product(result[i]) == product(result[j]))
94      {//There are two equal results; throw away that one with the smaller size
95        if (size(result[i])>=size(result[j]))
96        {//result[i] has more entries
97          result = delete(result,j);
98          continue;
99        }//result[i] has more entries
100        else
101        {//result[j] has more entries
102          result = delete(result,i);
103          i--;
104          break;
105        }//result[j] has more entries
106      }//There are two equal results; throw away that one with the smaller size
107    }//comparing with the other elements
108  }//Iterating over all elements in result
109  return(result);
110}//proc delete_dublicates_eval
111
112
113//==================================================*
114//given a list of factors g and a desired size nof, the following
115//procedure combines the factors, such that we recieve a
116//list of the length nof.
117static proc combinekfinlf(list g, int nof, intvec limits) //nof stands for "number of factors"
118{//Procedure combinekfinlf
119  list result;
120  int i; int j; int k; //iteration variables
121  list fc; //fc stands for "factors combined"
122  list temp; //a temporary store for factors
123  def nofgl = size(g); //nofgl stands for "number of factors of the given list"
124  if (nofgl == 0)
125  {//g was the empty list
126    return(result);
127  }//g was the empty list
128  if (nof <= 0)
129  {//The user wants to recieve a negative number or no element as a result
130    return(result);
131  }//The user wants to recieve a negative number or no element as a result
132  if (nofgl == nof)
133  {//There are no factors to combine
134    if (limitcheck(g,limits))
135    {
136      result = result + list(g);
137    }
138    return(result);
139  }//There are no factors to combine
140  if (nof == 1)
141  {//User wants to get just one factor
142    if (limitcheck(list(product(g)),limits))
143    {
144      result = result + list(list(product(g)));
145    }
146    return(result);
147  }//User wants to get just one factor
148  for (i = nof; i > 1; i--)
149  {//computing the possibilities that have at least one original factor from g
150    for (j = i; j>=1; j--)
151    {//shifting the window of combinable factors to the left
152      //fc below stands for "factors combined"
153      fc = combinekfinlf(list(g[(j)..(j+nofgl - i)]),nof - i + 1,limits);
154      for (k = 1; k<=size(fc); k++)
155      {//iterating over the different solutions of the smaller problem
156        if (j>1)
157        {//There are g_i before the combination
158          if (j+nofgl -i < nofgl)
159          {//There are g_i after the combination
160            temp = list(g[1..(j-1)]) + fc[k] + list(g[(j+nofgl-i+1)..nofgl]);
161          }//There are g_i after the combination
162          else
163          {//There are no g_i after the combination
164            temp = list(g[1..(j-1)]) + fc[k];
165          }//There are no g_i after the combination
166        }//There are g_i before the combination
167        if (j==1)
168        {//There are no g_i before the combination
169          if (j+ nofgl -i <nofgl)
170          {//There are g_i after the combination
171            temp = fc[k]+ list(g[(j + nofgl - i +1)..nofgl]);
172          }//There are g_i after the combination
173        }//There are no g_i before the combination
174        if (limitcheck(temp,limits))
175        {
176          result = result + list(temp);
177        }
178      }//iterating over the different solutions of the smaller problem
179    }//shifting the window of combinable factors to the left
180  }//computing the possibilities that have at least one original factor from g
181  for (i = 2; i<=nofgl/nof;i++)
182  {//getting the other possible results
183    result = result + combinekfinlf(list(product(list(g[1..i])))+list(g[(i+1)..nofgl]),nof,limits);
184  }//getting the other possible results
185  result = delete_dublicates_noteval(result);
186  return(result);
187}//Procedure combinekfinlf
188
189
190//==================================================*
191//merges two sets of factors ignoring common
192//factors
193static proc merge_icf(list l1, list l2, intvec limits)
194{//proc merge_icf
195  list g;
196  list f;
197  int i; int j;
198  if (size(l1)==0)
199  {
200    return(list());
201  }
202  if (size(l2)==0)
203  {
204    return(list());
205  }
206  if (size(l2)<=size(l1))
207  {//l1 will be our g, l2 our f
208    g = l1;
209    f = l2;
210  }//l1 will be our g, l2 our f
211  else
212  {//l1 will be our f, l2 our g
213    g = l2;
214    f = l1;
215  }//l1 will be our f, l2 our g
216  def result = combinekfinlf(g,size(f),limits);
217  for (i = 1 ; i<= size(result); i++)
218  {//Adding the factors of f to every possibility listed in temp
219    for (j = 1; j<= size(f); j++)
220    {
221      result[i][j] = result[i][j]+f[j];
222    }
223    if(!limitcheck(result[i],limits))
224    {
225      result = delete(result,i);
226      i--;
227    }
228  }//Adding the factors of f to every possibility listed in temp
229  return(result);
230}//proc merge_icf
231
232//==================================================*
233//merges two sets of factors with respect to the occurrence
234//of common factors
235static proc merge_cf(list l1, list l2, intvec limits)
236{//proc merge_cf
237  list g;
238  list f;
239  int i; int j;
240  list pre;
241  list post;
242  list candidate;
243  list temp;
244  int temppos;
245  if (size(l1)==0)
246  {//the first list is empty
247    return(list());
248  }//the first list is empty
249  if(size(l2)==0)
250  {//the second list is empty
251    return(list());
252  }//the second list is empty
253  if (size(l2)<=size(l1))
254  {//l1 will be our g, l2 our f
255    g = l1;
256    f = l2;
257  }//l1 will be our g, l2 our f
258  else
259  {//l1 will be our f, l2 our g
260    g = l2;
261    f = l1;
262  }//l1 will be our f, l2 our g
263  list M;
264  for (i = 2; i<size(f); i++)
265  {//finding common factors of f and g...
266    for (j=2; j<size(g);j++)
267    {//... with g
268      if (f[i] == g[j])
269      {//we have an equal pair
270        M = M + list(list(i,j));
271      }//we have an equal pair
272    }//... with g
273  }//finding common factors of f and g...
274  if (g[1]==f[1])
275  {//Checking for the first elements to be equal
276    M = M + list(list(1,1));
277  }//Checking for the first elements to be equal
278  if (g[size(g)]==f[size(f)])
279  {//Checking for the last elements to be equal
280    M = M + list(list(size(f),size(g)));
281  }//Checking for the last elements to be equal
282  list result;//= list(list());
283  while(size(M)>0)
284  {//set of equal pairs is not empty
285    temp = M[1];
286    temppos = 1;
287    for (i = 2; i<=size(M); i++)
288    {//finding the minimal element of M
289      if (M[i][1]<=temp[1])
290      {//a possible candidate that is smaller than temp could have been found
291        if (M[i][1]==temp[1])
292        {//In this case we must look at the second number
293          if (M[i][2]< temp[2])
294          {//the candidate is smaller
295            temp = M[i];
296            temppos = i;
297          }//the candidate is smaller
298        }//In this case we must look at the second number
299        else
300        {//The candidate is definately smaller
301          temp = M[i];
302          temppos = i;
303        }//The candidate is definately smaller
304      }//a possible candidate that is smaller than temp could have been found
305    }//finding the minimal element of M
306    M = delete(M, temppos);
307    if(temp[1]>1)
308    {//There are factors to combine before the equal factor
309      if (temp[1]<size(f))
310      {//The most common case
311        //first the combinations ignoring common factors
312        pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
313        post = merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
314        for (i = 1; i <= size(pre); i++)
315        {//all possible pre's...
316          for (j = 1; j<= size(post); j++)
317          {//...combined with all possible post's
318            candidate = pre[i]+list(f[temp[1]])+post[j];
319            if (limitcheck(candidate,limits))
320            {
321              result = result + list(candidate);
322            }
323          }//...combined with all possible post's
324        }//all possible pre's...
325        //Now the combinations with respect to common factors
326        post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
327        if (size(post)>0)
328        {//There are factors to combine
329          for (i = 1; i <= size(pre); i++)
330          {//all possible pre's...
331            for (j = 1; j<= size(post); j++)
332            {//...combined with all possible post's
333              candidate= pre[i]+list(f[temp[1]])+post[j];
334              if (limitcheck(candidate,limits))
335              {
336                result = result + list(candidate);
337              }
338            }//...combined with all possible post's
339          }//all possible pre's...
340        }//There are factors to combine
341      }//The most common case
342      else
343      {//the last factor is the common one
344        pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
345        for (i = 1; i<= size(pre); i++)
346        {//iterating over the possible pre-factors
347          candidate = pre[i]+list(f[temp[1]]);
348          if (limitcheck(candidate,limits))
349          {
350            result = result + list(candidate);
351          }
352        }//iterating over the possible pre-factors
353      }//the last factor is the common one
354    }//There are factors to combine before the equal factor
355    else
356    {//There are no factors to combine before the equal factor
357      if (temp[1]<size(f))
358      {//Just a check for security
359        //first without common factors
360        post=merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
361        for (i = 1; i<=size(post); i++)
362        {
363          candidate = list(f[temp[1]])+post[i];
364          if (limitcheck(candidate,limits))
365          {
366            result = result + list(candidate);
367          }
368        }
369        //Now with common factors
370        post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
371        if(size(post)>0)
372        {//we could find other combinations
373          for (i = 1; i<=size(post); i++)
374          {
375            candidate = list(f[temp[1]])+post[i];
376            if (limitcheck(candidate,limits))
377            {
378              result = result + list(candidate);
379            }
380          }
381        }//we could find other combinations
382      }//Just a check for security
383    }//There are no factors to combine before the equal factor
384  }//set of equal pairs is not empty
385  return(result);
386}//proc merge_cf
387
388
389//==================================================*
390//merges two sets of factors
391
392static proc mergence(list l1, list l2, intvec limits)
393{//Procedure mergence
394  list g;
395  list f;
396  int l; int k;
397  list F;
398  if (size(l2)<=size(l1))
399  {//l1 will be our g, l2 our f
400    g = l1;
401    f = l2;
402  }//l1 will be our g, l2 our f
403  else
404  {//l1 will be our f, l2 our g
405    g = l2;
406    f = l1;
407  }//l1 will be our f, l2 our g
408  list result;
409  for (l = size(f); l>=1; l--)
410  {//all possibilities to combine the factors of f
411    F = combinekfinlf(f,l,limits);
412    for (k = 1; k<= size(F);k++)
413    {//for all possibilities of combinations of the factors of f
414      result = result + merge_cf(F[k],g,limits);
415      result = result + merge_icf(F[k],g,limits);
416    }//for all possibilities of combinations of the factors of f
417  }//all possibilities to combine the factors of f
418  return(result);
419}//Procedure mergence
420
421
422//==================================================
423//Checks, whether a list of factors doesn't exceed the given limits
424static proc limitcheck(list g, intvec limits)
425{//proc limitcheck
426  int i;
427  if (size(limits)!=3)
428  {//check the input
429    return(0);
430  }//check the input
431  if(size(g)==0)
432  {
433    return(0);
434  }
435  def prod = product(g);
436  intvec iv11 = intvec(1,1);
437  intvec iv10 = intvec(1,0);
438  intvec iv01 = intvec(0,1);
439  def limg = intvec(deg(prod,iv11) ,deg(prod,iv10),deg(prod,iv01));
440  for (i = 1; i<=size(limg);i++)
441  {//the final check
442    if(limg[i]>limits[i])
443    {
444      return(0);
445    }
446  }//the final check
447  return(1);
448}//proc limitcheck
449
450
451//==================================================*
452//one factorization of a homogeneous polynomial
453//in the first Weyl Algebra
454static proc homogfacFirstWeyl(poly h)
455"USAGE: homogfacFirstWeyl(h); h is a homogeneous polynomial in the
456        first Weyl algebra with respect to the weight vector [-1,1]
457RETURN: list
458PURPOSE: Computes a factorization of a homogeneous polynomial h with
459         respect to the weight vector [-1,1] in the first Weyl algebra
460THEORY: @code{homogfacFirstWeyl} returns a list with a factorization of the given,
461        [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with
462        k positive, the last k entries in the output list are the second
463        variable. If k is positive, the last k entries will be x. The other
464        entries will be irreducible polynomials of degree zero or 1 resp. -1.
465SEE ALSO: homogfacFirstWeyl_all
466"{//proc homogfacFirstWeyl
467  int p = printlevel-voice+2;//for dbprint
468  def r = basering;
469  poly hath;
470  int i; int j;
471  intvec ivm11 = intvec(-1,1);
472  if (!homogwithorder(h,ivm11))
473  {//The given polynomial is not homogeneous
474    ERROR("Given polynomial was not [-1,1]-homogeneous");
475    return(list());
476  }//The given polynomial is not homogeneous
477  if (h==0)
478  {
479    return(list(0));
480  }
481  list result;
482  int m = deg(h,ivm11);
483  dbprint(p,"==> Splitting the polynomial in A_0 and A_k-Part");
484  if (m!=0)
485  {//The degree is not zero
486    if (m <0)
487    {//There are more x than y
488      hath = lift(var(1)^(-m),h)[1,1];
489      for (i = 1; i<=-m; i++)
490      {
491        result = result + list(var(1));
492      }
493    }//There are more x than y
494    else
495    {//There are more y than x
496      hath = lift(var(2)^m,h)[1,1];
497      for (i = 1; i<=m;i++)
498      {
499        result = result + list(var(2));
500      }
501    }//There are more y than x
502  }//The degree is not zero
503  else
504  {//The degree is zero
505    hath = h;
506  }//The degree is zero
507  dbprint(p,"==> Done");
508  //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1)
509  list mons;
510  dbprint(p,"==> Putting the monomials in the A_0-part in a list.");
511  for(i = 1; i<=size(hath);i++)
512  {//Putting the monomials in a list
513    mons = mons+list(hath[i]);
514  }//Putting the monomials in a list
515  dbprint(p,"==> Done");
516  dbprint(p,"==> Mapping this monomials to K[theta]");
517  ring tempRing = 0,(x,y,theta),dp;
518  setring tempRing;
519  map thetamap = r,x,y;
520  list mons = thetamap(mons);
521  poly entry;
522  for (i = 1; i<=size(mons);i++)
523  {//transforming the monomials as monomials in theta
524    entry = leadcoef(mons[i]);
525    for (j = 0; j<leadexp(mons[i])[2];j++)
526    {
527      entry = entry * (theta-j);
528    }
529    mons[i] = entry;
530  }//transforming the monomials as monomials in theta
531  dbprint(p,"==> Done");
532  dbprint(p,"==> Factorize the A_0-Part in K[theta]");
533  list azeroresult = factorize(sum(mons));
534  dbprint(p,"==> Successful");
535  list azeroresult_return_form;
536  for (i = 1; i<=size(azeroresult[1]);i++)
537  {//rewrite the result of the commutative factorization
538    for (j = 1; j <= azeroresult[2][i];j++)
539    {
540      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
541    }
542  }//rewrite the result of the commutative factorization
543  dbprint(p,"==> Mapping back to A_0.");
544  setring(r);
545  map finalmap = tempRing,var(1),var(2),var(1)*var(2);
546  list tempresult = finalmap(azeroresult_return_form);
547  dbprint(p,"Successful.");
548  for (i = 1; i<=size(tempresult);i++)
549  {//factorizations of theta resp. theta +1
550    if(tempresult[i]==var(1)*var(2))
551    {
552      tempresult = insert(tempresult,var(1),i-1);
553      i++;
554      tempresult[i]=var(2);
555    }
556    if(tempresult[i]==var(2)*var(1))
557    {
558      tempresult = insert(tempresult,var(2),i-1);
559      i++;
560      tempresult[i]=var(1);
561    }
562  }//factorizations of theta resp. theta +1
563  result = tempresult+result;
564  return(result);
565}//proc homogfacFirstWeyl
566/* example */
567/* { */
568/*      "EXAMPLE:";echo=2; */
569/*      ring R = 0,(x,y),Ws(-1,1); */
570/*      def r = nc_algebra(1,1); */
571/*      setring(r); */
572/*      poly h = (x^2*y^2+1)*(x^4); */
573/*      homogfacFirstWeyl(h); */
574/* } */
575
576//==================================================
577//Computes all possible homogeneous factorizations
578static proc homogfacFirstWeyl_all(poly h)
579"USAGE: homogfacFirstWeyl_all(h); h is a homogeneous polynomial in the first Weyl algebra
580        with respect to the weight vector [-1,1]
581RETURN: list
582PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect
583         to the weight vector [-1,1] in the first Weyl algebra
584THEORY: @code{homogfacFirstWeyl} returns a list with all factorization of the given,
585        homogeneous polynomial. It uses the output of homogfacFirstWeyl and permutes
586        its entries with respect to the commutation rule. Furthermore, if a
587        factor of degree zero is irreducible in K[\theta], but reducible in
588        the first Weyl algebra, the permutations of this element with the other
589        entries will also be computed.
590SEE ALSO: homogfacFirstWeyl
591"{//proc HomogfacFirstWeylAll
592  int p=printlevel-voice+2;//for dbprint
593  intvec iv11= intvec(1,1);
594  if (deg(h,iv11) <= 0 )
595  {//h is a constant
596    dbprint(p,"Given polynomial was not homogeneous");
597    return(list(list(h)));
598  }//h is a constant
599  def r = basering;
600  list one_hom_fac; //stands for one homogeneous factorization
601  int i; int j; int k;
602  intvec ivm11 = intvec(-1,1);
603  dbprint(p,"==> Calculate one homogeneous factorization using homogfacFirstWeyl");
604  //Compute again a homogeneous factorization
605  one_hom_fac = homogfacFirstWeyl(h);
606  dbprint(p,"Successful");
607  if (size(one_hom_fac) == 0)
608  {//there is no homogeneous factorization or the polynomial was not homogeneous
609    return(list());
610  }//there is no homogeneous factorization or the polynomial was not homogeneous
611  //divide list in A0-Part and a list of x's resp. y's
612  list list_not_azero = list();
613  list list_azero;
614  list k_factor;
615  int is_list_not_azero_empty = 1;
616  int is_list_azero_empty = 1;
617  k_factor = list(one_hom_fac[1]);
618  if (absValue(deg(h,ivm11))<size(one_hom_fac)-1)
619  {//There is a nontrivial A_0-part
620    list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,ivm11)))];
621    is_list_azero_empty = 0;
622  }//There is a nontrivial A_0 part
623  dbprint(p,"==> Combine x,y to xy in the factorization again.");
624  for (i = 1; i<=size(list_azero)-1;i++)
625  {//in homogfacFirstWeyl, we factorized theta, and this will be made undone
626    if (list_azero[i] == var(1))
627    {
628      if (list_azero[i+1]==var(2))
629      {
630        list_azero[i] = var(1)*var(2);
631        list_azero = delete(list_azero,i+1);
632      }
633    }
634    if (list_azero[i] == var(2))
635    {
636      if (list_azero[i+1]==var(1))
637      {
638        list_azero[i] = var(2)*var(1);
639        list_azero = delete(list_azero,i+1);
640      }
641    }
642  }//in homogfacFirstWeyl, we factorized theta, and this will be made undone
643  dbprint(p,"==> Done");
644  if(deg(h,ivm11)!=0)
645  {//list_not_azero is not empty
646    list_not_azero =
647      one_hom_fac[(size(one_hom_fac)-absValue(deg(h,ivm11))+1)..size(one_hom_fac)];
648    is_list_not_azero_empty = 0;
649  }//list_not_azero is not empty
650  //Map list_azero in K[theta]
651  dbprint(p,"==> Map list_azero to K[theta]");
652  ring tempRing = 0,(x,y,theta), dp;
653  setring(tempRing);
654  poly entry;
655  map thetamap = r,x,y;
656  if(!is_list_not_azero_empty)
657  {//Mapping in Singular is only possible, if the list before
658    //contained at least one element of the other ring
659    list list_not_azero = thetamap(list_not_azero);
660  }//Mapping in Singular is only possible, if the list before
661  //contained at least one element of the other ring
662  if(!is_list_azero_empty)
663  {//Mapping in Singular is only possible, if the list before
664    //contained at least one element of the other ring
665    list list_azero= thetamap(list_azero);
666  }//Mapping in Singular is only possible, if the list before
667  //contained at least one element of the other ring
668  list k_factor = thetamap(k_factor);
669  list tempmons;
670  dbprint(p,"==> Done");
671  for(i = 1; i<=size(list_azero);i++)
672  {//rewrite the polynomials in A1 as polynomials in K[theta]
673    tempmons = list();
674    for (j = 1; j<=size(list_azero[i]);j++)
675    {
676      tempmons = tempmons + list(list_azero[i][j]);
677    }
678    for (j = 1 ; j<=size(tempmons);j++)
679    {
680      entry = leadcoef(tempmons[j]);
681      for (k = 0; k < leadexp(tempmons[j])[2];k++)
682      {
683        entry = entry*(theta-k);
684      }
685      tempmons[j] = entry;
686    }
687    list_azero[i] = sum(tempmons);
688  }//rewrite the polynomials in A1 as polynomials in K[theta]
689  //Compute all permutations of the A0-part
690  dbprint(p,"==> Compute all permutations of the A_0-part with the first resp. the snd. variable");
691  list result;
692  int shift_sign;
693  int shift;
694  poly shiftvar;
695  if (size(list_not_azero)!=0)
696  {//Compute all possibilities to permute the x's resp. the y's in the list
697    if (list_not_azero[1] == x)
698    {//h had a negative weighted degree
699      shift_sign = 1;
700      shiftvar = x;
701    }//h had a negative weighted degree
702    else
703    {//h had a positive weighted degree
704      shift_sign = -1;
705      shiftvar = y;
706    }//h had a positive weighted degree
707    result = permpp(list_azero + list_not_azero);
708    for (i = 1; i<= size(result); i++)
709    {//adjust the a_0-parts
710      shift = 0;
711      for (j=1; j<=size(result[i]);j++)
712      {
713        if (result[i][j]==shiftvar)
714        {
715          shift = shift + shift_sign;
716        }
717        else
718        {
719          result[i][j] = subst(result[i][j],theta,theta + shift);
720        }
721      }
722    }//adjust the a_0-parts
723  }//Compute all possibilities to permute the x's resp. the y's in the list
724  else
725  {//The result is just all the permutations of the a_0-part
726    result = permpp(list_azero);
727  }//The result is just all the permutations of the a_0 part
728  if (size(result)==0)
729  {
730    return(result);
731  }
732  dbprint(p,"==> Done");
733  dbprint(p,"==> Searching for theta resp. theta + 1 in the list and factorize them");
734  //Now we are going deeper and search for theta resp. theta + 1, substitute
735  //them by xy resp. yx and go on permuting
736  int found_theta;
737  int thetapos;
738  list leftpart;
739  list rightpart;
740  list lparts;
741  list rparts;
742  list tempadd;
743  for (i = 1; i<=size(result) ; i++)
744  {//checking every entry of result for theta or theta +1
745    found_theta = 0;
746    for(j=1;j<=size(result[i]);j++)
747    {
748      if (result[i][j]==theta)
749      {//the jth entry is theta and can be written as x*y
750        thetapos = j;
751        result[i]= insert(result[i],x,j-1);
752        j++;
753        result[i][j] = y;
754        found_theta = 1;
755        break;
756      }//the jth entry is theta and can be written as x*y
757      if(result[i][j] == theta +1)
758      {
759        thetapos = j;
760        result[i] = insert(result[i],y,j-1);
761        j++;
762        result[i][j] = x;
763        found_theta = 1;
764        break;
765      }
766    }
767    if (found_theta)
768    {//One entry was theta resp. theta +1
769      leftpart = result[i];
770      leftpart = leftpart[1..thetapos];
771      rightpart = result[i];
772      rightpart = rightpart[(thetapos+1)..size(rightpart)];
773      lparts = list(leftpart);
774      rparts = list(rightpart);
775      //first deal with the left part
776      if (leftpart[thetapos] == x)
777      {
778        shift_sign = 1;
779        shiftvar = x;
780      }
781      else
782      {
783        shift_sign = -1;
784        shiftvar = y;
785      }
786      for (j = size(leftpart); j>1;j--)
787      {//drip x resp. y
788        if (leftpart[j-1]==shiftvar)
789        {//commutative
790          j--;
791          continue;
792        }//commutative
793        if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
794        {//stop here
795          break;
796        }//stop here
797        //Here, we can only have a a0- part
798        leftpart[j] = subst(leftpart[j-1],theta, theta + shift_sign);
799        leftpart[j-1] = shiftvar;
800        lparts = lparts + list(leftpart);
801      }//drip x resp. y
802      //and now deal with the right part
803      if (rightpart[1] == x)
804      {
805        shift_sign = 1;
806        shiftvar = x;
807      }
808      else
809      {
810        shift_sign = -1;
811        shiftvar = y;
812      }
813      for (j = 1 ; j < size(rightpart); j++)
814      {
815        if (rightpart[j+1] == shiftvar)
816        {
817          j++;
818          continue;
819        }
820        if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
821        {
822          break;
823        }
824        rightpart[j] = subst(rightpart[j+1], theta, theta - shift_sign);
825        rightpart[j+1] = shiftvar;
826        rparts = rparts + list(rightpart);
827      }
828      //And now, we put all possibilities together
829      tempadd = list();
830      for (j = 1; j<=size(lparts); j++)
831      {
832        for (k = 1; k<=size(rparts);k++)
833        {
834          tempadd = tempadd + list(lparts[j]+rparts[k]);
835        }
836      }
837      tempadd = delete(tempadd,1); // The first entry is already in the list
838      result = result + tempadd;
839      continue; //We can may be not be done already with the ith entry
840    }//One entry was theta resp. theta +1
841  }//checking every entry of result for theta or theta +1
842  dbprint(p,"==> Done");
843  //map back to the basering
844  dbprint(p,"==> Mapping back everything to the basering");
845  setring(r);
846  map finalmap = tempRing, var(1), var(2),var(1)*var(2);
847  list result = finalmap(result);
848  for (i=1; i<=size(result);i++)
849  {//adding the K factor
850    result[i] = k_factor + result[i];
851  }//adding the k-factor
852  dbprint(p,"==> Done");
853  dbprint(p,"==> Delete double entries in the list.");
854  result = delete_dublicates_noteval(result);
855  dbprint(p,"==> Done");
856  return(result);
857}//proc HomogfacFirstWeylAll
858/* example */
859/* { */
860/*      "EXAMPLE:";echo=2; */
861/*      ring R = 0,(x,y),Ws(-1,1); */
862/*      def r = nc_algebra(1,1); */
863/*      setring(r); */
864/*      poly h = (x^2*y^2+1)*(x^4); */
865/*      homogfacFirstWeyl_all(h); */
866/* } */
867
868//==================================================*
869//Computes all permutations of a given list
870static proc perm(list l)
871{//proc perm
872  int i; int j;
873  list tempresult;
874  list result;
875  if (size(l)==0)
876  {
877    return(list());
878  }
879  if (size(l)==1)
880  {
881    return(list(l));
882  }
883  for (i = 1; i<=size(l); i++ )
884  {
885    tempresult = perm(delete(l,i));
886    for (j = 1; j<=size(tempresult);j++)
887    {
888      tempresult[j] = list(l[i])+tempresult[j];
889    }
890    result = result+tempresult;
891  }
892  return(result);
893}//proc perm
894
895//==================================================
896//computes all permutations of a given list by
897//ignoring equal entries (faster than perm)
898static proc permpp(list l)
899{//proc permpp
900  int i; int j;
901  list tempresult;
902  list l_without_double;
903  list l_without_double_pos;
904  int double_entry;
905  list result;
906  if (size(l)==0)
907  {
908    return(list());
909  }
910  if (size(l)==1)
911  {
912    return(list(l));
913  }
914  for (i = 1; i<=size(l);i++)
915  {//Filling the list with unique entries
916    double_entry = 0;
917    for (j = 1; j<=size(l_without_double);j++)
918    {
919      if (l_without_double[j] == l[i])
920      {
921        double_entry = 1;
922        break;
923      }
924    }
925    if (!double_entry)
926    {
927      l_without_double = l_without_double + list(l[i]);
928      l_without_double_pos = l_without_double_pos + list(i);
929    }
930  }//Filling the list with unique entries
931  for (i = 1; i<=size(l_without_double); i++ )
932  {
933    tempresult = permpp(delete(l,l_without_double_pos[i]));
934    for (j = 1; j<=size(tempresult);j++)
935    {
936      tempresult[j] = list(l_without_double[i])+tempresult[j];
937    }
938    result = result+tempresult;
939  }
940  return(result);
941}//proc permpp
942
943//==================================================
944//factorization of the first Weyl Algebra
945
946//The following procedure just serves the purpose to
947//transform the input into an appropriate input for
948//the procedure sfacwa, where the ring must contain the
949//variables in a certain order.
950proc facFirstWeyl(poly h)
951"USAGE: facFirstWeyl(h); h a polynomial in the first Weyl algebra
952RETURN: list
953PURPOSE: compute all factorizations of a polynomial in the first Weyl algebra
954THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle
955ASSUME: basering in the first Weyl algebra
956NOTE: Every entry of the output list is a list with factors for one possible factorization.
957The first factor is always a constant (1, if no nontrivial constant could be excluded).
958EXAMPLE: example facFirstWeyl; shows examples
959SEE ALSO: facSubWeyl, testNCfac, facFirstShift
960"{//proc facFirstWeyl
961  //Definition of printlevel variable
962  int p = printlevel-voice+2;
963  dbprint(p,"==> Checking if the given algebra is a Weyl algebra");
964  //Redefine the ring in my standard form
965  if (!isWeyl())
966  {//Our basering is not the Weyl algebra
967    ERROR("Ring was not the first Weyl algebra");
968    return(list());
969  }//Our basering is not the Weyl algebra
970  dbprint(p,"==> Successful");
971  dbprint(p,"==> Checking, if the given ring is the first Weyl algebra");
972  if(nvars(basering)!=2)
973  {//Our basering is the Weyl algebra, but not the first
974    ERROR("==>Ring is not the first Weyl algebra");
975    return(list());
976  }//Our basering is the Weyl algebra, but not the first
977  dbprint(p,"==> Successful");
978  list result = list();
979  int i;int j; int k; int l; //counter
980  if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why
981  {
982    dbprint(p,"==> positions of the variables have to be switched");
983    def r = basering;
984    ring tempRing = ringlist(r)[1][1],(x,y),Ws(-1,1); // very strange:
985    // setting Wp(-1,1) leads to SegFault; to clarify why!!!
986    def NTR = nc_algebra(1,1);
987    setring NTR ;
988    map transf = r, var(2), var(1);
989    dbprint(p,"==> Successful");
990    dbprint(p, "==> factorization of the polynomial with the routine sfacwa");
991    list resulttemp = sfacwa(transf(h));
992    dbprint(p,"==> Done");
993    setring(r);
994    map transfback = NTR, var(2),var(1);
995    result = transfback(resulttemp);
996  }
997  else
998  {
999    dbprint(p, "==> factorization of the polynomial with the routine sfacwa");
1000    result = sfacwa(h);
1001    dbprint(p,"==> Done");
1002  }
1003  dbprint(p,"==> recursively check factors for irreducibility");
1004  list recursivetemp;
1005  for(i = 1; i<=size(result);i++)
1006  {//recursively factorize factors
1007    if(size(result[i])>2)
1008    {//Nontrivial factorization
1009      for (j=2;j<=size(result[i]);j++)
1010      {//Factorize every factor
1011        recursivetemp = facFirstWeyl(result[i][j]);
1012        if(size(recursivetemp)>1)
1013        {//we have a nontrivial factorization
1014          for(k=1; k<=size(recursivetemp);k++)
1015          {//insert factorized factors
1016            if(size(recursivetemp[k])>2)
1017            {//nontrivial
1018              result = insert(result,result[i],i);
1019              for(l = size(recursivetemp[k]);l>=2;l--)
1020              {
1021                result[i+1] = insert(result[i+1],recursivetemp[k][l],j);
1022              }
1023              result[i+1] = delete(result[i+1],j);
1024            }//nontrivial
1025          }//insert factorized factors
1026        }//we have a nontrivial factorization
1027      }//Factorize every factor
1028    }//Nontrivial factorization
1029  }//recursively factorize factors
1030  dbprint(p,"==> Done");
1031  if (size(result)==0)
1032  {//only the trivial factorization could be found
1033    result = list(list(1,h));
1034  }//only the trivial factorization could be found
1035  //now, refine the possible redundant list
1036  return( delete_dublicates_noteval(result) );
1037}//proc facFirstWeyl
1038example
1039{
1040  "EXAMPLE:";echo=2;
1041  ring R = 0,(x,y),dp;
1042  def r = nc_algebra(1,1);
1043  setring(r);
1044  poly h = (x^2*y^2+x)*(x+1);
1045  facFirstWeyl(h);
1046}
1047
1048//This is the main program
1049static proc sfacwa(poly h)
1050"USAGE: sfacwa(h); h is a polynomial in the first Weyl algebra
1051RETURN: list
1052PURPOSE: Computes a factorization of a polynomial h in the first Weyl algebra
1053THEORY: @code{sfacwa} returns a list with some factorizations of the given
1054        polynomial. The possibilities of the factorization of the highest
1055        homogeneous part and those of the lowest will be merged. If during this
1056        procedure a factorization of the polynomial occurs, it will be added to
1057        the output list. For a more detailed description visit
1058        @url{http://www.math.rwth-aachen.de/\~Albert.Heinle}
1059SEE ALSO: homogfacFirstWeyl_all, homogfacFirstWeyl
1060"{//proc sfacwa
1061  int p=printlevel-voice+2; // for dbprint
1062  intvec ivm11 = intvec(-1,1);
1063  intvec iv11 = intvec(1,1);
1064  intvec iv10 = intvec(1,0);
1065  intvec iv01 = intvec(0,1);
1066  intvec iv1m1 = intvec(1,-1);
1067  if(homogwithorder(h,ivm11))
1068  {
1069    dbprint(p,"==> Given polynomial is -1,1 homogeneous. Starting homog. fac. and ret. its result");
1070    return(homogfacFirstWeyl_all(h));
1071  }
1072  def r = basering;
1073  map invo = basering,-var(1),var(2);
1074  int i; int j; int k;
1075  dbprint(p,"==> Computing the degree-limits of the factorization");
1076  intvec limits = deg(h,iv11) ,deg(h,iv10),deg(h,iv01);
1077  def prod;
1078  //end finding the limits
1079  dbprint(p,"==> Computing the maximal and the minimal homogeneous part of the give polynomial");
1080  poly maxh = jet(h,deg(h,ivm11),ivm11)-jet(h,deg(h,ivm11)-1,ivm11);
1081  poly minh = jet(h,deg(h,iv1m1),iv1m1)-jet(h,deg(h,iv1m1)-1,iv1m1);
1082  list result;
1083  list temp;
1084  list homogtemp;
1085  def invhath;
1086  def invideal;
1087  list M; list hatM;
1088  dbprint(p,"==> Factorize the maximal and the minimal homog. part of the given polynomial");
1089  list f1 = homogfacFirstWeyl_all(maxh);
1090  list f2 = homogfacFirstWeyl_all(minh);
1091  int is_equal;
1092  poly hath;
1093  dbprint(p,"==> Merging the factorizations");
1094  for (i = 1; i<=size(f1);i++)
1095  {//Merge all combinations
1096    for (j = 1; j<=size(f2); j++)
1097    {
1098      M = M + mergence(f1[i],f2[j],limits);
1099    }
1100  }//Merge all combinations
1101  dbprint(p,"==> Filtering invalid combinations");
1102  for (i = 1 ; i<= size(M); i++)
1103  {//filter valid combinations
1104    if (product(M[i]) == h)
1105    {//We have one factorization
1106      result = result + list(M[i]);
1107      M = delete(M,i);
1108      continue;
1109    }//We have one factorization
1110    else
1111    {
1112      if (deg(h,ivm11)<=deg(h-product(M[i]),ivm11))
1113      {
1114        M = delete(M,i);
1115        continue;
1116      }
1117      if (deg(h,iv1m1)<=deg(h-product(M[i]),iv1m1))
1118      {
1119        M = delete(M,i);
1120        continue;
1121      }
1122    }
1123  }//filter valid combinations
1124  dbprint(p,"==> Deleting doublicates in the resulting list");
1125  M = delete_dublicates_eval(M);
1126  dbprint(p,"==> Iterating over all possible Combinations");
1127  while(size(M)>0)
1128  {//work on the elements of M
1129    hatM = list();
1130    for(i = 1; i<=size(M); i++)
1131    {//iterate over all elements of M
1132      hath = h-product(M[i]);
1133      temp = list();
1134      //First check for common inhomogeneous factors between hath and h
1135      invhath = involution(hath,invo);
1136      invideal = std(involution(ideal(M[i][1]),invo));
1137      if (involution(NF(invhath,invideal) ,invo)==0)
1138      {//hath and h have a common factor on the left
1139        j = 1;
1140        f1 = M[i];
1141        if (j+1<=size(f1))
1142        {//Checking for more than one common factor
1143          invideal = std(involution(ideal(product(f1[1..(j+1)])),invo));
1144          while(involution(NF(invhath,invideal),invo)==0)
1145          {
1146            if (j+1<size(f1))
1147            {
1148              j++;
1149            }
1150            else
1151            {
1152              break;
1153            }
1154          }
1155        }//Checking for more than one common factor
1156        invideal = lift(involution(product(f1[1..j]),invo),invhath);
1157        f2 = list(f1[1..j])+list(involution(invideal[1,1],invo));
1158        temp = temp + merge_cf(f2,f1,limits);
1159      }//hath and h have a common factor on the left
1160      if (reduce(hath, std(ideal(M[i][size(M[i])])))==0)
1161      {//hath and h have a common factor on the right
1162        j = size(M[i]);
1163        f1 = M[i];
1164        if (j-1>0)
1165        {//Checking for more than one factor
1166          while(reduce(hath,std(ideal(product(f1[(j-1)..size(f1)]))))==0)
1167          {
1168            if (j-1>1)
1169            {
1170              j--;
1171            }
1172            else
1173            {
1174              break;
1175            }
1176          }
1177        }//Checking for more than one factor
1178        f2 = list(lift(product(f1[j..size(f1)]),hath)[1,1])+list(f1[j..size(f1)]);
1179        temp = temp + merge_cf(f2,M[i],limits);
1180      }//hath and h have a common factor on the right
1181      //and now the homogeneous
1182      maxh = jet(hath,deg(hath,ivm11),ivm11)-jet(hath,deg(hath,ivm11)-1,ivm11);
1183      minh = jet(hath,deg(hath,iv1m1),iv1m1)-jet(hath,deg(hath,iv1m1)-1,iv1m1);
1184      f1 = homogfacFirstWeyl_all(maxh);
1185      f2 = homogfacFirstWeyl_all(minh);
1186      for (j = 1; j<=size(f1);j++)
1187      {
1188        for (k=1; k<=size(f2);k++)
1189        {
1190          homogtemp = mergence(f1[j],f2[k],limits);
1191        }
1192      }
1193      for (j = 1; j<= size(homogtemp); j++)
1194      {
1195        temp = temp + mergence(homogtemp[j],M[i],limits);
1196      }
1197      for (j = 1; j<=size(temp); j++)
1198      {//filtering invalid entries in temp
1199        if(product(temp[j])==h)
1200        {//This is already a result
1201          result = result + list(temp[j]);
1202          temp = delete(temp,j);
1203          continue;
1204        }//This is already a result
1205        if (deg(hath,ivm11)<=deg(hath-product(temp[j]),ivm11))
1206        {
1207          temp = delete(temp,j);
1208          continue;
1209        }
1210      }//filtering invalid entries in temp
1211      hatM = hatM + temp;
1212    }//iterate over all elements of M
1213    M = hatM;
1214    for (i = 1; i<=size(M);i++)
1215    {//checking for complete factorizations
1216      if (h == product(M[i]))
1217      {
1218        result = result + list(M[i]);
1219        M = delete(M,i);
1220        continue;
1221      }
1222    }//checking for complete factorizations
1223    M = delete_dublicates_eval(M);
1224  }//work on the elements of M
1225  //In the case, that there is none, write a constant factor before the factor of interest.
1226  dbprint(p,"==> Done");
1227  for (i = 1 ; i<=size(result);i++)
1228  {//add a constant factor
1229    if (deg(result[i][1],iv11)!=0)
1230    {
1231      result[i] = insert(result[i],1);
1232    }
1233  }//add a constant factor
1234  dbprint(p,"==> Deleting doublicates in the output list");
1235  result = delete_dublicates_noteval(result);
1236  return(result);
1237}//proc sfacwa
1238
1239
1240//==================================================
1241/*Singular has no way implemented to test polynomials
1242  for homogenity with respect to a weight vector.
1243  The following procedure does exactly this*/
1244static proc homogwithorder(poly h, intvec weights)
1245{//proc homogwithorder
1246  if(size(weights) != nvars(basering))
1247  {//The user does not know how many variables the current ring has
1248    return(0);
1249  }//The user does not know how many variables the current ring has
1250  int i;
1251  int dofp = deg(h,weights); //degree of polynomial
1252  for (i = 1; i<=size(h);i++)
1253  {
1254    if (deg(h[i],weights)!=dofp)
1255    {
1256      return(0);
1257    }
1258  }
1259  return(1);
1260}//proc homogwithorder
1261
1262//==================================================
1263//Testfac: Given a list with different factorizations of
1264// one polynomial, the following procedure checks
1265// whether they all refer to the same polynomial.
1266// If they do, the output will be a list, that contains
1267// the product of each factorization. If not, the empty
1268// list will be returned.
1269// If the optional argument # is given (i.e. the polynomial
1270// which is factorized by the elements of the given list),
1271// then we look, if the entries are factorizations of p
1272// and if not, a list with the products subtracted by p
1273// will be returned
1274proc testNCfac(list l, list #)
1275"USAGE: testNCfac(l[,p,b]); l is a list, p is an optional poly, b is 1 or 0
1276RETURN: Case 1: No optional argument. In this case the output is 1, if the
1277                entries in the given list represent the same polynomial or 0
1278                otherwise.
1279        Case 2: One optional argument p is given. In this case it returns 1,
1280                if all the entries in l are factorizations of p, otherwise 0.
1281        Case 3: Second optional b is given. In this case a list is returned
1282                containing the difference between the product of each entry in
1283                l and p.
1284ASSUME: basering is the first Weyl algebra, the entries of l are polynomials
1285PURPOSE: Checks whether a list of factorizations contains factorizations of
1286         the same element in the first Weyl algebra
1287THEORY: @code{testNCfac} multiplies out each factorization and checks whether
1288        each factorization was a factorization of the same element.
1289@* - if there is only a list given, the output will be 0, if it
1290     does not contain factorizations of the same element. Otherwise the output
1291     will be 1.
1292@* - if there is a polynomial in the second argument, then the procedure checks
1293     whether the given list contains factorizations of this polynomial. If it
1294     does, then the output depends on the third argument. If it is not given,
1295     the procedure will check whether the factorizations in the list
1296     l are associated to this polynomial and return either 1 or 0, respectively.
1297     If the third argument is given, the output will be a list with
1298     the length of the given one and in each entry is the product of one
1299     entry in l subtracted by the polynomial.
1300EXAMPLE: example testNCfac; shows examples
1301SEE ALSO: facFirstWeyl, facSubWeyl, facFirstShift
1302"{//proc testfac
1303  int p = printlevel - voice + 2;
1304  dbprint(p,"==> Checking the input");
1305  if (size(l)==0)
1306  {//The empty list is given
1307    dbprint(p,"==> Given list was empty");
1308    return(list());
1309  }//The empty list is given
1310  if (size(#)>2)
1311  {//We want max. two optional arguments
1312    dbprint(p,"==> More than two optional arguments");
1313    return(list());
1314  }//We want max. two optional arguments
1315  dbprint(p,"==> Done");
1316  list result;
1317  int i; int j;
1318  if (size(#)==0)
1319  {//No optional argument is given
1320    dbprint(p,"==> No optional arguments");
1321    int valid = 1;
1322    for (i = size(l);i>=1;i--)
1323    {//iterate over the elements of the given list
1324      if (size(result)>0)
1325      {
1326        if (product(l[i])!=result[size(l)-i])
1327        {
1328          valid = 0;
1329          break;
1330        }
1331      }
1332      result = insert(result, product(l[i]));
1333    }//iterate over the elements of the given list
1334    return(valid);
1335  }//No optional argument is given
1336  else
1337  {
1338    dbprint(p,"==> Optional arguments are given.");
1339    int valid = 1;
1340    for (i = size(l);i>=1;i--)
1341    {//iterate over the elements of the given list
1342      if (product(l[i])!=#[1])
1343      {
1344        valid = 0;
1345      }
1346      result = insert(result, product(l[i])-#[1]);
1347    }//iterate over the elements of the given list
1348    if(size(#)==2)
1349    {
1350      dbprint(p,"==> A third argument is given. Output is a list now.");
1351      return(result);
1352    }
1353    return(valid);
1354  }
1355}//proc testfac
1356example
1357{
1358  "EXAMPLE:";echo=2;
1359  ring r = 0,(x,y),dp;
1360  def R = nc_algebra(1,1);
1361  setring R;
1362  poly h = (x^2*y^2+1)*(x^2);
1363  def t1 = facFirstWeyl(h);
1364  //fist a correct list
1365  testNCfac(t1);
1366  //now a correct list with the factorized polynomial
1367  testNCfac(t1,h);
1368  //now we put in an incorrect list without a polynomial
1369  t1[3][3] = y;
1370  testNCfac(t1);
1371  // take h as additional input
1372  testNCfac(t1,h);
1373  // take h as additional input and output list of differences
1374  testNCfac(t1,h,1);
1375}
1376//==================================================
1377//Procedure facSubWeyl:
1378//This procedure serves the purpose to compute a
1379//factorization of a given polynomial in a ring, whose subring
1380//is the first Weyl algebra. The polynomial must only contain
1381//the two arguments, which are also given by the user.
1382
1383proc facSubWeyl(poly h, X, D)
1384"USAGE:  facSubWeyl(h,x,y); h, X, D polynomials
1385RETURN: list
1386ASSUME: X and D are variables of a basering, which satisfy DX = XD +1.
1387@* That is,  they generate the copy of the first Weyl algebra in a basering.
1388@* Moreover, h is a polynomial in X and D only.
1389PURPOSE: compute factorizations of the polynomial, which depends on X and D.
1390EXAMPLE: example facSubWeyl; shows examples
1391SEE ALSO: facFirstWeyl, testNCfac, facFirstShift
1392"{
1393  int p = printlevel - voice + 2;
1394  dbprint(p,"==> Start initial Checks of the input.");
1395    // basering can be anything having a Weyl algebra as subalgebra
1396  def @r = basering;
1397  //We begin to check the input for assumptions
1398  // which are: X,D are vars of the basering,
1399  if ( (isVar(X)!=1) || (isVar(D)!=1) || (size(X)>1) || (size(D)>1) ||
1400       (leadcoef(X) != number(1)) || (leadcoef(D) != number(1)) )
1401  {
1402    ERROR("expected pure variables as generators of a subalgebra");
1403  }
1404  // Weyl algebra:
1405  poly w = D*X-X*D-1; // [D,X]=1
1406  poly u = D*X-X*D+1; // [X,D]=1
1407  if (u*w!=0)
1408  {
1409    // that is no combination gives Weyl
1410    ERROR("2nd and 3rd argument do not generate a Weyl algebra");
1411  }
1412  // one of two is correct
1413  int isReverted = 0; // Reverted Weyl if dx=xd-1 holds
1414  if (u==0)
1415  {
1416    isReverted = 1;
1417  }
1418  // else: do nothing
1419  // DONE with assumptions, Input successfully checked
1420  dbprint(p,"==> Successful");
1421  intvec lexpofX = leadexp(X);
1422  intvec lexpofD = leadexp(D);
1423  int varnumX=1;
1424  int varnumD=1;
1425  while(lexpofX[varnumX] != 1)
1426  {
1427    varnumX++;
1428  }
1429  while(lexpofD[varnumD] != 1)
1430  {
1431    varnumD++;
1432  }
1433  /* VL : to add printlevel stuff */
1434  dbprint(p,"==> Change positions of the two variables in the list, if needed");
1435  if (isReverted)
1436  {
1437    ring firstweyl = 0,(var(varnumD),var(varnumX)),dp;
1438    def Firstweyl = nc_algebra(1,1);
1439    setring Firstweyl;
1440    ideal M = 0:nvars(@r);
1441    M[varnumX]=var(2);
1442    M[varnumD]=var(1);
1443    map Q = @r,M;
1444    poly h= Q(h);
1445  }
1446  else
1447  { // that is unReverted
1448    ring firstweyl = 0,(var(varnumX),var(varnumD)),dp;
1449    def Firstweyl = nc_algebra(1,1);
1450    setring Firstweyl;
1451    poly h= imap(@r,h);
1452  }
1453  dbprint(p,"==> Done!");
1454  list result = facFirstWeyl(h);
1455  setring @r;
1456  list result;
1457  if (isReverted)
1458  {
1459    // map swap back
1460    ideal M; M[1] = var(varnumD); M[2] = var(varnumX);
1461    map S = Firstweyl, M;
1462    result = S(result);
1463  }
1464  else
1465  {
1466    // that is unReverted
1467    result = imap(Firstweyl,result);
1468  }
1469  return(result);
1470}//proc facSubWeyl
1471example
1472{
1473  "EXAMPLE:";echo=2;
1474  ring r = 0,(x,y,z),dp;
1475  matrix D[3][3]; D[1,3]=-1;
1476  def R = nc_algebra(1,D); // x,z generate Weyl subalgebra
1477  setring R;
1478  poly h = (x^2*z^2+x)*x;
1479  list fact1 = facSubWeyl(h,x,z);
1480  // compare with facFirstWeyl:
1481  ring s = 0,(z,x),dp;
1482  def S = nc_algebra(1,1); setring S;
1483  poly h = (x^2*z^2+x)*x;
1484  list fact2 = facFirstWeyl(h);
1485  map F = R,x,0,z;
1486  list fact1 = F(fact1); // it is identical to list fact2
1487  testNCfac(fact1); // check the correctness again
1488}
1489//==================================================
1490
1491//==================================================
1492//************From here: Shift-Algebra**************
1493//==================================================
1494//==================================================*
1495//one factorization of a homogeneous polynomial
1496//in the first Shift Algebra
1497static proc homogfacFirstShift(poly h)
1498{//proc homogfacFirstShift
1499  int p=printlevel-voice+2; //for dbprint
1500  def r = basering;
1501  poly hath;
1502  intvec iv01 = intvec(0,1);
1503  int i; int j;
1504  if (!homogwithorder(h,iv01))
1505  {//The given polynomial is not homogeneous
1506    ERROR("The given polynomial is not homogeneous.");
1507    return(list());
1508  }//The given polynomial is not homogeneous
1509  if (h==0)
1510  {
1511    return(list(0));
1512  }
1513  list result;
1514  int m = deg(h,iv01);
1515  dbprint(p,"==> exclude the homogeneous part of deg. 0");
1516  if (m>0)
1517  {//The degree is not zero
1518    hath = lift(var(2)^m,h)[1,1];
1519    for (i = 1; i<=m;i++)
1520    {
1521      result = result + list(var(2));
1522    }
1523  }//The degree is not zero
1524  else
1525  {//The degree is zero
1526    hath = h;
1527  }//The degree is zero
1528  ring tempRing = 0,(x),dp;
1529  setring tempRing;
1530  map thetamap = r,x,1;
1531  poly hath = thetamap(hath);
1532  dbprint(p,"==> Factorize it using commutative factorization.");
1533  list azeroresult = factorize(hath);
1534  list azeroresult_return_form;
1535  for (i = 1; i<=size(azeroresult[1]);i++)
1536  {//rewrite the result of the commutative factorization
1537    for (j = 1; j <= azeroresult[2][i];j++)
1538    {
1539      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
1540    }
1541  }//rewrite the result of the commutative factorization
1542  setring(r);
1543  map finalmap = tempRing,var(1);
1544  list tempresult = finalmap(azeroresult_return_form);
1545  result = tempresult+result;
1546  return(result);
1547}//proc homogfacFirstShift
1548
1549//==================================================
1550//Computes all possible homogeneous factorizations
1551static proc homogfacFirstShift_all(poly h)
1552{//proc HomogfacFirstShiftAll
1553  int p=printlevel-voice+2; //for dbprint
1554  intvec iv11 = intvec(1,1);
1555  if (deg(h,iv11) <= 0 )
1556  {//h is a constant
1557    return(list(list(h)));
1558  }//h is a constant
1559  def r = basering;
1560  list one_hom_fac; //stands for one homogeneous factorization
1561  int i; int j; int k;
1562  int shiftcounter;
1563  //Compute again a homogeneous factorization
1564  dbprint(p,"==> Computing one homog. factorization of the polynomial");
1565  one_hom_fac = homogfacFirstShift(h);
1566  one_hom_fac = delete(one_hom_fac,1);
1567  if (size(one_hom_fac) == 0)
1568  {//there is no homogeneous factorization or the polynomial was not homogeneous
1569    return(list());
1570  }//there is no homogeneous factorization or the polynomial was not homogeneous
1571  dbprint(p,"==> Permuting the 0-homogeneous part with the s");
1572  list result = permpp(one_hom_fac);
1573  for (i = 1; i<=size(result);i++)
1574  {
1575    shiftcounter = 0;
1576    for (j = 1; j<=size(result[i]); j++)
1577    {
1578      if (result[i][j]==var(2))
1579      {
1580        shiftcounter++;
1581      }
1582      else
1583      {
1584        result[i][j] = subst(result[i][j], var(1), var(1)-shiftcounter);
1585      }
1586    }
1587    result[i] = insert(result[i],1);
1588  }
1589  dbprint(p,"==> Deleting double entries in the resulting list");
1590  result = delete_dublicates_noteval(result);
1591  return(result);
1592}//proc HomogfacFirstShiftAll
1593
1594//==================================================
1595//factorization of the first Shift Algebra
1596proc facFirstShift(poly h)
1597"USAGE: facFirstShift(h); h a polynomial in the first shift algebra
1598RETURN: list
1599PURPOSE: compute all factorizations of a polynomial in the first shift algebra
1600THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle
1601ASSUME: basering in the first shift algebra
1602NOTE: Every entry of the output list is a list with factors for one possible factorization.
1603EXAMPLE: example facFirstShift; shows examples
1604SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl
1605"{//facFirstShift
1606  int p = printlevel - voice + 2;
1607  dbprint(p,"==> Checking the input.");
1608  if(nvars(basering)!=2)
1609  {//Our basering is the Shift algebra, but not the first
1610    ERROR("Basering is not the first shift algebra");
1611    return(list());
1612  }//Our basering is the Shift algebra, but not the first
1613  def r = basering;
1614  setring r;
1615  list LR = ringlist(r);
1616  number @n = leadcoef(LR[5][1,2]);
1617  poly @p = LR[6][1,2];
1618  if  ( @n!=number(1) )
1619  {
1620    ERROR("Basering is not the first shift algebra");
1621    return(list());
1622  }
1623  dbprint(p,"==> Done");
1624  list result = list();
1625  int i;int j; int k; int l; //counter
1626  // create a ring with the ordering which makes shift algebra
1627  // graded
1628  // def r = basering; // done before
1629  ring tempRing = LR[1][1],(x,s),(a(0,1),Dp);
1630  def tempRingnc = nc_algebra(1,s);
1631  setring r;
1632  // information on relations
1633  if (@p == -var(1)) // reverted shift algebra
1634  {
1635    dbprint(p,"==> Reverted shift algebra. Swaping variables in Ringlist");
1636    setring(tempRingnc);
1637    map transf = r, var(2), var(1);
1638    setring(r);
1639    map transfback = tempRingnc, var(2),var(1);
1640    //    result = transfback(resulttemp);
1641  }
1642  else
1643  {
1644    if ( @p == var(2)) // usual shift algebra
1645    {
1646      setring(tempRingnc);
1647      map transf = r, var(1), var(2);
1648      //    result = facshift(h);
1649      setring(r);
1650      map transfback = tempRingnc, var(1),var(2);
1651    }
1652    else
1653    {
1654      ERROR("Basering is not the first shift algebra");
1655      return(list());
1656    }
1657  }
1658  // main calls
1659  setring(tempRingnc);
1660  dbprint(p,"==> Factorize the given polynomial with the subroutine facshift");
1661  list resulttemp = facshift(transf(h));
1662  dbprint(p,"==> Successful");
1663  setring(r);
1664  result = transfback(resulttemp);
1665  dbprint(p,"==> Recursively check the found factors for reducibility.");
1666  list recursivetemp;
1667  for(i = 1; i<=size(result);i++)
1668  {//recursively factorize factors
1669    if(size(result[i])>2)
1670    {//Nontrivial factorization
1671      for (j=2;j<=size(result[i]);j++)
1672      {//Factorize every factor
1673        recursivetemp = facFirstShift(result[i][j]);
1674        if(size(recursivetemp)>1)
1675        {//we have a nontrivial factorization
1676          for(k=1; k<=size(recursivetemp);k++)
1677          {//insert factorized factors
1678            if(size(recursivetemp[k])>2)
1679            {//nontrivial
1680              result = insert(result,result[i],i);
1681              for(l = size(recursivetemp[k]);l>=2;l--)
1682              {
1683                result[i+1] = insert(result[i+1],recursivetemp[k][l],j);
1684              }
1685              result[i+1] = delete(result[i+1],j);
1686            }//nontrivial
1687          }//insert factorized factors
1688        }//we have a nontrivial factorization
1689      }//Factorize every factor
1690    }//Nontrivial factorization
1691  }//recursively factorize factors
1692  //now, refine the possible redundant list
1693  dbprint(p,"==> Done");
1694  return( delete_dublicates_noteval(result) );
1695}//facFirstShift
1696example
1697{
1698  "EXAMPLE:";echo=2;
1699  ring R = 0,(x,s),dp;
1700  def r = nc_algebra(1,s);
1701  setring(r);
1702  poly h = (s^2*x+x)*s;
1703  facFirstShift(h);
1704}
1705
1706static proc facshift(poly h)
1707"USAGE: facshift(h); h is a polynomial in the first Shift algebra
1708RETURN: list
1709PURPOSE: Computes a factorization of a polynomial h in the first Shift algebra
1710THEORY: @code{facshift} returns a list with some factorizations of the given
1711        polynomial. The possibilities of the factorization of the highest
1712        homogeneous part and those of the lowest will be merged. If during this
1713        procedure a factorization of the polynomial occurs, it will be added to
1714        the output list. For a more detailled description visit
1715        @url{http://www.math.rwth-aachen.de/\~Albert.Heinle}
1716SEE ALSO: homogfacFirstShift_all, homogfacFirstShift
1717"{//proc facshift
1718  int p=printlevel-voice+2; // for dbprint
1719  intvec iv01 = intvec(0,1);
1720  intvec iv11 = intvec(1,1);
1721  intvec iv10 = intvec(1,0);
1722  intvec iv0m1 = intvec(0,-1);
1723  if(homogwithorder(h,iv01))
1724  {
1725    dbprint(p,"==> Given polynomial was 1,0-homogeneous. Calling subroutine for homog. polys.");
1726    return(homogfacFirstShift_all(h));
1727  }
1728  def r = basering;
1729  map invo = basering,-var(1),-var(2);
1730  int i; int j; int k;
1731  intvec limits = deg(h,iv11) ,deg(h,iv10),deg(h,iv01);
1732  def prod;
1733  //end finding the limits
1734  poly maxh = jet(h,deg(h,iv01),iv01)-jet(h,deg(h,iv01)-1,iv01);
1735  poly minh = jet(h,deg(h,iv0m1),iv0m1)-jet(h,deg(h,iv0m1)-1,iv0m1);
1736  list result;
1737  list temp;
1738  list homogtemp;
1739  list M; list hatM;
1740  def invhath;
1741  def invideal;
1742  dbprint(p,"==> Factorizing the highest and the lowest homogeneous part");
1743  list f1 = homogfacFirstShift_all(maxh);
1744  list f2 = homogfacFirstShift_all(minh);
1745  int is_equal;
1746  poly hath;
1747  dbprint(p,"==> Merging their factors");
1748  for (i = 1; i<=size(f1);i++)
1749  {//Merge all combinations
1750    for (j = 1; j<=size(f2); j++)
1751    {
1752      M = M + mergence(f1[i],f2[j],limits);
1753    }
1754  }//Merge all combinations
1755  dbprint(p,"==> Filtering invalid combinations");
1756  for (i = 1 ; i<= size(M); i++)
1757  {//filter valid combinations
1758    if (product(M[i]) == h)
1759    {//We have one factorization
1760      result = result + list(M[i]);
1761      M = delete(M,i);
1762      continue;
1763    }//We have one factorization
1764    else
1765    {
1766      if (deg(h,iv01)<=deg(h-product(M[i]),iv01))
1767      {
1768        M = delete(M,i);
1769        continue;
1770      }
1771      if (deg(h,iv0m1)<=deg(h-product(M[i]),iv0m1))
1772      {
1773        M = delete(M,i);
1774        continue;
1775      }
1776    }
1777  }//filter valid combinations
1778  dbprint(p,"==> Iterating over all elements in M");
1779  M = delete_dublicates_eval(M);
1780  while(size(M)>0)
1781  {//work on the elements of M
1782    hatM = list();
1783    for(i = 1; i<=size(M); i++)
1784    {//iterate over all elements of M
1785      hath = h-product(M[i]);
1786      temp = list();
1787      //First check for common inhomogeneous factors between hath and h
1788      invhath = involution(hath,invo);
1789      invideal = std(involution(ideal(M[i][1]),invo));
1790      if (involution(NF(invhath, invideal ),invo)==0)
1791      {//hath and h have a common factor on the left
1792        j = 1;
1793        f1 = M[i];
1794        if (j+1<=size(f1))
1795        {//Checking for more than one common factor
1796          invideal = std(involution(ideal(product(f1[1..(j+1)])),invo));
1797          while(involution(NF(invhath,invideal),invo)==0)
1798          {
1799            if (j+1<size(f1))
1800            {
1801              j++;
1802            }
1803            else
1804            {
1805              break;
1806            }
1807          }
1808        }//Checking for more than one common factor
1809        if (deg(product(f1[1..j]),iv11)!=0)
1810        {
1811          invideal = lift(involution(product(f1[1..j]),invo),invhath);
1812          f2 = list(f1[1..j])+list(involution(invideal[1,1],invo));
1813        }
1814        else
1815        {
1816          invideal = lift(product(f1[1..j]),involution(hath,invo));
1817          f2 = list(f1[1..j])+list(involution(invideal[1,1],invo));
1818        }
1819        temp = temp + merge_cf(f2,f1,limits);
1820      }//hath and h have a common factor on the left
1821      if (reduce(hath, std(ideal(M[i][size(M[i])])))==0)
1822      {//hath and h have a common factor on the right
1823        j = size(M[i]);
1824        f1 = M[i];
1825        if (j-1>0)
1826        {//Checking for more than one factor
1827          while(reduce(hath,std(ideal(product(f1[(j-1)..size(f1)]))))==0)
1828          {
1829            if (j-1>1)
1830            {
1831              j--;
1832            }
1833            else
1834            {
1835              break;
1836            }
1837          }
1838        }//Checking for more than one factor
1839        f2 = list(lift(product(f1[j..size(f1)]),hath)[1,1])+list(f1[j..size(f1)]);
1840        temp = temp + merge_cf(f2,M[i],limits);
1841      }//hath and h have a common factor on the right
1842      //and now the homogeneous
1843      maxh = jet(hath,deg(hath,iv01),iv01)-jet(hath,deg(hath,iv01)-1,iv01);
1844      minh = jet(hath,deg(hath,iv0m1),iv0m1)-jet(hath,deg(hath,iv0m1)-1,iv0m1);
1845      f1 = homogfacFirstShift_all(maxh);
1846      f2 = homogfacFirstShift_all(minh);
1847      for (j = 1; j<=size(f1);j++)
1848      {
1849        for (k=1; k<=size(f2);k++)
1850        {
1851          homogtemp = mergence(f1[j],f2[k],limits);
1852        }
1853      }
1854      for (j = 1; j<= size(homogtemp); j++)
1855      {
1856        temp = temp + mergence(homogtemp[j],M[i],limits);
1857      }
1858      for (j = 1; j<=size(temp); j++)
1859      {//filtering invalid entries in temp
1860        if(product(temp[j])==h)
1861        {//This is already a result
1862          result = result + list(temp[j]);
1863          temp = delete(temp,j);
1864          continue;
1865        }//This is already a result
1866        if (deg(hath,iv01)<=deg(hath-product(temp[j]),iv01))
1867        {
1868          temp = delete(temp,j);
1869          continue;
1870        }
1871      }//filtering invalid entries in temp
1872      hatM = hatM + temp;
1873    }//iterate over all elements of M
1874    M = hatM;
1875    for (i = 1; i<=size(M);i++)
1876    {//checking for complete factorizations
1877      if (h == product(M[i]))
1878      {
1879        result = result + list(M[i]);
1880        M = delete(M,i);
1881        continue;
1882      }
1883    }//checking for complete factorizations
1884    M = delete_dublicates_eval(M);
1885  }//work on the elements of M
1886  //In the case, that there is none, write a constant factor before the factor of interest.
1887  for (i = 1 ; i<=size(result);i++)
1888  {//add a constant factor
1889    if (deg(result[i][1],iv11)!=0)
1890    {
1891      result[i] = insert(result[i],1);
1892    }
1893  }//add a constant factor
1894  dbprint(p,"==> Deleting double entries in the resulting list.");
1895  result = delete_dublicates_noteval(result);
1896  if (size(result)==0)
1897  {//only the trivial factorization could be found
1898    result = list(list(1,h));
1899  }//only the trivial factorization could be found
1900  return(result);
1901}//proc facshift
1902
1903static proc refineFactList(list L)
1904{
1905  // assume: list L is an output of factorization proc
1906  // doing: remove doubled entries
1907  int s = size(L); int sm;
1908  int i,j,k,cnt;
1909  list M, U, A, B;
1910  A = L;
1911  k = 0;
1912  cnt  = 1;
1913  for (i=1; i<=s; i++)
1914  {
1915    if (size(A[i]) != 0)
1916    {
1917      M = A[i];
1918      //      "probing with"; M; i;
1919      B[cnt] = M; cnt++;
1920      for (j=i+1; j<=s; j++)
1921      {
1922        if ( isEqualList(M,A[j]) )
1923        {
1924          k++;
1925        // U consists of intvecs with equal pairs
1926          U[k] = intvec(i,j);
1927          A[j] = 0;
1928        }
1929      }
1930    }
1931  }
1932  kill A,U,M;
1933  return(B);
1934}
1935example
1936{
1937  "EXAMPLE:";echo=2;
1938  ring R = 0,(x,s),dp;
1939  def r = nc_algebra(1,1);
1940  setring(r);
1941  list l,m;
1942  l = list(1,s2+1,x,s,x+s);
1943  m = l,list(1,s,x,s,x),l;
1944  refineFactList(m);
1945}
1946
1947static proc isEqualList(list L, list M)
1948{
1949  // int boolean: 1=yes, 0 =no : test whether two lists are identical
1950  int s = size(L);
1951  if (size(M)!=s) { return(0); }
1952  int j=1;
1953  while ( (L[j]==M[j]) && (j<s) )
1954  {
1955    j++;
1956  }
1957  if (L[j]==M[j])
1958  {
1959    return(1);
1960  }
1961  return(0);
1962}
1963example
1964{
1965  "EXAMPLE:";echo=2;
1966  ring R = 0,(x,s),dp;
1967  def r = nc_algebra(1,1);
1968  setring(r);
1969  list l,m;
1970  l = list(1,s2+1,x,s,x+s);
1971  m = l;
1972  isEqualList(m,l);
1973}
1974
1975
1976//////////////////////////////////////////////////
1977// Q-WEYL-SECTION
1978//////////////////////////////////////////////////
1979
1980//==================================================
1981//A function to get the i'th triangular number
1982proc triangNum(int n)
1983{
1984     if (n == 0)
1985     {
1986          return(0);
1987     }
1988     return (n*(n+1)/2);
1989}
1990
1991//==================================================*
1992//one factorization of a homogeneous polynomial
1993//in the first Q Weyl Algebra
1994proc homogfacFirstQWeyl(poly h)
1995"USAGE: homogfacFirstQWeyl(h); h is a homogeneous polynomial in the
1996        first q-Weyl algebra with respect to the weight vector [-1,1]
1997RETURN: list
1998PURPOSE: Computes a factorization of a homogeneous polynomial h with
1999         respect to the weight vector [-1,1] in the first q-Weyl algebra
2000THEORY: @code{homogfacFirstQWeyl} returns a list with a factorization of the given,
2001        [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with
2002        k positive, the last k entries in the output list are the second
2003        variable. If k is positive, the last k entries will be x. The other
2004        entries will be irreducible polynomials of degree zero or 1 resp. -1.
2005SEE ALSO: homogfacFirstWeyl
2006"{//proc homogfacFirstQWeyl
2007  int p = printlevel-voice+2;//for dbprint
2008  def r = basering;
2009  poly hath;
2010  int i; int j;
2011  intvec ivm11 = intvec(-1,1);
2012  if (!homogwithorder(h,ivm11))
2013  {//The given polynomial is not homogeneous
2014    ERROR("Given polynomial was not [-1,1]-homogeneous");
2015    return(list());
2016  }//The given polynomial is not homogeneous
2017  if (h==0)
2018  {
2019    return(list(0));
2020  }
2021  list result;
2022  int m = deg(h,ivm11);
2023  dbprint(p,"==> Splitting the polynomial in A_0 and A_k-Part");
2024  if (m!=0)
2025  {//The degree is not zero
2026    if (m <0)
2027    {//There are more x than y
2028      hath = lift(var(1)^(-m),h)[1,1];
2029      for (i = 1; i<=-m; i++)
2030      {
2031        result = result + list(var(1));
2032      }
2033    }//There are more x than y
2034    else
2035    {//There are more y than x
2036      hath = lift(var(2)^m,h)[1,1];
2037      for (i = 1; i<=m;i++)
2038      {
2039        result = result + list(var(2));
2040      }
2041    }//There are more y than x
2042  }//The degree is not zero
2043  else
2044  {//The degree is zero
2045    hath = h;
2046  }//The degree is zero
2047  dbprint(p,"==> Done");
2048  //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1)
2049  list mons;
2050  dbprint(p,"==> Putting the monomials in the A_0-part in a list.");
2051  for(i = 1; i<=size(hath);i++)
2052  {//Putting the monomials in a list
2053    mons = mons+list(hath[i]);
2054  }//Putting the monomials in a list
2055  dbprint(p,"==> Done");
2056  dbprint(p,"==> Mapping this monomials to K(q)[theta]");
2057  def characteristic = ringlist(r)[1][1];
2058  def qparameter      = ringlist(r)[1][2][1];
2059  ring tempRing = (characteristic,q),(x,y,theta),dp; //TODO: How to map a parameter?
2060  setring tempRing;
2061  map thetamap = r,x,y;
2062  list mons = thetamap(mons);
2063  poly entry;
2064  poly tempSummand;
2065  for (i = 1; i<=size(mons);i++)
2066  {//transforming the monomials as monomials in theta
2067       entry = 1; //leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
2068    for (j = 0; j<leadexp(mons[i])[2];j++)
2069    {
2070//"j:";j;
2071         tempSummand = (q^j-1)/(q-1);
2072         entry = entry * theta-tempSummand*entry;
2073    }
2074//~;   
2075 mons[i] = entry*leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
2076  }//transforming the monomials as monomials in theta
2077  dbprint(p,"==> Done");
2078  dbprint(p,"==> Factorize the A_0-Part in K[theta]");
2079  list azeroresult = factorize(sum(mons));
2080  dbprint(p,"==> Successful");
2081  list azeroresult_return_form;
2082  for (i = 1; i<=size(azeroresult[1]);i++)
2083  {//rewrite the result of the commutative factorization
2084    for (j = 1; j <= azeroresult[2][i];j++)
2085    {
2086      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
2087    }
2088  }//rewrite the result of the commutative factorization
2089  dbprint(p,"==> Mapping back to A_0.");
2090  setring(r);
2091  map finalmap = tempRing,var(1),var(2),var(1)*var(2);
2092  list tempresult = finalmap(azeroresult_return_form);
2093  dbprint(p,"Successful.");
2094  for (i = 1; i<=size(tempresult);i++)
2095  {//factorizations of theta resp. theta +1
2096    if(tempresult[i]==var(1)*var(2))
2097    {
2098      tempresult = insert(tempresult,var(1),i-1);
2099      i++;
2100      tempresult[i]=var(2);
2101    }
2102    if(tempresult[i]==var(2)*var(1))
2103    {
2104      tempresult = insert(tempresult,var(2),i-1);
2105      i++;
2106      tempresult[i]=var(1);
2107    }
2108  }//factorizations of theta resp. theta +1
2109  result = tempresult+result;
2110  //Correction of the result in the special q-Case:
2111  for (j = 2 ; j<= size(result);j++)
2112  {//Divide the whole Term by the leading coefficient and multiply it to the first entry in result[i]
2113       result[1] = result[1] * leadcoef(result[j]);
2114       result[j] = 1/leadcoef(result[j]) * result[j];
2115  }//Divide the whole Term by the leading coefficient and multiply it to the first entry in result[i]
2116  return(result);
2117}//proc homogfacFirstQWeyl
2118
2119
2120
2121//==================================================
2122//Computes all possible homogeneous factorizations for an element in the first Q-Weyl Algebra
2123proc homogfacFirstQWeyl_all(poly h)
2124"USAGE: homogfacFirstWWeyl_all(h); h is a homogeneous polynomial in the first q-Weyl algebra
2125        with respect to the weight vector [-1,1]
2126RETURN: list
2127PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect
2128         to the weight vector [-1,1] in the first q-Weyl algebra
2129THEORY: @code{homogfacFirstQWeyl} returns a list with all factorization of the given,
2130        homogeneous polynomial. It uses the output of homogfacFirstQWeyl and permutes
2131        its entries with respect to the commutation rule. Furthermore, if a
2132        factor of degree zero is irreducible in K[\theta], but reducible in
2133        the first q-Weyl algebra, the permutations of this element with the other
2134        entries will also be computed.
2135SEE ALSO: homogfacFirstWeyl
2136"{//proc HomogfacFirstQWeylAll
2137  int p=printlevel-voice+2;//for dbprint
2138  intvec iv11= intvec(1,1);
2139  if (deg(h,iv11) <= 0 )
2140  {//h is a constant
2141    dbprint(p,"Given polynomial was not homogeneous");
2142    return(list(list(h)));
2143  }//h is a constant
2144  def r = basering;
2145  list one_hom_fac; //stands for one homogeneous factorization
2146  int i; int j; int k;
2147  intvec ivm11 = intvec(-1,1);
2148  dbprint(p,"==> Calculate one homogeneous factorization using homogfacFirstQWeyl");
2149  //Compute again a homogeneous factorization
2150  one_hom_fac = homogfacFirstQWeyl(h);
2151  dbprint(p,"Successful");
2152  if (size(one_hom_fac) == 0)
2153  {//there is no homogeneous factorization or the polynomial was not homogeneous
2154    return(list());
2155  }//there is no homogeneous factorization or the polynomial was not homogeneous
2156  //divide list in A0-Part and a list of x's resp. y's
2157  list list_not_azero = list();
2158  list list_azero;
2159  list k_factor;
2160  int is_list_not_azero_empty = 1;
2161  int is_list_azero_empty = 1;
2162  k_factor = list(one_hom_fac[1]);
2163  if (absValue(deg(h,ivm11))<size(one_hom_fac)-1)
2164  {//There is a nontrivial A_0-part
2165    list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,ivm11)))];
2166    is_list_azero_empty = 0;
2167  }//There is a nontrivial A_0 part
2168  dbprint(p,"==> Combine x,y to xy in the factorization again.");
2169  for (i = 1; i<=size(list_azero)-1;i++)
2170  {//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
2171    if (list_azero[i] == var(1))
2172    {
2173      if (list_azero[i+1]==var(2))
2174      {
2175        list_azero[i] = var(1)*var(2);
2176        list_azero = delete(list_azero,i+1);
2177      }
2178    }
2179    if (list_azero[i] == var(2))
2180    {
2181      if (list_azero[i+1]==var(1))
2182      {
2183        list_azero[i] = var(2)*var(1);
2184        list_azero = delete(list_azero,i+1);
2185      }
2186    }
2187  }//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
2188  dbprint(p,"==> Done");
2189  if(deg(h,ivm11)!=0)
2190  {//list_not_azero is not empty
2191    list_not_azero =
2192      one_hom_fac[(size(one_hom_fac)-absValue(deg(h,ivm11))+1)..size(one_hom_fac)];
2193    is_list_not_azero_empty = 0;
2194  }//list_not_azero is not empty
2195  //Map list_azero in K[theta]
2196  dbprint(p,"==> Map list_azero to K[theta]");
2197  def characteristic = ringlist(r)[1][1];
2198  def qparameter      = ringlist(r)[1][2][1];
2199  ring tempRing = (characteristic,q),(x,y,theta),dp; //TODO: How to map a parameter?
2200  setring(tempRing);
2201  poly entry;
2202  map thetamap = r,x,y;
2203  if(!is_list_not_azero_empty)
2204  {//Mapping in Singular is only possible, if the list before
2205    //contained at least one element of the other ring
2206    list list_not_azero = thetamap(list_not_azero);
2207  }//Mapping in Singular is only possible, if the list before
2208  //contained at least one element of the other ring
2209  if(!is_list_azero_empty)
2210  {//Mapping in Singular is only possible, if the list before
2211    //contained at least one element of the other ring
2212    list list_azero= thetamap(list_azero);
2213  }//Mapping in Singular is only possible, if the list before
2214  //contained at least one element of the other ring
2215  list k_factor = thetamap(k_factor);
2216  list tempmons;
2217  dbprint(p,"==> Done");
2218  for(i = 1; i<=size(list_azero);i++)
2219  {//rewrite the polynomials in A1 as polynomials in K[theta]
2220    tempmons = list();
2221    for (j = 1; j<=size(list_azero[i]);j++)
2222    {
2223      tempmons = tempmons + list(list_azero[i][j]);
2224    }
2225    for (j = 1 ; j<=size(tempmons);j++)
2226    {
2227         //entry = leadcoef(tempmons[j]);
2228         entry = leadcoef(tempmons[j]) * q^(-triangNum(leadexp(tempmons[j])[2]-1));
2229      for (k = 0; k < leadexp(tempmons[j])[2];k++)
2230      {
2231           entry = entry*(theta-(q^k-1)/(q-1));
2232      }
2233      tempmons[j] = entry;
2234    }
2235    list_azero[i] = sum(tempmons);
2236  }//rewrite the polynomials in A1 as polynomials in K[theta]
2237  //Compute all permutations of the A0-part
2238  dbprint(p,"==> Compute all permutations of the A_0-part with the first resp. the snd. variable");
2239  list result;
2240  int shift_sign;
2241  int shift;
2242  poly shiftvar;
2243  if (size(list_not_azero)!=0)
2244  {//Compute all possibilities to permute the x's resp. the y's in the list
2245    if (list_not_azero[1] == x)
2246    {//h had a negative weighted degree
2247      shift_sign = 1;
2248      shiftvar = x;
2249    }//h had a negative weighted degree
2250    else
2251    {//h had a positive weighted degree
2252      shift_sign = -1;
2253      shiftvar = y;
2254    }//h had a positive weighted degree
2255    result = permpp(list_azero + list_not_azero);
2256    for (i = 1; i<= size(result); i++)
2257    {//adjust the a_0-parts
2258      shift = 0;
2259      for (j=1; j<=size(result[i]);j++)
2260      {
2261        if (result[i][j]==shiftvar)
2262        {
2263          shift = shift + shift_sign;
2264        }
2265        else
2266        {
2267             if (shift < 0)
2268             {//We have two distict formulas for x and y. In this case use formula for y
2269                  if (shift == -1)
2270                  {
2271                       result[i][j] = subst(result[i][j],theta,1/q*(theta - 1));
2272                  }
2273                  else
2274                  {
2275                       result[i][j] = subst(result[i][j],theta,1/q*((theta - 1)/q^(absValue(shift)-1) - (q^(shift +2)-q)/(1-q)));
2276                  }
2277             }//We have two distict formulas for x and y. In this case use formula for y
2278             if (shift > 0)
2279             {//We have two distict formulas for x and y. In this case use formula for x
2280                  if (shift == 1)
2281                  {
2282                       result[i][j] = subst(result[i][j],theta,q*theta + 1);
2283                  }
2284                  else
2285                  {
2286                       result[i][j] = subst(result[i][j],theta,q^shift*theta+(q^shift-1)/(q-1));
2287                  }
2288             }//We have two distict formulas for x and y. In this case use formula for x
2289        }
2290      }
2291    }//adjust the a_0-parts
2292  }//Compute all possibilities to permute the x's resp. the y's in the list
2293  else
2294  {//The result is just all the permutations of the a_0-part
2295    result = permpp(list_azero);
2296  }//The result is just all the permutations of the a_0 part
2297  if (size(result)==0)
2298  {
2299    return(result);
2300  }
2301  dbprint(p,"==> Done");
2302  dbprint(p,"==> Searching for theta resp. theta + 1 in the list and factorize them");
2303  //Now we are going deeper and search for theta resp. theta + 1, substitute
2304  //them by xy resp. yx and go on permuting
2305  int found_theta;
2306  int thetapos;
2307  list leftpart;
2308  list rightpart;
2309  list lparts;
2310  list rparts;
2311  list tempadd;
2312  for (i = 1; i<=size(result) ; i++)
2313  {//checking every entry of result for theta or theta +1
2314    found_theta = 0;
2315    for(j=1;j<=size(result[i]);j++)
2316    {
2317      if (result[i][j]==theta)
2318      {//the jth entry is theta and can be written as x*y
2319        thetapos = j;
2320        result[i]= insert(result[i],x,j-1);
2321        j++;
2322        result[i][j] = y;
2323        found_theta = 1;
2324        break;
2325      }//the jth entry is theta and can be written as x*y
2326      if(result[i][j] == q*theta +1)
2327      {
2328        thetapos = j;
2329        result[i] = insert(result[i],y,j-1);
2330        j++;
2331        result[i][j] = x;
2332        found_theta = 1;
2333        break;
2334      }
2335    }
2336    if (found_theta)
2337    {//One entry was theta resp. theta +1
2338      leftpart = result[i];
2339      leftpart = leftpart[1..thetapos];
2340      rightpart = result[i];
2341      rightpart = rightpart[(thetapos+1)..size(rightpart)];
2342      lparts = list(leftpart);
2343      rparts = list(rightpart);
2344      //first deal with the left part
2345      if (leftpart[thetapos] == x)
2346      {
2347        shift_sign = 1;
2348        shiftvar = x;
2349      }
2350      else
2351      {
2352        shift_sign = -1;
2353        shiftvar = y;
2354      }
2355      for (j = size(leftpart); j>1;j--)
2356      {//drip x resp. y
2357        if (leftpart[j-1]==shiftvar)
2358        {//commutative
2359          j--;
2360          continue;
2361        }//commutative
2362        if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
2363        {//stop here
2364          break;
2365        }//stop here
2366        //Here, we can only have a a0- part
2367        if (shift_sign<0)
2368        {
2369             leftpart[j] = subst(leftpart[j-1],theta, 1/q*(theta +shift_sign));
2370        }
2371        if (shift_sign>0)
2372        {
2373             leftpart[j] = subst(leftpart[j-1],theta, q*theta + shift_sign);
2374        }
2375        leftpart[j-1] = shiftvar;
2376        lparts = lparts + list(leftpart);
2377      }//drip x resp. y
2378      //and now deal with the right part
2379      if (rightpart[1] == x)
2380      {
2381        shift_sign = 1;
2382        shiftvar = x;
2383      }
2384      else
2385      {
2386        shift_sign = -1;
2387        shiftvar = y;
2388      }
2389      for (j = 1 ; j < size(rightpart); j++)
2390      {
2391        if (rightpart[j+1] == shiftvar)
2392        {
2393          j++;
2394          continue;
2395        }
2396        if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
2397        {
2398          break;
2399        }
2400        if (shift_sign<0)
2401        {
2402             rightpart[j] = subst(rightpart[j+1], theta, q*theta - shift_sign);
2403        }
2404        if (shift_sign>0)
2405        {
2406             rightpart[j] = subst(rightpart[j+1], theta, 1/q*(theta - shift_sign));
2407        }
2408        rightpart[j+1] = shiftvar;
2409        rparts = rparts + list(rightpart);
2410      }
2411      //And now, we put all possibilities together
2412      tempadd = list();
2413      for (j = 1; j<=size(lparts); j++)
2414      {
2415        for (k = 1; k<=size(rparts);k++)
2416        {
2417          tempadd = tempadd + list(lparts[j]+rparts[k]);
2418        }
2419      }
2420      tempadd = delete(tempadd,1); // The first entry is already in the list
2421      result = result + tempadd;
2422      continue; //We can may be not be done already with the ith entry
2423    }//One entry was theta resp. theta +1
2424  }//checking every entry of result for theta or theta +1
2425  dbprint(p,"==> Done");
2426  //map back to the basering
2427  dbprint(p,"==> Mapping back everything to the basering");
2428  setring(r);
2429  map finalmap = tempRing, var(1), var(2),var(1)*var(2);
2430  list result = finalmap(result);
2431  for (i=1; i<=size(result);i++)
2432  {//adding the K factor
2433    result[i] = k_factor + result[i];
2434  }//adding the k-factor
2435  dbprint(p,"==> Done");
2436  dbprint(p,"==> Delete double entries in the list.");
2437  result = delete_dublicates_noteval(result);
2438  dbprint(p,"==> Done");
2439  return(result);
2440}//proc HomogfacFirstQWeylAll
2441
2442//TODO: FirstQWeyl check the parameters...
2443
2444/*
2445  Example polynomials where one can find factorizations: K<x,y |yx=xy+1>
2446  (x^2+y)*(x^2+y);
2447  (x^2+x)*(x^2+y);
2448  (x^3+x+1)*(x^4+y*x+2);
2449  (x^2*y+y)*(y+x*y);
2450  y^3+x*y^3+2*y^2+2*(x+1)*y^2+y+(x+2)*y; //Example 5 Grigoriev-Schwarz.
2451  (y+1)*(y+1)*(y+x*y); //Landau Example projected to the first dimension.
2452*/
2453
2454
2455/* very hard things from Martin Lee:
2456// ex1, ex2
2457ring s = 0,(z,x),Ws(-1,1);
2458def S = nc_algebra(1,1); setring S;
2459poly a = 10z5x4+26z4x5+47z5x2-97z4x3; //Abgebrochen nach einer Stunde; yes, it takes long
2460def l= facFirstWeyl (a); l;
2461kill l;
2462poly b = -5328z8x5-5328z7x6+720z9x2+720z8x3-16976z7x4-38880z6x5-5184z7x3-5184z6x4-3774z5x5+2080z8x+5760z7x2-6144z6x3-59616z5x4+3108z3x6-4098z6x2-25704z5x3-21186z4x4+8640z6x-17916z4x3+22680z2x5+2040z5x-4848z4x2-9792z3x3+3024z2x4-10704z3x2-3519z2x3+34776zx4+12096zx3+2898x4-5040z2x+8064x3+6048x2; //Abgebrochen nach 1.5 Stunden; seems to be very complicated
2463def l= facFirstWeyl (b); l;
2464
2465// ex3: there was difference in answers => fixed
2466LIB "ncfactor.lib";
2467ring r = 0,(x,y,z),dp;
2468matrix D[3][3]; D[1,3]=-1;
2469def R = nc_algebra(1,D);
2470setring R;
2471poly g= 7*z4*x+62*z3+26*z;
2472def l1= facSubWeyl (g, x, z);
2473l1;
2474//---- other ring
2475ring s = 0,(x,z),dp;
2476def S = nc_algebra(1,-1); setring S;
2477poly g= 7*z4*x+62*z3+26*z;
2478def l2= facFirstWeyl (g);
2479l2;
2480map F = R,x,0,z;
2481list l1 = F(l1);
2482l1;
2483//---- so the answers look different, check them!
2484testNCfac(l2); // ok
2485testNCfac(l1); // was not ok, but now it's been fixed!!!
2486
2487// selbst D und X so vertauschen dass sie erfuellt ist : ist gemacht
2488
2489*/
2490
2491/*
2492// bug from M Lee
2493LIB "ncfactor.lib";
2494ring s = 0,(z,x),dp;
2495def S = nc_algebra(1,1); setring S;
2496poly f= -60z4x2-54z4-56zx3-59z2x-64;
2497def l= facFirstWeyl (f);
2498l; // before: empty list; after fix: 1 entry, f is irreducible
2499poly g = 75z3x2+92z3+24;
2500def l= facFirstWeyl (g);
2501l; //before: empty list, now: correct
2502*/
2503
2504/* more things from Martin Lee; fixed
2505ring R = 0,(x,s),dp;
2506def r = nc_algebra(1,s);
2507setring(r);
2508poly h = (s2*x+x)*s;
2509h= h* (x+s);
2510def l= facFirstShift(h);
2511l; // contained doubled entries: not anymore, fixed!
2512
2513ring R = 0,(x,s),dp;
2514def r = nc_algebra(1,-1);
2515setring(r);
2516poly h = (s2*x+x)*s;
2517h= h* (x+s);
2518def l= facFirstWeyl(h);
2519l; // contained doubled entries: not anymore, fixed!
2520
2521*/
Note: See TracBrowser for help on using the repository browser.