source: git/Singular/LIB/ncfactor.lib @ e1b97a7

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