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

spielwiese
Last change on this file since ee03c0 was ee03c0, checked in by Martin Lee <martinlee84@…>, 11 years ago
fixed version of ncfactor.lib courtesy of aheinle@uwaterloo.ca
  • Property mode set to 100644
File size: 153.1 KB
Line 
1//////////////////////////////////////////////////////////
2version="version ncfactor.lib 4.0.0.0 Jun_2013 "; // $Id$
3category="Noncommutative";
4info="
5LIBRARY: ncfactor.lib  Tools for factorization in some noncommutative algebras
6AUTHORS: Albert Heinle,     aheinle@uwaterloo.ca
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[,1]]);      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  homogfacFirstQWeyl(h);     [-1,1]-homogeneous factorization in the first Q-Weyl algebra
27  homogfacFirstQWeyl_all(h); [-1,1] homogeneous factorization(complete) in the first Q-Weyl algebra
28  tst_ncfactor();            Runs the examples of all contained not static functions. Test thing.
29";
30
31LIB "general.lib";
32LIB "nctools.lib";
33LIB "involut.lib";
34LIB "freegb.lib"; // for isVar
35LIB "crypto.lib"; //for introot
36LIB "matrix.lib"; //for submatrix
37LIB "solve.lib"; //right now not needed
38LIB "poly.lib"; //for content
39
40proc tst_ncfactor()
41"
42A little test if the library works correct.
43Runs simply all examples of non-static functions.
44"
45{
46  example facFirstWeyl;
47  example facFirstShift;
48  example facSubWeyl;
49  example testNCfac;
50  example homogfacFirstQWeyl;
51  example homogfacFirstQWeyl_all;
52}
53example
54{
55  "EXAMPLE:";echo=2;
56   tst_ncfactor();
57}
58
59/////////////////////////////////////////////////////
60//==================================================*
61//deletes double-entries in a list of factorization
62//without evaluating the product.
63static proc delete_dublicates_noteval(list l)
64"
65INPUT: A list of lists; Output same as e.g. FacFirstWeyl. Containing different factorizations
66       of a polynomial
67OUTPUT: If there are dublicates in this list, this procedure deletes them and returns the list
68 without double entries
69"
70{//proc delete_dublicates_noteval
71  list result= l;
72  int j; int k; int i;
73  int deleted = 0;
74  int is_equal;
75  for (i = 1; i<= size(l); i++)
76  {//Iterate over the different factorizations
77    for (j = i+1; j<= size(l); j++)
78    {//Compare the i'th factorization to the j'th
79      if (size(l[i])!= size(l[j]))
80      {//different sizes => not equal
81        j++;
82        continue;
83      }//different sizes => not equal
84      is_equal = 1;
85      for (k = 1; k <= size(l[i]);k++)
86      {//Compare every entry
87        if (l[i][k]!=l[j][k])
88        {
89          is_equal = 0;
90          break;
91        }
92      }//Compare every entry
93      if (is_equal == 1)
94      {//Delete this entry, because there is another equal one int the list
95        result = delete(result, i-deleted);
96        deleted = deleted+1;
97        break;
98      }//Delete this entry, because there is another equal one int the list
99    }//Compare the i'th factorization to the j'th
100  }//Iterate over the different factorizations
101  return(result);
102}//proc delete_dublicates_noteval
103
104//==================================================
105//deletes the double-entries in a list with
106//evaluating the products
107static proc delete_dublicates_eval(list l)
108"
109DEPRECATED
110"
111{//proc delete_dublicates_eval
112  list result=l;
113  int j; int k; int i;
114  int deleted = 0;
115  int is_equal;
116  for (i = 1; i<= size(result); i++)
117  {//Iterating over all elements in result
118    for (j = i+1; j<= size(result); j++)
119    {//comparing with the other elements
120      if (product(result[i]) == product(result[j]))
121      {//There are two equal results; throw away that one with the smaller size
122        if (size(result[i])>=size(result[j]))
123        {//result[i] has more entries
124          result = delete(result,j);
125          continue;
126        }//result[i] has more entries
127        else
128        {//result[j] has more entries
129          result = delete(result,i);
130          i--;
131          break;
132        }//result[j] has more entries
133      }//There are two equal results; throw away that one with the smaller size
134    }//comparing with the other elements
135  }//Iterating over all elements in result
136  return(result);
137}//proc delete_dublicates_eval
138
139
140//==================================================*
141
142static proc combinekfinlf(list g, int nof) //nof stands for "number of factors"
143"
144given a list of factors g and a desired size nof, this
145procedure combines the factors, such that we recieve a
146list of the length nof.
147INPUT: A list of containing polynomials or any type where the *-operator is existent
148OUTPUT: All possibilities (without permutation of the given list) to combine the polynomials
149 into nof polynomials given by the user.
150"
151{//Procedure combinekfinlf
152  list result;
153  int i; int j; int k; //iteration variables
154  list fc; //fc stands for "factors combined"
155  list temp; //a temporary store for factors
156  def nofgl = size(g); //nofgl stands for "number of factors of the given list"
157  if (nofgl == 0)
158  {//g was the empty list
159    return(result);
160  }//g was the empty list
161  if (nof <= 0)
162  {//The user wants to recieve a negative number or no element as a result
163    return(result);
164  }//The user wants to recieve a negative number or no element as a result
165  if (nofgl == nof)
166  {//There are no factors to combine
167    result = result + list(g);
168    return(result);
169  }//There are no factors to combine
170  if (nof == 1)
171  {//User wants to get just one factor
172    result = result + list(list(product(g)));
173    return(result);
174  }//User wants to get just one factor
175  for (i = nof; i > 1; i--)
176  {//computing the possibilities that have at least one original factor from g
177    for (j = i; j>=1; j--)
178    {//shifting the window of combinable factors to the left
179      //fc below stands for "factors combined"
180      fc = combinekfinlf(list(g[(j)..(j+nofgl - i)]),nof - i + 1);
181      for (k = 1; k<=size(fc); k++)
182      {//iterating over the different solutions of the smaller problem
183        if (j>1)
184        {//There are g_i before the combination
185          if (j+nofgl -i < nofgl)
186          {//There are g_i after the combination
187            temp = list(g[1..(j-1)]) + fc[k] + list(g[(j+nofgl-i+1)..nofgl]);
188          }//There are g_i after the combination
189          else
190          {//There are no g_i after the combination
191            temp = list(g[1..(j-1)]) + fc[k];
192          }//There are no g_i after the combination
193        }//There are g_i before the combination
194        if (j==1)
195        {//There are no g_i before the combination
196          if (j+ nofgl -i <nofgl)
197          {//There are g_i after the combination
198            temp = fc[k]+ list(g[(j + nofgl - i +1)..nofgl]);
199          }//There are g_i after the combination
200        }//There are no g_i before the combination
201        result = result + list(temp);
202      }//iterating over the different solutions of the smaller problem
203    }//shifting the window of combinable factors to the left
204  }//computing the possibilities that have at least one original factor from g
205  for (i = 2; i<=nofgl div nof;i++)
206  {//getting the other possible results
207    result = result + combinekfinlf(list(product(list(g[1..i])))+list(g[(i+1)..nofgl]),nof);
208  }//getting the other possible results
209  result = delete_dublicates_noteval(result);
210  return(result);
211}//Procedure combinekfinlf
212
213
214//==================================================*
215//merges two sets of factors ignoring common
216//factors
217static proc merge_icf(list l1, list l2, intvec limits)
218"
219DEPRECATED
220"
221{//proc merge_icf
222  list g;
223  list f;
224  int i; int j;
225  if (size(l1)==0)
226  {
227    return(list());
228  }
229  if (size(l2)==0)
230  {
231    return(list());
232  }
233  if (size(l2)<=size(l1))
234  {//l1 will be our g, l2 our f
235    g = l1;
236    f = l2;
237  }//l1 will be our g, l2 our f
238  else
239  {//l1 will be our f, l2 our g
240    g = l2;
241    f = l1;
242  }//l1 will be our f, l2 our g
243  def result = combinekfinlf(g,size(f),limits);
244  for (i = 1 ; i<= size(result); i++)
245  {//Adding the factors of f to every possibility listed in temp
246    for (j = 1; j<= size(f); j++)
247    {
248      result[i][j] = result[i][j]+f[j];
249    }
250    if(!limitcheck(result[i],limits))
251    {
252      result = delete(result,i);
253      continue;
254    }
255    for (j = 1; j<=size(f);j++)
256    {//Delete entry if there is a zero or an integer as a factor
257      if (deg(result[i][j]) <= 0)
258      {//found one
259        result = delete(result,i);
260        i--;
261        break;
262      }//found one
263    }//Delete entry if there is a zero as factor
264  }//Adding the factors of f to every possibility listed in temp
265  return(result);
266}//proc merge_icf
267
268//==================================================*
269//merges two sets of factors with respect to the occurrence
270//of common factors
271static proc merge_cf(list l1, list l2, intvec limits)
272"
273DEPRECATED
274"
275{//proc merge_cf
276  list g;
277  list f;
278  int i; int j;
279  list pre;
280  list post;
281  list candidate;
282  list temp;
283  int temppos;
284  if (size(l1)==0)
285  {//the first list is empty
286    return(list());
287  }//the first list is empty
288  if(size(l2)==0)
289  {//the second list is empty
290    return(list());
291  }//the second list is empty
292  if (size(l2)<=size(l1))
293  {//l1 will be our g, l2 our f
294    g = l1;
295    f = l2;
296  }//l1 will be our g, l2 our f
297  else
298  {//l1 will be our f, l2 our g
299    g = l2;
300    f = l1;
301  }//l1 will be our f, l2 our g
302  list M;
303  for (i = 2; i<size(f); i++)
304  {//finding common factors of f and g...
305    for (j=2; j<size(g);j++)
306    {//... with g
307      if (f[i] == g[j])
308      {//we have an equal pair
309        M = M + list(list(i,j));
310      }//we have an equal pair
311    }//... with g
312  }//finding common factors of f and g...
313  if (g[1]==f[1])
314  {//Checking for the first elements to be equal
315    M = M + list(list(1,1));
316  }//Checking for the first elements to be equal
317  if (g[size(g)]==f[size(f)])
318  {//Checking for the last elements to be equal
319    M = M + list(list(size(f),size(g)));
320  }//Checking for the last elements to be equal
321  list result;//= list(list());
322  while(size(M)>0)
323  {//set of equal pairs is not empty
324    temp = M[1];
325    temppos = 1;
326    for (i = 2; i<=size(M); i++)
327    {//finding the minimal element of M
328      if (M[i][1]<=temp[1])
329      {//a possible candidate that is smaller than temp could have been found
330        if (M[i][1]==temp[1])
331        {//In this case we must look at the second number
332          if (M[i][2]< temp[2])
333          {//the candidate is smaller
334            temp = M[i];
335            temppos = i;
336          }//the candidate is smaller
337        }//In this case we must look at the second number
338        else
339        {//The candidate is definately smaller
340          temp = M[i];
341          temppos = i;
342        }//The candidate is definately smaller
343      }//a possible candidate that is smaller than temp could have been found
344    }//finding the minimal element of M
345    M = delete(M, temppos);
346    if(temp[1]>1)
347    {//There are factors to combine before the equal factor
348      if (temp[1]<size(f))
349      {//The most common case
350        //first the combinations ignoring common factors
351        pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
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(pre); i++)
354        {//all possible pre's...
355          for (j = 1; j<= size(post); j++)
356          {//...combined with all possible post's
357            candidate = pre[i]+list(f[temp[1]])+post[j];
358            if (limitcheck(candidate,limits))
359            {
360              result = result + list(candidate);
361            }
362          }//...combined with all possible post's
363        }//all possible pre's...
364        //Now the combinations with respect to common factors
365        post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
366        if (size(post)>0)
367        {//There are factors to combine
368          for (i = 1; i <= size(pre); i++)
369          {//all possible pre's...
370            for (j = 1; j<= size(post); j++)
371            {//...combined with all possible post's
372              candidate= pre[i]+list(f[temp[1]])+post[j];
373              if (limitcheck(candidate,limits))
374              {
375                result = result + list(candidate);
376              }
377            }//...combined with all possible post's
378          }//all possible pre's...
379        }//There are factors to combine
380      }//The most common case
381      else
382      {//the last factor is the common one
383        pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
384        for (i = 1; i<= size(pre); i++)
385        {//iterating over the possible pre-factors
386          candidate = pre[i]+list(f[temp[1]]);
387          if (limitcheck(candidate,limits))
388          {
389            result = result + list(candidate);
390          }
391        }//iterating over the possible pre-factors
392      }//the last factor is the common one
393    }//There are factors to combine before the equal factor
394    else
395    {//There are no factors to combine before the equal factor
396      if (temp[1]<size(f))
397      {//Just a check for security
398        //first without common factors
399        post=merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
400        for (i = 1; i<=size(post); i++)
401        {
402          candidate = list(f[temp[1]])+post[i];
403          if (limitcheck(candidate,limits))
404          {
405            result = result + list(candidate);
406          }
407        }
408        //Now with common factors
409        post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
410        if(size(post)>0)
411        {//we could find other combinations
412          for (i = 1; i<=size(post); i++)
413          {
414            candidate = list(f[temp[1]])+post[i];
415            if (limitcheck(candidate,limits))
416            {
417              result = result + list(candidate);
418            }
419          }
420        }//we could find other combinations
421      }//Just a check for security
422    }//There are no factors to combine before the equal factor
423  }//set of equal pairs is not empty
424  for (i = 1; i <= size(result); i++)
425  {//delete those combinations, who have an entry with degree less or equal 0
426    for (j = 1; j<=size(result[i]);j++)
427    {//Delete entry if there is a zero or an integer as a factor
428      if (deg(result[i][j]) <= 0)
429      {//found one
430        result = delete(result,i);
431        i--;
432        break;
433      }//found one
434    }//Delete entry if there is a zero as factor
435  }//delete those combinations, who have an entry with degree less or equal 0
436  return(result);
437}//proc merge_cf
438
439
440//==================================================*
441//merges two sets of factors
442
443static proc mergence(list l1, list l2, intvec limits)
444"
445DEPRECATED
446"
447{//Procedure mergence
448  list g;
449  list f;
450  int k; int i; int j;
451  list F = list();
452  list G = list();
453  list tempEntry;
454  list comb;
455  if (size(l2)<=size(l1))
456  {//l1 will be our g, l2 our f
457    g = l1;
458    f = l2;
459  }//l1 will be our g, l2 our f
460  else
461  {//l1 will be our f, l2 our g
462    g = l2;
463    f = l1;
464  }//l1 will be our f, l2 our g
465  if (size(f)==1 or size(g)==1)
466  {//One of them just has one entry
467    if (size(f)== 1) {f = list(1) + f;}
468    if (size(g) == 1) {g = list(1) + g;}
469  }//One of them just has one entry
470  //first, we need to add some latent -1's to the list f and to the list g in order
471  //to get really all possibilities of combinations later
472  for (i=1;i<=size(f)-1;i++)
473  {//first iterator
474    for (j=i+1;j<=size(f);j++)
475    {//second iterator
476      tempEntry = f;
477      tempEntry[i] = (-1)*tempEntry[i];
478      tempEntry[j] = (-1)*tempEntry[j];
479      F = F + list(tempEntry);
480    }//secont iterator
481  }//first iterator
482  F = F + list(f);
483  //And now same game with g
484  for (i=1;i<=size(g)-1;i++)
485  {//first iterator
486    for (j=i+1;j<=size(g);j++)
487    {//second iterator
488      tempEntry = g;
489      tempEntry[i] = (-1)*tempEntry[i];
490      tempEntry[j] = (-1)*tempEntry[j];
491      G = G + list(tempEntry);
492    }//secont iterator
493  }//first iterator
494  G = G + list(g);
495  //Done with that
496
497  list result;
498  for (i = 1; i<=size(F); i++)
499  {//Iterate over all entries in F
500    for (j = 1;j<=size(G);j++)
501    {//Same with G
502      comb = combinekfinlf(F[i],2,limits);
503      for (k = 1; k<= size(comb);k++)
504      {//for all possibilities of combinations of the factors of f
505        result = result + merge_cf(comb[k],G[j],limits);
506        result = result + merge_icf(comb[k],G[j],limits);
507        result = delete_dublicates_noteval(result);
508      }//for all possibilities of combinations of the factors of f
509    }//Same with G
510  }//Iterate over all entries in F
511  return(result);
512}//Procedure mergence
513
514
515//==================================================
516//Checks, whether a list of factors doesn't exceed the given limits
517static proc limitcheck(list g, intvec limits)
518"
519DEPRECATED
520"
521{//proc limitcheck
522  int i;
523  if (size(limits)!=3)
524  {//check the input
525    return(0);
526  }//check the input
527  if(size(g)==0)
528  {
529    return(0);
530  }
531  def prod = product(g);
532  intvec iv11 = intvec(1,1);
533  intvec iv10 = intvec(1,0);
534  intvec iv01 = intvec(0,1);
535  def limg = intvec(deg(prod,iv11) ,deg(prod,iv10),deg(prod,iv01));
536  for (i = 1; i<=size(limg);i++)
537  {//the final check
538    if(limg[i]>limits[i])
539    {
540      return(0);
541    }
542  }//the final check
543  return(1);
544}//proc limitcheck
545
546
547//==================================================*
548//one factorization of a homogeneous polynomial
549//in the first Weyl Algebra
550static proc homogfacFirstWeyl(poly h)
551"USAGE: homogfacFirstWeyl(h); h is a homogeneous polynomial in the
552 first Weyl algebra with respect to the weight vector [-1,1]
553RETURN: list
554PURPOSE: Computes a factorization of a homogeneous polynomial h with
555  respect to the weight vector [-1,1] in the first Weyl algebra
556THEORY: @code{homogfacFirstWeyl} returns a list with a factorization of the given,
557 [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with
558 k positive, the last k entries in the output list are the second
559 variable. If k is positive, the last k entries will be x. The other
560 entries will be irreducible polynomials of degree zero or 1 resp. -1.
561SEE ALSO: homogfacFirstWeyl_all
562"{//proc homogfacFirstWeyl
563  int p = printlevel-voice+2;//for dbprint
564  def r = basering;
565  poly hath;
566  int i; int j;
567  string dbprintWhitespace = "";
568  for (i = 1; i<=voice;i++)
569  {dbprintWhitespace = dbprintWhitespace + " ";}
570  intvec ivm11 = intvec(-1,1);
571  if (!homogwithorder(h,ivm11))
572  {//The given polynomial is not homogeneous
573    ERROR("Given polynomial was not [-1,1]-homogeneous");
574    return(list());
575  }//The given polynomial is not homogeneous
576  if (h==0)
577  {
578    return(list(0));
579  }
580  list result;
581  int m = deg(h,ivm11);
582  dbprint(p,dbprintWhitespace +" Splitting the polynomial in A_0 and A_k-Part");
583  if (m!=0)
584  {//The degree is not zero
585    if (m <0)
586    {//There are more x than y
587      hath = lift(var(1)^(-m),h)[1,1];
588      for (i = 1; i<=-m; i++)
589      {
590        result = result + list(var(1));
591      }
592    }//There are more x than y
593    else
594    {//There are more y than x
595      hath = lift(var(2)^m,h)[1,1];
596      for (i = 1; i<=m;i++)
597      {
598        result = result + list(var(2));
599      }
600    }//There are more y than x
601  }//The degree is not zero
602  else
603  {//The degree is zero
604    hath = h;
605  }//The degree is zero
606  dbprint(p,dbprintWhitespace+" Done");
607  //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1)
608  list mons;
609  dbprint(p,dbprintWhitespace+" Putting the monomials in the A_0-part in a list.");
610  for(i = 1; i<=size(hath);i++)
611  {//Putting the monomials in a list
612    mons = mons+list(hath[i]);
613  }//Putting the monomials in a list
614  dbprint(p,dbprintWhitespace+" Done");
615  dbprint(p,dbprintWhitespace+" Mapping this monomials to K[theta]");
616  ring tempRing = 0,(x,y,theta),dp;
617  setring tempRing;
618  map thetamap = r,x,y;
619  list mons = thetamap(mons);
620  poly entry;
621  for (i = 1; i<=size(mons);i++)
622  {//transforming the monomials as monomials in theta
623    entry = leadcoef(mons[i]);
624    for (j = 0; j<leadexp(mons[i])[2];j++)
625    {
626      entry = entry * (theta-j);
627    }
628    mons[i] = entry;
629  }//transforming the monomials as monomials in theta
630  dbprint(p,dbprintWhitespace+" Done");
631  dbprint(p,dbprintWhitespace+" Factorize the A_0-Part in K[theta]");
632  list azeroresult = factorize(sum(mons));
633  dbprint(p,dbprintWhitespace+" Successful");
634  list azeroresult_return_form;
635  for (i = 1; i<=size(azeroresult[1]);i++)
636  {//rewrite the result of the commutative factorization
637    for (j = 1; j <= azeroresult[2][i];j++)
638    {
639      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
640    }
641  }//rewrite the result of the commutative factorization
642  dbprint(p,dbprintWhitespace+" Mapping back to A_0.");
643  setring(r);
644  map finalmap = tempRing,var(1),var(2),var(1)*var(2);
645  list tempresult = finalmap(azeroresult_return_form);
646  dbprint(p,dbprintWhitespace+"Successful.");
647  for (i = 1; i<=size(tempresult);i++)
648  {//factorizations of theta resp. theta +1
649    if(tempresult[i]==var(1)*var(2))
650    {
651      tempresult = insert(tempresult,var(1),i-1);
652      i++;
653      tempresult[i]=var(2);
654    }
655    if(tempresult[i]==var(2)*var(1))
656    {
657      tempresult = insert(tempresult,var(2),i-1);
658      i++;
659      tempresult[i]=var(1);
660    }
661  }//factorizations of theta resp. theta +1
662  result = tempresult+result;
663  return(result);
664}//proc homogfacFirstWeyl
665/* example */
666/* { */
667/*      "EXAMPLE:";echo=2; */
668/*      ring R = 0,(x,y),Ws(-1,1); */
669/*      def r = nc_algebra(1,1); */
670/*      setring(r); */
671/*      poly h = (x^2*y^2+1)*(x^4); */
672/*      homogfacFirstWeyl(h); */
673/* } */
674
675//==================================================
676//Computes all possible homogeneous factorizations
677static proc homogfacFirstWeyl_all(poly h)
678"USAGE: homogfacFirstWeyl_all(h); h is a homogeneous polynomial in the first Weyl algebra
679 with respect to the weight vector [-1,1]
680RETURN: list
681PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect
682  to the weight vector [-1,1] in the first Weyl algebra
683THEORY: @code{homogfacFirstWeyl} returns a list with all factorization of the given,
684 homogeneous polynomial. It uses the output of homogfacFirstWeyl and permutes
685 its entries with respect to the commutation rule. Furthermore, if a
686 factor of degree zero is irreducible in K[  heta], but reducible in
687 the first Weyl algebra, the permutations of this element with the other
688 entries will also be computed.
689SEE ALSO: homogfacFirstWeyl
690"{//proc HomogfacFirstWeylAll
691  int p=printlevel-voice+2;//for dbprint
692  intvec iv11= intvec(1,1);
693  if (deg(h,iv11) <= 0 )
694  {//h is a constant
695    dbprint(p,"Given polynomial was not homogeneous");
696    return(list(list(h)));
697  }//h is a constant
698  def r = basering;
699  list one_hom_fac; //stands for one homogeneous factorization
700  int i; int j; int k;
701  string dbprintWhitespace = "";
702  for (i = 1; i<=voice;i++)
703  {dbprintWhitespace = dbprintWhitespace + " ";}
704  intvec ivm11 = intvec(-1,1);
705  dbprint(p,dbprintWhitespace +" Calculate one homogeneous factorization using homogfacFirstWeyl");
706  //Compute again a homogeneous factorization
707  one_hom_fac = homogfacFirstWeyl(h);
708  dbprint(p,dbprintWhitespace +"Successful");
709  if (size(one_hom_fac) == 0)
710  {//there is no homogeneous factorization or the polynomial was not homogeneous
711    return(list());
712  }//there is no homogeneous factorization or the polynomial was not homogeneous
713  //divide list in A0-Part and a list of x's resp. y's
714  list list_not_azero = list();
715  list list_azero;
716  list k_factor;
717  int is_list_not_azero_empty = 1;
718  int is_list_azero_empty = 1;
719  k_factor = list(one_hom_fac[1]);
720  if (absValue(deg(h,ivm11))<size(one_hom_fac)-1)
721  {//There is a nontrivial A_0-part
722    list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,ivm11)))];
723    is_list_azero_empty = 0;
724  }//There is a nontrivial A_0 part
725  dbprint(p,dbprintWhitespace +" Combine x,y to xy in the factorization again.");
726  for (i = 1; i<=size(list_azero)-1;i++)
727  {//in homogfacFirstWeyl, we factorized theta, and this will be made undone
728    if (list_azero[i] == var(1))
729    {
730      if (list_azero[i+1]==var(2))
731      {
732        list_azero[i] = var(1)*var(2);
733        list_azero = delete(list_azero,i+1);
734      }
735    }
736    if (list_azero[i] == var(2))
737    {
738      if (list_azero[i+1]==var(1))
739      {
740        list_azero[i] = var(2)*var(1);
741        list_azero = delete(list_azero,i+1);
742      }
743    }
744  }//in homogfacFirstWeyl, we factorized theta, and this will be made undone
745  dbprint(p,dbprintWhitespace +" Done");
746  if(deg(h,ivm11)!=0)
747  {//list_not_azero is not empty
748    list_not_azero =
749      one_hom_fac[(size(one_hom_fac)-absValue(deg(h,ivm11))+1)..size(one_hom_fac)];
750    is_list_not_azero_empty = 0;
751  }//list_not_azero is not empty
752  //Map list_azero in K[theta]
753  dbprint(p,dbprintWhitespace +" Map list_azero to K[theta]");
754  ring tempRing = 0,(x,y,theta), dp;
755  setring(tempRing);
756  poly entry;
757  map thetamap = r,x,y;
758  if(!is_list_not_azero_empty)
759  {//Mapping in Singular is only possible, if the list before
760    //contained at least one element of the other ring
761    list list_not_azero = thetamap(list_not_azero);
762  }//Mapping in Singular is only possible, if the list before
763  //contained at least one element of the other ring
764  if(!is_list_azero_empty)
765  {//Mapping in Singular is only possible, if the list before
766    //contained at least one element of the other ring
767    list list_azero= thetamap(list_azero);
768  }//Mapping in Singular is only possible, if the list before
769  //contained at least one element of the other ring
770  list k_factor = thetamap(k_factor);
771  list tempmons;
772  dbprint(p,dbprintWhitespace +" Done");
773  for(i = 1; i<=size(list_azero);i++)
774  {//rewrite the polynomials in A1 as polynomials in K[theta]
775    tempmons = list();
776    for (j = 1; j<=size(list_azero[i]);j++)
777    {
778      tempmons = tempmons + list(list_azero[i][j]);
779    }
780    for (j = 1 ; j<=size(tempmons);j++)
781    {
782      entry = leadcoef(tempmons[j]);
783      for (k = 0; k < leadexp(tempmons[j])[2];k++)
784      {
785        entry = entry*(theta-k);
786      }
787      tempmons[j] = entry;
788    }
789    list_azero[i] = sum(tempmons);
790  }//rewrite the polynomials in A1 as polynomials in K[theta]
791  //Compute all permutations of the A0-part
792  dbprint(p,dbprintWhitespace +" Compute all permutations of the A_0-part with the first resp.
793the snd. variable");
794  list result;
795  int shift_sign;
796  int shift;
797  poly shiftvar;
798  if (size(list_not_azero)!=0)
799  {//Compute all possibilities to permute the x's resp. the y's in the list
800    if (list_not_azero[1] == x)
801    {//h had a negative weighted degree
802      shift_sign = 1;
803      shiftvar = x;
804    }//h had a negative weighted degree
805    else
806    {//h had a positive weighted degree
807      shift_sign = -1;
808      shiftvar = y;
809    }//h had a positive weighted degree
810    result = permpp(list_azero + list_not_azero);
811    for (i = 1; i<= size(result); i++)
812    {//adjust the a_0-parts
813      shift = 0;
814      for (j=1; j<=size(result[i]);j++)
815      {
816        if (result[i][j]==shiftvar)
817        {
818          shift = shift + shift_sign;
819        }
820        else
821        {
822          result[i][j] = subst(result[i][j],theta,theta + shift);
823        }
824      }
825    }//adjust the a_0-parts
826  }//Compute all possibilities to permute the x's resp. the y's in the list
827  else
828  {//The result is just all the permutations of the a_0-part
829    result = permpp(list_azero);
830  }//The result is just all the permutations of the a_0 part
831  if (size(result)==0)
832  {
833    return(result);
834  }
835  dbprint(p,dbprintWhitespace +" Done");
836  dbprint(p,dbprintWhitespace +" Searching for theta resp. theta+1 in the list and fact. them");
837  //Now we are going deeper and search for theta resp. theta + 1, substitute
838  //them by xy resp. yx and go on permuting
839  int found_theta;
840  int thetapos;
841  list leftpart;
842  list rightpart;
843  list lparts;
844  list rparts;
845  list tempadd;
846  for (i = 1; i<=size(result) ; i++)
847  {//checking every entry of result for theta or theta +1
848    found_theta = 0;
849    for(j=1;j<=size(result[i]);j++)
850    {
851      if (result[i][j]==theta)
852      {//the jth entry is theta and can be written as x*y
853        thetapos = j;
854        result[i]= insert(result[i],x,j-1);
855        j++;
856        result[i][j] = y;
857        found_theta = 1;
858        break;
859      }//the jth entry is theta and can be written as x*y
860      if(result[i][j] == theta +1)
861      {
862        thetapos = j;
863        result[i] = insert(result[i],y,j-1);
864        j++;
865        result[i][j] = x;
866        found_theta = 1;
867        break;
868      }
869    }
870    if (found_theta)
871    {//One entry was theta resp. theta +1
872      leftpart = result[i];
873      leftpart = leftpart[1..thetapos];
874      rightpart = result[i];
875      rightpart = rightpart[(thetapos+1)..size(rightpart)];
876      lparts = list(leftpart);
877      rparts = list(rightpart);
878      //first deal with the left part
879      if (leftpart[thetapos] == x)
880      {
881        shift_sign = 1;
882        shiftvar = x;
883      }
884      else
885      {
886        shift_sign = -1;
887        shiftvar = y;
888      }
889      for (j = size(leftpart); j>1;j--)
890      {//drip x resp. y
891        if (leftpart[j-1]==shiftvar)
892        {//commutative
893          j--;
894          continue;
895        }//commutative
896        if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
897        {//stop here
898          break;
899        }//stop here
900        //Here, we can only have a a0- part
901        leftpart[j] = subst(leftpart[j-1],theta, theta + shift_sign);
902        leftpart[j-1] = shiftvar;
903        lparts = lparts + list(leftpart);
904      }//drip x resp. y
905      //and now deal with the right part
906      if (rightpart[1] == x)
907      {
908        shift_sign = 1;
909        shiftvar = x;
910      }
911      else
912      {
913        shift_sign = -1;
914        shiftvar = y;
915      }
916      for (j = 1 ; j < size(rightpart); j++)
917      {
918        if (rightpart[j+1] == shiftvar)
919        {
920          j++;
921          continue;
922        }
923        if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
924        {
925          break;
926        }
927        rightpart[j] = subst(rightpart[j+1], theta, theta - shift_sign);
928        rightpart[j+1] = shiftvar;
929        rparts = rparts + list(rightpart);
930      }
931      //And now, we put all possibilities together
932      tempadd = list();
933      for (j = 1; j<=size(lparts); j++)
934      {
935        for (k = 1; k<=size(rparts);k++)
936        {
937          tempadd = tempadd + list(lparts[j]+rparts[k]);
938        }
939      }
940      tempadd = delete(tempadd,1); // The first entry is already in the list
941      result = result + tempadd;
942      continue; //We can may be not be done already with the ith entry
943    }//One entry was theta resp. theta +1
944  }//checking every entry of result for theta or theta +1
945  dbprint(p,dbprintWhitespace +" Done");
946  //map back to the basering
947  dbprint(p,dbprintWhitespace +" Mapping back everything to the basering");
948  setring(r);
949  map finalmap = tempRing, var(1), var(2),var(1)*var(2);
950  list result = finalmap(result);
951  for (i=1; i<=size(result);i++)
952  {//adding the K factor
953    result[i] = k_factor + result[i];
954  }//adding the k-factor
955  dbprint(p,dbprintWhitespace +" Done");
956  dbprint(p,dbprintWhitespace +" Delete double entries in the list.");
957  result = delete_dublicates_noteval(result);
958  dbprint(p,dbprintWhitespace +" Done");
959  return(result);
960}//proc HomogfacFirstWeylAll
961/* example */
962/* { */
963/*      "EXAMPLE:";echo=2; */
964/*      ring R = 0,(x,y),Ws(-1,1); */
965/*      def r = nc_algebra(1,1); */
966/*      setring(r); */
967/*      poly h = (x^2*y^2+1)*(x^4); */
968/*      homogfacFirstWeyl_all(h); */
969/* } */
970
971//==================================================*
972//Computes all permutations of a given list
973static proc perm(list l)
974"
975DEPRECATED
976"
977{//proc perm
978  int i; int j;
979  list tempresult;
980  list result;
981  if (size(l)==0)
982  {
983    return(list());
984  }
985  if (size(l)==1)
986  {
987    return(list(l));
988  }
989  for (i = 1; i<=size(l); i++ )
990  {
991    tempresult = perm(delete(l,i));
992    for (j = 1; j<=size(tempresult);j++)
993    {
994      tempresult[j] = list(l[i])+tempresult[j];
995    }
996    result = result+tempresult;
997  }
998  return(result);
999}//proc perm
1000
1001//==================================================
1002//computes all permutations of a given list by
1003//ignoring equal entries (faster than perm)
1004static proc permpp(list l)
1005"
1006INPUT: A list with entries of a type, where the ==-operator is defined
1007OUTPUT: A list with all permutations of this given list.
1008"
1009{//proc permpp
1010  int i; int j;
1011  list tempresult;
1012  list l_without_double;
1013  list l_without_double_pos;
1014  int double_entry;
1015  list result;
1016  if (size(l)==0)
1017  {
1018    return(list());
1019  }
1020  if (size(l)==1)
1021  {
1022    return(list(l));
1023  }
1024  for (i = 1; i<=size(l);i++)
1025  {//Filling the list with unique entries
1026    double_entry = 0;
1027    for (j = 1; j<=size(l_without_double);j++)
1028    {
1029      if (l_without_double[j] == l[i])
1030      {
1031        double_entry = 1;
1032        break;
1033      }
1034    }
1035    if (!double_entry)
1036    {
1037      l_without_double = l_without_double + list(l[i]);
1038      l_without_double_pos = l_without_double_pos + list(i);
1039    }
1040  }//Filling the list with unique entries
1041  for (i = 1; i<=size(l_without_double); i++ )
1042  {
1043    tempresult = permpp(delete(l,l_without_double_pos[i]));
1044    for (j = 1; j<=size(tempresult);j++)
1045    {
1046      tempresult[j] = list(l_without_double[i])+tempresult[j];
1047    }
1048    result = result+tempresult;
1049  }
1050  return(result);
1051}//proc permpp
1052
1053//==================================================
1054//factorization of the first Weyl Algebra
1055
1056//The following procedure just serves the purpose to
1057//transform the input into an appropriate input for
1058//the procedure sfacwa, where the ring must contain the
1059//variables in a certain order.
1060proc facFirstWeyl(poly h)
1061"USAGE: facFirstWeyl(h); h a polynomial in the first Weyl algebra
1062RETURN: list
1063PURPOSE: compute all factorizations of a polynomial in the first Weyl algebra
1064THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle
1065ASSUME: basering is the first Weyl algebra
1066NOTE: Every entry of the output list is a list with factors for one possible factorization.
1067The first factor is always a constant (1, if no nontrivial constant could be excluded).
1068EXAMPLE: example facFirstWeyl; shows examples
1069SEE ALSO: facSubWeyl, testNCfac, facFirstShift
1070"{//proc facFirstWeyl
1071  //Definition of printlevel variable
1072  int p = printlevel-voice+2;
1073  int i;
1074  string dbprintWhitespace = "";
1075  for (i = 1; i<=voice;i++)
1076  {dbprintWhitespace = dbprintWhitespace + " ";}
1077  dbprint(p,dbprintWhitespace +" Checking if the given algebra is a Weyl algebra");
1078  //Redefine the ring in my standard form
1079  if (!isWeyl())
1080  {//Our basering is not the Weyl algebra
1081    ERROR("Ring was not the first Weyl algebra");
1082    return(list());
1083  }//Our basering is not the Weyl algebra
1084  dbprint(p,dbprintWhitespace +" Successful");
1085  dbprint(p,dbprintWhitespace +" Checking, if the given ring is the first Weyl algebra");
1086  if(nvars(basering)!=2)
1087  {//Our basering is the Weyl algebra, but not the first
1088    ERROR("Ring is not the first Weyl algebra");
1089    return(list());
1090  }//Our basering is the Weyl algebra, but not the first
1091
1092  //A last check before we start the real business: Is h already given as a polynomial just
1093  //in one variable?
1094  if (deg(h,intvec(1,0))== 0 or deg(h,intvec(0,1)) == 0)
1095  {//h is in K[x] or in K[d]
1096    if (deg(h,intvec(1,0))== 0 and deg(h,intvec(0,1)) == 0)
1097    {//We just have a constant
1098      return(list(list(h)));
1099    }//We just have a constant
1100    dbprint(p,dbprintWhitespace+"Polynomial was given in one variable.
1101Performing commutative factorization.");
1102    int theCommVar;
1103    if (deg(h,intvec(1,0)) == 0)
1104    {//The second variable is the variable to factorize
1105      theCommVar = 2;
1106    }//The second variable is the variable to factorize
1107    else{theCommVar = 1;}
1108    def r = basering;
1109    ring tempRing = 0,(var(theCommVar)),dp;
1110    if (theCommVar == 1){map mapToCommutative = r,var(1),1;}
1111    else {map mapToCommutative = r,1,var(1);}
1112    poly h = mapToCommutative(h);
1113    list tempResult = factorize(h);
1114    list result = list(list());
1115    int j;
1116    for (i = 1; i<=size(tempResult[1]); i++)
1117    {
1118      for (j = 1; j<=tempResult[2][i]; j++)
1119      {
1120        result[1] = result[1] + list(tempResult[1][i]);
1121      }
1122    }
1123    //mapping back
1124    setring(r);
1125    map mapBackFromCommutative = tempRing,var(theCommVar);
1126    def result = mapBackFromCommutative(result);
1127    dbprint(p,dbprintWhitespace+"result:");
1128    dbprint(p,result);
1129    dbprint(p,dbprintWhitespace+"Computing all permutations of this factorization");
1130    poly constantFactor = result[1][1];
1131    result[1] = delete(result[1],1);//Deleting the constant factor
1132    result=permpp(result[1]);
1133    for (i = 1; i<=size(result);i++)
1134    {//Insert constant factor
1135      result[i] = insert(result[i],constantFactor);
1136    }//Insert constant factor
1137    dbprint(p,dbprintWhitespace+"Done.");
1138    return(result);
1139  }//h is in K[x] or in K[d]
1140
1141
1142  dbprint(p,dbprintWhitespace +" Successful");
1143  list result = list();
1144  int j; int k; int l; //counter
1145  if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why
1146  {
1147    dbprint(p,dbprintWhitespace +" positions of the variables have to be switched");
1148    def r = basering;
1149    ring tempRing = ringlist(r)[1][1],(x,y),Ws(-1,1); // very strange:
1150    // setting Wp(-1,1) leads to SegFault; to clarify why!!!
1151    def NTR = nc_algebra(1,1);
1152    setring NTR ;
1153    map transf = r, var(2), var(1);
1154    dbprint(p,dbprintWhitespace +" Successful");
1155    list resulttemp = sfacwa(h);
1156    setring(r);
1157    map transfback = NTR, var(2),var(1);
1158    result = transfback(resulttemp);
1159  }
1160  else
1161  {
1162    dbprint(p, dbprintWhitespace +" factorization of the polynomial with the routine sfacwa");
1163    result = sfacwa(h);
1164    dbprint(p,dbprintWhitespace +" Done");
1165  }
1166  if (homogwithorder(h,intvec(-1,1)))
1167  {
1168    dbprint(p, dbprintWhitespace + " Polynomial was homogeneous, therefore we have
1169already a complete factorization and do not have to go through the factors recursively.");
1170    return(result);
1171  }
1172  dbprint(p,dbprintWhitespace + "We have the following intermediate list of inhomogeneous
1173factorizations:");
1174  dbprint(p,result);
1175  dbprint(p,dbprintWhitespace +" recursively check factors for irreducibility");
1176  list recursivetemp;
1177  int changedSomething;
1178  for(i = 1; i<=size(result);i++)
1179  {//recursively factorize factors
1180    if(size(result[i])>2)
1181    {//Nontrivial factorization
1182      for (j=2;j<=size(result[i]);j++)
1183      {//Factorize every factor
1184        recursivetemp = facFirstWeyl(result[i][j]);
1185        //if(size(recursivetemp)>1)
1186        //{//we have a nontrivial factorization
1187        changedSomething = 0;
1188        for(k=1; k<=size(recursivetemp);k++)
1189        {//insert factorized factors
1190          if(size(recursivetemp[k])>2)
1191          {//nontrivial
1192            changedSomething = 1;
1193            result = insert(result,result[i],i);
1194            for(l = size(recursivetemp[k]);l>=2;l--)
1195            {
1196              result[i+1] = insert(result[i+1],recursivetemp[k][l],j);
1197            }
1198            result[i+1] = delete(result[i+1],j);
1199          }//nontrivial
1200        }//insert factorized factors
1201        if (changedSomething)
1202        {
1203          result = delete(result,i);
1204        }
1205        //}//we have a nontrivial factorization
1206      }//Factorize every factor
1207    }//Nontrivial factorization
1208  }//recursively factorize factors
1209  dbprint(p,dbprintWhitespace +" Done");
1210  if (size(result)==0)
1211  {//only the trivial factorization could be found
1212    result = list(list(1,h));
1213  }//only the trivial factorization could be found
1214  list resultWithInterchanges;
1215  dbprint(p,dbprintWhitespace+ "And the result without interchanges with homogeneous factors is:");
1216  dbprint(p,result);
1217  for (i = 1; i <= size(result) ; i++)
1218  {//applying the interchanges to result
1219    resultWithInterchanges = resultWithInterchanges +
1220                             checkForHomogInhomogInterchangability(result[i],2,size(result[i]));
1221  }//applying the interchanges to result
1222  dbprint(p,dbprintWhitespace + "With interchanges, the result is:");
1223  dbprint(p,resultWithInterchanges);
1224  //now, refine the possible redundant list
1225  return( delete_dublicates_noteval(resultWithInterchanges) );
1226}//proc facFirstWeyl
1227example
1228{
1229  "EXAMPLE:";echo=2;
1230  ring R = 0,(x,y),dp;
1231  def r = nc_algebra(1,1);
1232  setring(r);
1233  poly h = (x^2*y^2+x)*(x+1);
1234  facFirstWeyl(h);
1235}
1236
1237//////////////////////////////////////////////////
1238/////BRANDNEW!!!!////////////////////
1239//////////////////////////////////////////////////
1240
1241static proc checkForHomogInhomogInterchangability(list factors, posLeft, posRight)
1242"
1243INPUT:  A list consisting of factors of a certain polynomial in the first Weyl
1244        algebra, factors, and a position from the left and the right, where the last swap was done.
1245OUTPUT: A list containing lists consisting of factors of a certain polynomial in the first Weyl
1246        algebra.
1247The purpose of this function is to check whether we can interchange certain inhomogeneous factors
1248with homogeneous ones. If it is possible, this function returns a list of lists
1249of possible interchanges.
1250
1251The idea came because of an example, where we need an extra swap in the end, otherwise we would
1252not capture all factorizations. The example was
1253h = x4d7+11x3d6+x2d7+x2d6+x3d4+29x2d5+xd6+8xd5+d6+5x2d3+14xd4+13d4+5xd2+d3+d;
1254
1255ASSUMPTIONS:
1256
1257- All factors are irreducible
1258"
1259{//checkForHomogInhomogInterchangability
1260  int p = printlevel-voice+2;
1261  string dbprintWhitespace = "";
1262  int i; int j; int k;
1263  for (i = 1; i<=voice;i++)
1264  {dbprintWhitespace = dbprintWhitespace + " ";}
1265  if (size(factors) <= 2 || posLeft >= posRight - 1)
1266  {//easiest case: There is nothing to swap
1267    return (list(factors));
1268  }//easiest case: There is nothing to swap
1269  list result = list(factors);
1270  list tempResultEntries;
1271  list tempSwaps;
1272  list tempSwapsTempEntry;
1273  list attemptToSwap;
1274  intvec ivm11 = intvec(-1,1);
1275  dbprint(p, dbprintWhitespace+"We try to swap elements in the following list:");
1276  dbprint(p, factors);
1277  for (i = posLeft; i < posRight; i++)
1278  {//checking within the window posLeft <--> posRight, if there are interchanges possible
1279    if (homogwithorder(factors[i],ivm11) && !homogwithorder(factors[i+1],ivm11))
1280    {//position i is homogeneous, position i+1 is not ==> trying to swap
1281      attemptToSwap = extractHomogeneousDivisorsRight(factors[i]*factors[i+1]);
1282      if (size(attemptToSwap[1])>1)
1283      {//Bingo, we were able to swap this one element
1284        dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1));
1285        dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:");
1286        dbprint(p,attemptToSwap);
1287        tempSwapsTempEntry = list();
1288        for (j = size(factors); j >=1; j--)
1289        {//creating a new entry for the resulting list, replacing the swap in factors
1290          if (j==i+1)
1291          {
1292            for (k = size(attemptToSwap[1]); k >=1 ; k--)
1293            {
1294              tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][k]);
1295            }
1296            j--; //Because we changed entry i+1 and i
1297          }
1298          else
1299          {
1300            tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]);
1301          }
1302        }//creating a new entry for the resulting list, replacing the swap in factors
1303        tempSwaps = insert(tempSwaps,list(list(i+1,posRight),tempSwapsTempEntry));
1304      }//Bingo, we were able to swap this one element
1305    }//position i is homogeneous, position i+1 is not ==> trying to swap
1306    else
1307    {
1308      if(!homogwithorder(factors[i],ivm11) && homogwithorder(factors[i+1],ivm11))
1309      {//position i+1 is homogeneous, position i is not ==> trying to swap
1310        attemptToSwap = extractHomogeneousDivisorsLeft(factors[i]*factors[i+1]);
1311        if (size(attemptToSwap[1])>1)
1312        {//Bingo, we were able to swap this one element
1313          dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1));
1314          dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:");
1315          dbprint(p,attemptToSwap);
1316          tempSwapsTempEntry = list();
1317          for (j = size(factors); j >=1; j--)
1318          {//creating a new entry for the resulting list, replacing the swap in factors
1319            if (j==i+1)
1320            {
1321              for (k = size(attemptToSwap[1]); k >=1 ; k--)
1322              {
1323                tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][k]);
1324              }
1325              j--; //Because we changed entry i+1 and i
1326            }
1327            else
1328            {
1329              tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]);
1330            }
1331          }//creating a new entry for the resulting list, replacing the swap in factors
1332          tempSwaps = insert(tempSwaps,list(list(posLeft,i),tempSwapsTempEntry));
1333        }//Bingo, we were able to swap this one element
1334      }//position i+1 is homogeneous, position i is not ==> trying to swap
1335    }
1336  }//checking within the window posLeft <--> posRight, if there are interchanges possible
1337  //Now we will recursively call the function for all swapped entries.
1338  dbprint(p,dbprintWhitespace+ "Our list of different factorizations is now:");
1339  dbprint(p,tempSwaps);
1340  for (i = 1; i<=size(tempSwaps);i++)
1341  {//recursive call to all formerly attempted swaps.
1342    tempResultEntries=checkForHomogInhomogInterchangability(tempSwaps[i][2],
1343                                                            tempSwaps[i][1][1],tempSwaps[i][1][2]);
1344    result = result + tempResultEntries;
1345  }//recursive call to all formerly attempted swaps.
1346  result = delete_dublicates_noteval(result);
1347  return(result);
1348}//checkForHomogInhomogInterchangability
1349
1350static proc sfacwa(poly h)
1351"INPUT: A polynomial h in the first Weyl algebra
1352OUTPUT: A list of factorizations, where the factors might still be reducible.
1353ASSUMPTIONS:
1354- Our basering is the first Weyl algebra; the x is the first variable,
1355  the differential operator the second.
1356"
1357{//proc sfacwa
1358  int i; int j; int k;
1359  int p = printlevel-voice+2;
1360  string dbprintWhitespace = "";
1361  number commonCoefficient = content(h);
1362  for (i = 1; i<=voice;i++)
1363  {dbprintWhitespace = dbprintWhitespace + " ";}
1364  dbprint(p,dbprintWhitespace + " Extracting homogeneous left and right factors");
1365  if(homogwithorder(h,intvec(-1,1)))
1366  {//we are already dealing with a -1,1 homogeneous poly
1367    dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous. Start homog.
1368fac. and ret. its result");
1369    return(homogfacFirstWeyl_all(h));
1370  }//we are already dealing with a -1,1 homogeneous poly
1371  list resulttemp = extractHomogeneousDivisors(h/commonCoefficient);
1372  //resulttemp = resulttemp + list(list(h/commonCoefficient));
1373  list inhomogeneousFactorsToFactorize;
1374  int isAlreadyInInhomogList;
1375  dbprint(p,dbprintWhitespace +" Done");
1376  dbprint(p,dbprintWhitespace +" Making Set of inhomogeneous polynomials we have to factorize.");
1377  for (i = 1; i<=size(resulttemp); i++)
1378  {//Going through all different kinds of factorizations where we extracted homogeneous factors
1379    for (j = 1;j<=size(resulttemp[i]);j++)
1380    {//searching for the inhomogeneous factor
1381      if (!homogwithorder(resulttemp[i][j],intvec(-1,1)))
1382      {//We have found our candidate
1383        isAlreadyInInhomogList = 0;
1384        for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++)
1385        {//Checking if our candidate is already in our tofactorize-list
1386          if (inhomogeneousFactorsToFactorize[k]==resulttemp[i][j])
1387          {//The candidate was already in the list
1388            isAlreadyInInhomogList = 1;
1389            break;
1390          }//The candidate was already in the list
1391        }//Checking if our candidate is already in our tofactorize-list
1392        if (!isAlreadyInInhomogList)
1393        {
1394          inhomogeneousFactorsToFactorize=inhomogeneousFactorsToFactorize + list(resulttemp[i][j]);
1395        }
1396      }//We have found our candidate
1397    }//searching for the inhomogeneous factor
1398  }//Going through all different kinds of factorizations where we extracted homogeneous factors
1399  dbprint(p,dbprintWhitespace +" Done");
1400  dbprint(p,dbprintWhitespace + "The set is:");
1401  dbprint(p,inhomogeneousFactorsToFactorize);
1402  dbprint(p,dbprintWhitespace+ "Factorizing the different occuring inhomogeneous factors");
1403  for (i = 1; i<= size(inhomogeneousFactorsToFactorize); i++)
1404  {//Factorizing all kinds of inhomogeneous factors
1405    inhomogeneousFactorsToFactorize[i] = sfacwa2(inhomogeneousFactorsToFactorize[i]);
1406    for (j = 1; j<=size(inhomogeneousFactorsToFactorize[i]);j++)
1407    {//Deleting the leading coefficient since we don't need him
1408      if (deg(inhomogeneousFactorsToFactorize[i][j][1],intvec(1,1))==0)
1409      {
1410        inhomogeneousFactorsToFactorize[i][j] = delete(inhomogeneousFactorsToFactorize[i][j],1);
1411      }
1412    }//Deleting the leading coefficient since we don't need him
1413  }//Factorizing all kinds of inhomogeneous factors
1414  dbprint(p,dbprintWhitespace +" Done");
1415  dbprint(p,dbprintWhitespace +" Putting the factorizations in the lists");
1416  list result;
1417  int posInhomogPoly;
1418  int posInhomogFac;
1419  for (i = 1; i<=size(resulttemp); i++)
1420  {//going through all by now calculated factorizations
1421    for (j = 1;j<=size(resulttemp[i]); j++)
1422    {//Finding the inhomogeneous factor
1423      if (!homogwithorder(resulttemp[i][j],intvec(-1,1)))
1424      {//Found it
1425        posInhomogPoly = j;
1426        break;
1427      }//Found it
1428    }//Finding the inhomogeneous factor
1429    for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++)
1430    {//Finding the matching inhomogeneous factorization we already determined
1431      if(product(inhomogeneousFactorsToFactorize[k][1]) == resulttemp[i][j])
1432      {//found it
1433        posInhomogFac = k;
1434        break;
1435      }//Found it
1436    }//Finding the matching inhomogeneous factorization we already determined
1437    for (j = 1; j <= size(inhomogeneousFactorsToFactorize[posInhomogFac]); j++)
1438    {
1439      result = insert(result, resulttemp[i]);
1440      result[1] = delete(result[1],posInhomogPoly);
1441      for (k =size(inhomogeneousFactorsToFactorize[posInhomogFac][j]);k>=1; k--)
1442      {//Inserting factorizations
1443        result[1] = insert(result[1],inhomogeneousFactorsToFactorize[posInhomogFac][j][k],
1444                           posInhomogPoly-1);
1445      }//Inserting factorizations
1446      dbprint(p,dbprintWhitespace + "Added a factorization to result, namely:");
1447      dbprint(p, result[1]);
1448    }
1449  }//going through all by now calculated factorizations
1450  dbprint(p,dbprintWhitespace +" Done");
1451  result = delete_dublicates_noteval(result);
1452  for (i = 1; i<=size(result);i++)
1453  {//Putting the content everywhere
1454    result[i] = insert(result[i],commonCoefficient);
1455  }//Putting the content everywhere
1456  return(result);
1457}//proc sfacwa
1458
1459static proc sfacwa2(poly h)
1460"
1461Subprocedure of sfacwa
1462Assumptions:
1463- h is not in K[x] or in K[d], or even in K. These cases are caught by the input
1464- The coefficients are integer values and the gcd of the coefficients is 1
1465"
1466{//proc sfacwa2
1467  int p=printlevel-voice+2; // for dbprint
1468  int i;
1469  string dbprintWhitespace = "";
1470  for (i = 1; i<=voice;i++)
1471  {dbprintWhitespace = dbprintWhitespace + " ";}
1472  intvec ivm11 = intvec(-1,1);
1473  intvec iv11 = intvec(1,1);
1474  intvec iv10 = intvec(1,0);
1475  intvec iv01 = intvec(0,1);
1476  intvec iv1m1 = intvec(1,-1);
1477  poly p_max; poly p_min; poly q_max; poly q_min;
1478  map invo = basering,-var(1),var(2);
1479  list calculatedRightFactors;
1480  if(homogwithorder(h,ivm11))
1481  {//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here
1482    dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous.
1483Start homog. fac. and ret. its result");
1484    return(homogfacFirstWeyl_all(h));
1485  }//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here
1486  list result = list();
1487  int j; int k; int l;
1488  dbprint(p,dbprintWhitespace+" Computing the degree-limits of the factorization");
1489  //end finding the limits
1490  dbprint(p,dbprintWhitespace+" Computing the maximal and the minimal
1491homogeneous part of the given polynomial");
1492  list M = computeCombinationsMinMaxHomog(h);
1493  dbprint(p,dbprintWhitespace+" Done.");
1494  dbprint(p,dbprintWhitespace+" Filtering invalid combinations in M.");
1495  for (i = 1 ; i<= size(M); i++)
1496  {//filter valid combinations
1497    if (product(M[i]) == h)
1498    {//We have one factorization
1499      result = result + divides(M[i][1],h,invo,1);
1500      dbprint(p,dbprintWhitespace+"Result list updated:");
1501      dbprint(p,dbprintWhitespace+string(result));
1502      M = delete(M,i);
1503      continue;
1504    }//We have one factorization
1505  }//filter valid combinations
1506  dbprint(p,dbprintWhitespace+"Done.");
1507  dbprint(p,dbprintWhitespace+"The size of M is "+string(size(M)));
1508  for (i = 1; i<=size(M); i++)
1509  {//Iterate over all first combinations (p_max + p_min)(q_max + q_min)
1510    dbprint(p,dbprintWhitespace+" Combination No. "+string(i)+" in M:" );
1511    p_max = jet(M[i][1],deg(M[i][1],ivm11),ivm11)-jet(M[i][1],deg(M[i][1],ivm11)-1,ivm11);
1512    p_min = jet(M[i][1],deg(M[i][1],iv1m1),iv1m1)-jet(M[i][1],deg(M[i][1],iv1m1)-1,iv1m1);
1513    q_max = jet(M[i][2],deg(M[i][2],ivm11),ivm11)-jet(M[i][2],deg(M[i][2],ivm11)-1,ivm11);
1514    q_min = jet(M[i][2],deg(M[i][2],iv1m1),iv1m1)-jet(M[i][2],deg(M[i][2],iv1m1)-1,iv1m1);
1515    dbprint(p,dbprintWhitespace+" pmax = "+string(p_max));
1516    dbprint(p,dbprintWhitespace+" pmin = "+string(p_min));
1517    dbprint(p,dbprintWhitespace+" qmax = "+string(q_max));
1518    dbprint(p,dbprintWhitespace+" qmin = "+string(q_min));
1519    //Check, whether p_max + p_min or q_max and q_min are already left or right divisors.
1520    if (divides(p_min + p_max,h,invo))
1521    {
1522      dbprint(p,dbprintWhitespace+" Got one result.");
1523      result = result + divides(p_min + p_max,h,invo,1);
1524    }
1525    else
1526    {
1527      if (divides(q_min + q_max,h,invo))
1528      {
1529        dbprint(p,dbprintWhitespace+" Got one result.");
1530        result = result + divides(q_min + q_max, h , invo, 1);
1531      }
1532    }
1533    //Now the check, if deg(p_max) = deg(p_min)+1 (and the same with q_max and q_min)
1534
1535    if (deg(p_max, ivm11) == deg(p_min, ivm11) +1 or  deg(q_max, ivm11) == deg(q_min, ivm11) +1 )
1536    {//Therefore, p_max + p_min must be a left factor or we can dismiss the combination
1537      dbprint(p,dbprintWhitespace+" There are no homogeneous parts we can put between
1538pmax and pmin resp. qmax and qmin.");
1539      //TODO: Prove, that then also a valid right factor is not possible
1540      M = delete(M,i);
1541      continue;
1542    }//Therefore, p_max + p_min must be a left factor or we can dismiss the combination
1543
1544    //Done with the Check
1545
1546    //If we come here, there are still homogeneous parts to be added to p_max + p_min
1547    //AND to q_max and q_min in
1548    //order to obtain a real factor
1549    //We use the procedure determineRestOfHomogParts to find our q.
1550    dbprint(p,dbprintWhitespace+" Solving for the other homogeneous parts in q");
1551    calculatedRightFactors = determineRestOfHomogParts(p_max,p_min,q_max,q_min,h);
1552    dbprint(p,dbprintWhitespace+" Done with it. Found "+string(size(calculatedRightFactors))
1553            +" solutions.");
1554    for (j = 1; j<=size(calculatedRightFactors);j++)
1555    {//Check out whether we really have right factors of h in calculatedRightFactors
1556      if (divides(calculatedRightFactors[j],h,invo))
1557      {
1558        result = result + divides(calculatedRightFactors[j],h,invo,1);
1559      }
1560      else
1561      {
1562        dbprint(p,"Solution for max and min homog found, but not a divisor of h");
1563        //TODO: Proof, why this can happen.
1564      }
1565    }//Check out whether we really have right factors of h in calculatedRightFactors
1566  }//Iterate over all first combinations (p_max + p_min)(q_max + q_min)
1567
1568
1569  result = delete_dublicates_noteval(result);
1570  //print(M);
1571  if (size(result) == 0)
1572  {//no factorization found
1573    result = list(list(h));
1574  }//no factorization found
1575  return(result);
1576}//proc sfacwa2
1577
1578static proc determineRestOfHomogParts(poly pmax, poly pmin, poly qmax, poly qmin, poly h)
1579"INPUT: Polynomials p_max, p_min, q_max, q_min and h. The maximum homogeneous part h_max of h is
1580 given by p_max*pmin, the minimum homogeneous part h_min of h is given by p_min*q_min.
1581OUTPUT: A list of right factors q of h that have q_max and q_min as their maximum respectively
1582 minimum homogeneous part. Empty list, if those elements are not existent
1583ASSUMPTIONS:
1584 - deg(p_max,intvec(-1,1))>deg(p_min,intvec(-1,1)) +1
1585 - deg(q_max,intvec(-1,1))>deg(q_min,intvec(-1,1)) +1
1586 - p_max*q_max = h_max
1587 - p_min*q_min = h_min
1588 - The basering is the first Weyl algebra
1589"
1590{//proc determineRestOfHomogParts
1591  int p=printlevel-voice+2; // for dbprint
1592  string dbprintWhitespace = "";
1593  int i;
1594  for (i = 1; i<=voice;i++)
1595  {dbprintWhitespace = dbprintWhitespace + " ";}
1596  int kappa = Min(intvec(deg(h,intvec(1,0)), deg(h,intvec(0,1))));
1597  def R = basering;
1598  int n1 = deg(pmax,intvec(-1,1));
1599  int nk = -deg(pmin,intvec(1,-1));
1600  int m1 = deg(qmax,intvec(-1,1));
1601  int ml = -deg(qmin,intvec(1,-1));
1602  int j; int k;
1603  ideal mons;
1604
1605  dbprint(p,dbprintWhitespace+" Extracting zero homog. parts of pmax, qmax, pmin, qmin and h.");
1606  //Extracting the zero homogeneous part of the given polynomials
1607  ideal pandqZero = pmax,pmin,qmax,qmin;
1608  if (n1 > 0){pandqZero[1] = lift(var(2)^n1,pmax)[1,1];}
1609  else{if (n1 < 0){pandqZero[1] = lift(var(1)^(-n1),pmax)[1,1];}
1610    else{pandqZero[1] = pmax;}}
1611  if (nk > 0){pandqZero[2] = lift(var(2)^nk,pmin)[1,1];}
1612  else{if (nk < 0){pandqZero[2] = lift(var(1)^(-nk),pmin)[1,1];}
1613    else{pandqZero[2] = pmin;}}
1614  if (m1 > 0){pandqZero[3] = lift(var(2)^m1,qmax)[1,1];}
1615  else{if (m1 < 0){pandqZero[3] = lift(var(1)^(-m1),qmax)[1,1];}
1616    else{pandqZero[3] = qmax;}}
1617  if (ml > 0){pandqZero[4] = lift(var(2)^ml,qmin)[1,1];}
1618  else{if (ml < 0){pandqZero[4] = lift(var(1)^(-ml),qmin)[1,1];}
1619    else{pandqZero[4] = qmin;}}
1620  list hZeroinR = homogDistribution(h);
1621  for (i = 1; i<=size(hZeroinR);i++)
1622  {//Extracting the zero homogeneous parts of the homogeneous summands of h
1623    if (hZeroinR[i][1] > 0){hZeroinR[i][2] = lift(var(2)^hZeroinR[i][1],hZeroinR[i][2])[1,1];}
1624    if (hZeroinR[i][1] < 0){hZeroinR[i][2] = lift(var(1)^(-hZeroinR[i][1]),hZeroinR[i][2])[1,1];}
1625  }//Extracting the zero homogeneous parts of the homogeneous summands of h
1626  dbprint(p,dbprintWhitespace+" Done!");
1627  //Moving everything into the ring K[theta]
1628  dbprint(p,dbprintWhitespace+" Moving everything into the ring K[theta]");
1629  ring KTheta = 0,(x,d,theta),dp;
1630  map thetamap = R, x, d;
1631  poly entry;
1632  ideal mons;
1633  ideal pandqZero;
1634  list hZeroinKTheta;
1635  setring(R);
1636
1637  //Starting with p and q
1638  for (k=1; k<=4; k++)
1639  {//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys
1640    mons = ideal();
1641    for(i = 1; i<=size(pandqZero[k]);i++)
1642    {//Putting the monomials in a list
1643      mons[size(mons)+1] = pandqZero[k][i];
1644    }//Putting the monomials in a list
1645    setring(KTheta);
1646    mons = thetamap(mons);
1647    for (i = 1; i<=size(mons);i++)
1648    {//transforming the monomials as monomials in theta
1649      entry = leadcoef(mons[i]);
1650      for (j = 0; j<leadexp(mons[i])[2];j++)
1651      {
1652        entry = entry * (theta-j);
1653      }
1654      mons[i] = entry;
1655    }//transforming the monomials as monomials in theta
1656    pandqZero[size(pandqZero)+1] = sum(mons);
1657    setring(R);
1658  }//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys
1659
1660  //Now hZero
1661  for (k = size(hZeroinR); k>= 1;k--)
1662  {//Transforming the different homogeneous parts of h into polys in K[theta]
1663    mons = ideal();
1664    for(i = 1; i<=size(hZeroinR[k][2]);i++)
1665    {//Putting the monomials in a list
1666      mons[size(mons)+1] = hZeroinR[k][2][i];
1667    }//Putting the monomials in a list
1668    setring(KTheta);
1669    mons = thetamap(mons);
1670    for (i = 1; i<=size(mons);i++)
1671    {//transforming the monomials as monomials in theta
1672      entry = leadcoef(mons[i]);
1673      for (j = 0; j<leadexp(mons[i])[2];j++)
1674      {
1675        entry = entry * (theta-j);
1676      }
1677      mons[i] = entry;
1678    }//transforming the monomials as monomials in theta
1679    hZeroinKTheta = hZeroinKTheta + list(sum(mons));
1680    setring(R);
1681  }//Transforming the different homogeneous parts of h into polys in K[theta]
1682  dbprint(p,dbprintWhitespace+" Done!");
1683  //Making the solutionRing
1684  ring solutionRing = 0,(theta,q(0..(kappa+1)*(m1-ml-1)-1)),lp;
1685  dbprint(p,dbprintWhitespace+" Our solution ring is given by "+ string(solutionRing));
1686  //mapping the different ps and qs and HZeros
1687  dbprint(p,dbprintWhitespace+" Setting up our solution system.");
1688  list ps;
1689  ideal pandqZero = imap(KTheta,pandqZero);
1690  ps[1] = list(n1,pandqZero[1]);
1691  ps[n1-nk+1] = list(nk,pandqZero[2]);
1692  for (i = 2; i<=n1-nk; i++)
1693  {
1694    ps[i] = list(n1-i+1,0);
1695  }
1696  list qs;
1697  qs[1] = list(m1,pandqZero[3]);
1698  qs[m1-ml+1] = list(ml,pandqZero[4]);
1699  for (i = 2; i<=m1-ml; i++)
1700  {
1701    qs[i] = list(m1-i+1,0);
1702    for (j = 0; j<=kappa; j++)
1703    {
1704      qs[i][2] = qs[i][2] + q((i-2)*(kappa+1)+j)*theta^j;
1705    }
1706  }
1707  list hZero = imap(KTheta,hZeroinKTheta);
1708  for (i = 1; i<=size(hZero); i++)
1709  {
1710    hZero[i] = list(n1+m1-i+1,hZero[i]);
1711  }
1712
1713  //writing and solving the system
1714  list lhs;
1715  lhs[1] = list(ps[1][2],1);
1716  list rhs;
1717  rhs[n1-nk+1] = list(ps[size(ps)][2],1);
1718  for (i = 2; i<= n1-nk+1;i++)
1719  {
1720    lhs[i] = list(hZero[i][2],subst(qs[1][2],theta,theta+n1-i+1)*gammaForTheta(ps[i][1],qs[1][1]));
1721    rhs[n1-nk-i+2] = list(hZero[size(hZero)-i+1][2],subst(qs[m1-ml+1][2],theta,
1722                          theta+nk+i-1)*gammaForTheta(ps[n1-nk-i+2][1],qs[m1-ml+1][1]));
1723    for (j = 1; j<i; j++)
1724    {
1725      for (k = 1; k<=size(qs);k++)
1726      {
1727        if(ps[j][1]+qs[k][1] == hZero[i][1])
1728        {
1729          lhs[i][1] = lhs[i][1]*lhs[j][2];
1730          lhs[i][2] = lhs[i][2]*lhs[j][2];
1731          lhs[i][1] = lhs[i][1]-lhs[j][1]*subst(qs[k][2],theta, theta + n1- j +1)
1732            *gammaForTheta(ps[j][1],qs[k][1]);
1733        }
1734        if(ps[n1-nk+2-j][1] + qs[m1-ml+2-k][1] ==hZero[size(hZero)-i+1][1])
1735        {
1736          rhs[n1-nk-i+2][1] = rhs[n1-nk-i+2][1]*rhs[n1-nk+2-j][2];
1737          rhs[n1-nk-i+2][2] = rhs[n1-nk-i+2][2]*rhs[n1-nk+2-j][2];
1738          rhs[n1-nk-i+2][1] = rhs[n1-nk-i+2][1]-rhs[n1-nk+2-j][1]*subst(qs[m1-ml+2-k][2],theta,
1739                             theta + nk -j+1)*gammaForTheta(ps[n1-nk+2-j][1],qs[m1-ml+2-k][1]);
1740        }
1741      }
1742    }
1743  }
1744  list eqs;
1745  poly tempgcd;
1746  poly templhscoeff;
1747  poly temprhscoeff;
1748  for (i = 2; i<=n1-nk;i++)
1749  {
1750    if (gcd(rhs[i][2],lhs[i][2]) == 1)
1751    {
1752      eqs = eqs + list(lhs[i][1]*rhs[i][2] -rhs[i][1]*lhs[i][2]);
1753    }
1754    else
1755    {
1756      tempgcd = gcd(rhs[i][2],lhs[i][2]);
1757      templhscoeff = quotient(rhs[i][2],tempgcd)[1];
1758      temprhscoeff = quotient(lhs[i][2],tempgcd)[1];
1759      eqs = eqs + list(lhs[i][1]*templhscoeff -rhs[i][1]*temprhscoeff);
1760    }
1761  }
1762  matrix tempCoefMatrix;
1763  ideal solutionSystemforqs;
1764  for (i= 1; i<=size(eqs); i++)
1765  {
1766    tempCoefMatrix = coef(eqs[i],theta);
1767    solutionSystemforqs = solutionSystemforqs + ideal(submat(tempCoefMatrix,intvec(2),
1768                                                             intvec(1..ncols(tempCoefMatrix))));
1769  }
1770  dbprint(p,dbprintWhitespace+" Solution system for the coefficients of q is given by:");
1771  dbprint(p,solutionSystemforqs);
1772  option(redSB);
1773  dbprint(p,dbprintWhitespace+" Calculating reduced Groebner Basis of that system.");
1774  solutionSystemforqs = slimgb(solutionSystemforqs);
1775  dbprint(p,dbprintWhitespace+" Done!, the solution for the system is:");
1776  dbprint(p,dbprintWhitespace+string(solutionSystemforqs));
1777  if(vdim(slimgb(solutionSystemforqs+theta))==0)
1778  {//No solution in this case. Return the empty list
1779    dbprint(p,dbprintWhitespace+"The Groebner Basis of the solution system was <1>.");
1780    setring(R);
1781    return(list());
1782  }//No solution in this case. Return the empty list
1783  if(vdim(slimgb(solutionSystemforqs+theta))==-1)
1784  {//My conjecture is that this would never happen
1785    //ERROR("This is an counterexample to your conjecture. We have infinitely many solutions");
1786    //TODO: See, what we would do here
1787    dbprint(p,dbprintWhitespace+"There are infinitely many solution to this system.
1788We will return the empty list.");
1789    setring(R);
1790    return(list());
1791  }//My conjecture is that this would never happen
1792  else
1793  {//We have finitely many solutions
1794    if(vdim(slimgb(solutionSystemforqs+theta))==1)
1795    {//exactly one solution
1796      for (i = 2; i<= size(qs)-1;i++)
1797      {
1798        qs[i][2] = NF(qs[i][2],solutionSystemforqs);
1799      }
1800      setring(R);
1801      map backFromSolutionRing = solutionRing,var(1)*var(2);
1802      list qs = backFromSolutionRing(qs);
1803      list result = list(0);
1804      for (i = 1; i<=size(qs); i++)
1805      {
1806        if (qs[i][1]>0){qs[i][2] = qs[i][2]*var(2)^qs[i][1];}
1807        if (qs[i][1]<0){qs[i][2] = qs[i][2]*var(1)^(-qs[i][1]);}
1808        result[1] = result[1] + qs[i][2];
1809      }
1810      dbprint(p,dbprintWhitespace+"Found one unique solution. Returning the result.");
1811      return(result);
1812    }//exactly one solution
1813    else
1814    {//We have more than one solution, but finitely many
1815      dbprint(p,dbprintWhitespace+"Finitely many, but more than one solution.
1816Right now the conjecture is that this cannot happen.");
1817      dbprint(p,dbprintWhitespace+"Yay to counterexample :-) We can have more than one solution.");
1818      //TODO: Some theoretical work why this can never happan
1819      return(list());
1820    }//We have more than one solution, but finitely many
1821  }//We have finitely many solutions
1822}//proc determineRestOfHomogParts
1823
1824static proc gammaForTheta(int j1,int j2)
1825"
1826INPUT: Two integers j1 and j2
1827OUTPUT: A polynomial in the first variable of the given ring. It calculates the following function:
1828   / 1,                                                if j1,j2>0 or j1,j2 <= 0
1829   | prod_{kappa = 0}^{|j1|-1}(var(1)-kappa),          if j1<0, j2>0, |j1|<=|j2|
1830gamma_{j1,j2}:= <  prod_{kappa = 0}^{j2-1}(var(1)-kappa-|j1|+|j2|),  if j1<0, j2>0, |j1|>|j2|
1831   | prod_{kappa = 1}^{j1}(var(1)+kappa),              if j1>0, j2<0, |j1|<=|j2|
1832   \ prod_{kappa = 1}^{|j2|}(\var(1)+kappa+|j1|-|j2|), if j1>0, j2<0, |j1|>|j2|
1833ASSUMPTION:
1834- Ring has at least one variable
1835"
1836{//gammaForTheta
1837  if (j1<=0 && j2 <=0) {return(1);}
1838  if (j1>=0 && j2 >=0){return(1);}
1839  poly result;
1840  int i;
1841  if (j1<0 && j2>0)
1842  {//case 2 or 3 from description above
1843    if (absValue(j1)<=absValue(j2))
1844    {//Case 2 holds here
1845      result = 1;
1846      for (i = 0;i<absValue(j1);i++)
1847      {
1848        result = result*(var(1)-i);
1849      }
1850      return(result);
1851    }//Case 2 holds here
1852    else
1853    {//Case 3 holds here
1854      result = 1;
1855      for (i = 0; i<j2; i++)
1856      {
1857        result = result*(var(1)-i-absValue(j1)+absValue(j2));
1858      }
1859      return(result);
1860    }//Case 3 holds here
1861  }//case 2 or 3 from description above
1862  else
1863  {//Case 4 or 5 from description above hold
1864    if (absValue(j1)<=absValue(j2))
1865    {//Case 4 holds
1866      result = 1;
1867      for (i = 1; i<=j1; i++)
1868      {
1869        result = result*(var(1)+i);
1870      }
1871      return(result);
1872    }//Case 4 holds
1873    else
1874    {//Case 5 holds
1875      result = 1;
1876      for (i = 1; i<=absValue(j2); i++)
1877      {
1878        result = result*(var(1)+i+absValue(j1)-absValue(j2));
1879      }
1880      return(result);
1881    }//Case 5 holds
1882  }//Case 4 or 5 from description above hold
1883}//gammaForTheta
1884
1885static proc extractHomogeneousDivisors(poly h)
1886"INPUT: A polynomial h in the first Weyl algebra
1887OUTPUT: If h is homogeneous, then all factorizations of h are returned.
1888 If h is inhomogeneous, then a list l is returned whose entries
1889 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h and there
1890 exists an i in {1,...,n}, such that k_i is inhomogeneous and there
1891 is no homogeneous polynomial that divides this k_i neither from the
1892 left nor from the right. All the other entries in k are homogeneous
1893 polynomials.
1894"
1895{//extractHomogeneousDivisors
1896  int p=printlevel-voice+2; // for dbprint
1897  string dbprintWhitespace = "";
1898  int i; int j; int k; int l;
1899  list result;
1900  for (i = 1; i<=voice;i++)
1901  {dbprintWhitespace = dbprintWhitespace + " ";}
1902  if (homogwithorder(h,intvec(-1,1)))
1903  {//given polynomial was homogeneous already
1904    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
1905    result = homogfacFirstWeyl_all(h);
1906    for (i = 1; i<=size(result);i++)
1907    {//removing the first entry (coefficient) from the list result
1908      result[i] = delete(result[i],1);
1909    }//removing the first entry (coefficient) from the list result
1910    return(result);
1911  }//given polynomial was homogeneous already
1912  dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous left divisors extracted");
1913  list leftDivisionPossibilities = extractHomogeneousDivisorsLeft(h);
1914  dbprint(p,dbprintWhitespace+"Done. The result is:");
1915  dbprint(p,leftDivisionPossibilities);
1916  dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous Right divisors extracted");
1917  list rightDivisionPossibilities = extractHomogeneousDivisorsRight(h);
1918  dbprint(p,dbprintWhitespace+"Done. The result is:");
1919  dbprint(p,rightDivisionPossibilities);
1920  list tempList;
1921  dbprint(p,dbprintWhitespace+"Calculating remaining right and left homogeneous divisors");
1922  for (i = 1; i<=size(leftDivisionPossibilities); i++)
1923  {//iterating through the list with extracted left divisors
1924    tempList = extractHomogeneousDivisorsRight
1925      (leftDivisionPossibilities[i][size(leftDivisionPossibilities[i])]);
1926    leftDivisionPossibilities[i] = delete(leftDivisionPossibilities[i],
1927                                          size(leftDivisionPossibilities[i]));
1928    for (j=1;j<=size(tempList);j++)
1929    {//Updating the list for Result
1930      tempList[j] = leftDivisionPossibilities[i] + tempList[j];
1931    }//Updating the list for Result
1932    result = result + tempList;
1933  }//iterating through the list with extracted left divisors
1934  for (i = 1; i<=size(rightDivisionPossibilities); i++)
1935  {//iterating through the list with extracted left divisors
1936    tempList = extractHomogeneousDivisorsLeft(rightDivisionPossibilities[i][1]);
1937    rightDivisionPossibilities[i] = delete(rightDivisionPossibilities[i],1);
1938    for (j=1;j<=size(tempList);j++)
1939    {//Updating the list for Result
1940      tempList[j] =  tempList[j]+rightDivisionPossibilities[i];
1941    }//Updating the list for Result
1942    result = result + tempList;
1943  }//iterating through the list with extracted left divisors
1944  dbprint(p,dbprintWhitespace+"Done");
1945  int posInhomog;
1946  poly hath = 1;
1947  list tempResult;
1948  dbprint(p,dbprintWhitespace+"Checking if we can swap left resp. right
1949divisors and updating result.");
1950  for (i = 1; i<= size(result); i++)
1951  {//Checking if we can swap left resp. right divisors
1952    for (j = 1; j<=size(result[i]);j++)
1953    {//finding the position of the inhomogeneous element in the list
1954      if(!homogwithorder(result[i][j],intvec(-1,1)))
1955      {
1956        posInhomog = j;
1957        break;
1958      }
1959    }//finding the position of the inhomogeneous element in the list
1960    hath = result[i][posInhomog];
1961    for(j=posInhomog-1;j>=1;j--)
1962    {
1963      hath = result[i][j]*hath;
1964      tempList = extractHomogeneousDivisorsRight(hath);
1965      if(size(tempList[1])==1)
1966      {//We could not swap this element to the right
1967        break;
1968      }//We could not swap this element to the right
1969      dbprint(p,dbprintWhitespace+"A swapping (left) of an element was possible");
1970      for(k = 1; k<=size(tempList);k++)
1971      {
1972        tempResult = insert(tempResult,result[i]);
1973        for (l = j;l<=posInhomog;l++)
1974        {
1975          tempResult[1] = delete(tempResult[1],j);
1976        }
1977        for (l = size(tempList[k]);l>=1;l--)
1978        {
1979          tempResult[1] = insert(tempResult[1],tempList[k][l],j-1);
1980        }
1981      }
1982    }
1983    hath = result[i][posInhomog];
1984    for(j=posInhomog+1;j<=size(result[i]);j++)
1985    {
1986      hath = hath*result[i][j];
1987      tempList = extractHomogeneousDivisorsLeft(hath);
1988      if(size(tempList[1])==1)
1989      {//We could not swap this element to the right
1990        break;
1991      }//We could not swap this element to the right
1992      dbprint(p,dbprintWhitespace+"A swapping (right) of an element was possible");
1993      for(k = 1; k<=size(tempList);k++)
1994      {
1995        tempResult = insert(tempResult,result[i]);
1996        for (l=posInhomog; l<=j;l++)
1997        {
1998          tempResult[1] = delete(tempResult[1],posInhomog);
1999        }
2000        for (l = size(tempList[k]);l>=1;l--)
2001        {
2002          tempResult[1] = insert(tempResult[1],tempList[k][l],posInhomog-1);
2003        }
2004      }
2005    }
2006  }//Checking if we can swap left resp. right divisors
2007  result = result + tempResult;
2008  result = delete_dublicates_noteval(result);
2009  return(result);
2010}//extractHomogeneousDivisors
2011
2012static proc extractHomogeneousDivisorsLeft(poly h)
2013"INPUT: A polynomial h in the first Weyl algebra
2014OUTPUT: If h is homogeneous, then all factorizations of h are returned.
2015 If h is inhomogeneous, then a list l is returned whose entries
2016 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h.
2017 The entry k_n is inhomogeneous and has no other homogeneous
2018 left divisors any more.
2019 All the other entries in k are homogeneous
2020 polynomials.
2021"
2022{//extractHomogeneousDivisorsLeft
2023  int p=printlevel-voice+2; // for dbprint
2024  string dbprintWhitespace = "";
2025  int i;int j; int k;
2026  list result;
2027  poly hath;
2028  list recResult;
2029  map invo = basering,-var(1),var(2);
2030  for (i = 1; i<=voice;i++)
2031  {dbprintWhitespace = dbprintWhitespace + " ";}
2032  if (homogwithorder(h,intvec(-1,1)))
2033  {//given polynomial was homogeneous already
2034    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
2035    result = homogfacFirstWeyl_all(h);
2036    for (i = 1; i<=size(result);i++)
2037    {//removing the first entry (coefficient) from the list result
2038      result[i] = delete(result[i],1);
2039    }//removing the first entry (coefficient) from the list result
2040    return(result);
2041  }//given polynomial was homogeneous already
2042  list hlist = homogDistribution(h);
2043  dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands.");
2044  for (i = 1; i<= size(hlist); i++)
2045  {
2046    hlist[i] = homogfacFirstWeyl_all(hlist[i][2]);
2047    if (size(hlist[i][1])==1)
2048    {//One homogeneous part just has a trivial factorization
2049      if(hlist[i][1][1] == 0)
2050      {
2051        hlist = delete(hlist,i);
2052        continue;
2053      }
2054      else
2055      {
2056        return(list(list(h)));
2057      }
2058    }//One homogeneous part just has a trivial factorization
2059  }
2060  dbprint(p,dbprintWhitespace+ " Done.");
2061  dbprint(p,dbprintWhitespace+ " Trying to find Left divisors");
2062  list alreadyConsideredCandidates;
2063  poly candidate;
2064  int isCandidate;
2065  for (i = 1; i<=size(hlist[1]);i++)
2066  {//Finding candidates for homogeneous left divisors of h
2067    candidate = hlist[1][i][2];
2068    isCandidate = 0;
2069    for (j=1;j<=size(alreadyConsideredCandidates);j++)
2070    {
2071      if(alreadyConsideredCandidates[j] == candidate)
2072      {
2073        isCandidate =1;
2074        break;
2075      }
2076    }
2077    if(isCandidate)
2078    {
2079      i++;
2080      continue;
2081    }
2082    else
2083    {
2084      alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate);
2085    }
2086    dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous left divisor");
2087    for (j = 2; j<=size(hlist);j++)
2088    {//Iterating through the other homogeneous parts
2089      isCandidate = 0;
2090      for(k=1; k<=size(hlist[j]);k++)
2091      {
2092        if(hlist[j][k][2]==candidate)
2093        {
2094          isCandidate = 1;
2095          break;
2096        }
2097      }
2098      if(!isCandidate)
2099      {
2100        break;
2101      }
2102    }//Iterating through the other homogeneous parts
2103    if(isCandidate)
2104    {//candidate was really a left divisor
2105      dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous left divisor");
2106      hath = involution(lift(involution(candidate,invo),involution(h,invo))[1,1],invo);
2107      recResult = extractHomogeneousDivisorsLeft(hath);
2108      for (j = 1; j<=size(recResult); j++)
2109      {
2110        recResult[j] = insert(recResult[j],candidate);
2111      }
2112      result = result + recResult;
2113    }//Candidate was really a left divisor
2114  }//Finding candidates for homogeneous left divisors of h
2115  if (size(result)==0)
2116  {
2117    return(list(list(h)));
2118  }
2119  return(result);
2120}//extractHomogeneousDivisorsLeft
2121
2122static proc extractHomogeneousDivisorsRight(poly h)
2123"INPUT: A polynomial h in the first Weyl algebra
2124OUTPUT: If h is homogeneous, then all factorizations of h are returned.
2125 If h is inhomogeneous, then a list l is returned whose entries
2126 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h.
2127 The entry k_1 is inhomogeneous and has no other homogeneous
2128 right divisors any more.
2129 All the other entries in k are homogeneous
2130 polynomials.
2131"
2132{//extractHomogeneousDivisorsRight
2133  int p=printlevel-voice+2; // for dbprint
2134  string dbprintWhitespace = "";
2135  int i;int j; int k;
2136  list result;
2137  poly hath;
2138  list recResult;
2139  for (i = 1; i<=voice;i++)
2140  {dbprintWhitespace = dbprintWhitespace + " ";}
2141  if (homogwithorder(h,intvec(-1,1)))
2142  {//given polynomial was homogeneous already
2143    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
2144    result = homogfacFirstWeyl_all(h);
2145    for (i = 1; i<=size(result);i++)
2146    {//removing the first entry (coefficient) from the list result
2147      result[i] = delete(result[i],1);
2148    }//removing the first entry (coefficient) from the list result
2149    return(result);
2150  }//given polynomial was homogeneous already
2151  list hlist = homogDistribution(h);
2152  dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands.");
2153  for (i = 1; i<= size(hlist); i++)
2154  {
2155    hlist[i] = homogfacFirstWeyl_all(hlist[i][2]);
2156    if (size(hlist[i][1])==1)
2157    {//One homogeneous part just has a trivial factorization
2158      if(hlist[i][1][1] == 0)
2159      {
2160        hlist = delete(hlist,i);
2161        continue;
2162      }
2163      else
2164      {
2165        return(list(list(h)));
2166      }
2167    }//One homogeneous part just has a trivial factorization
2168  }
2169  dbprint(p,dbprintWhitespace+ " Done.");
2170  dbprint(p,dbprintWhitespace+ " Trying to find right divisors");
2171  list alreadyConsideredCandidates;
2172  poly candidate;
2173  int isCandidate;
2174  for (i = 1; i<=size(hlist[1]);i++)
2175  {//Finding candidates for homogeneous left divisors of h
2176    candidate = hlist[1][i][size(hlist[1][i])];
2177    isCandidate = 0;
2178    for (j=1;j<=size(alreadyConsideredCandidates);j++)
2179    {
2180      if(alreadyConsideredCandidates[j] == candidate)
2181      {
2182        isCandidate =1;
2183        break;
2184      }
2185    }
2186    if(isCandidate)
2187    {
2188      i++;
2189      continue;
2190    }
2191    else
2192    {
2193      alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate);
2194    }
2195    dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous r-divisor");
2196    for (j = 2; j<=size(hlist);j++)
2197    {//Iterating through the other homogeneous parts
2198      isCandidate = 0;
2199      for(k=1; k<=size(hlist[j]);k++)
2200      {
2201        if(hlist[j][k][size(hlist[j][k])]==candidate)
2202        {
2203          isCandidate = 1;
2204          break;
2205        }
2206      }
2207      if(!isCandidate)
2208      {
2209        break;
2210      }
2211    }//Iterating through the other homogeneous parts
2212    if(isCandidate)
2213    {//candidate was really a left divisor
2214      dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous right divisor");
2215      hath = lift(candidate,h)[1,1];
2216      recResult = extractHomogeneousDivisorsRight(hath);
2217      for (j = 1; j<=size(recResult); j++)
2218      {
2219        recResult[j] = insert(recResult[j],candidate,size(recResult[j]));
2220      }
2221      result = result + recResult;
2222    }//Candidate was really a left divisor
2223  }//Finding candidates for homogeneous left divisors of h
2224  if (size(result)==0)
2225  {
2226    result = list(list(h));
2227  }
2228  return(result);
2229}//extractHomogeneousDivisorsRight
2230
2231static proc fromZeroHomogToThetaPoly(poly h)
2232"
2233//DEPRECATED
2234INPUT: A polynomial h in the first Weyl algebra, homogeneous of degree 0
2235 OUTPUT: The ring Ktheta with a polynomial result representing h as polynomial in
2236  theta
2237 ASSUMPTIONS:
2238- h is homogeneous of degree 0 with respect to the [-1,1] weight vector
2239- The basering is the first Weyl algebra
2240"
2241{//proc fromZeroHomogToThetaPoly
2242  int i; int j;
2243  list mons;
2244  if(!homogwithorder(h,intvec(-1,1)))
2245  {//Input not a homogeneous polynomial causes an error
2246    ERROR("The input was not a homogeneous polynomial");
2247  }//Input not a homogeneous polynomial causes an error
2248  if(deg(h,intvec(-1,1))!=0)
2249  {//Input does not have degree 0
2250    ERROR("The input did not have degree 0");
2251  }//Input does not have degree 0
2252  for(i = 1; i<=size(h);i++)
2253  {//Putting the monomials in a list
2254    mons = mons+list(h[i]);
2255  }//Putting the monomials in a list
2256  ring KTheta = 0,(x,y,theta),dp;
2257  setring KTheta;
2258  map thetamap = r,x,y;
2259  list mons = thetamap(mons);
2260  poly entry;
2261  for (i = 1; i<=size(mons);i++)
2262  {//transforming the monomials as monomials in theta
2263    entry = leadcoef(mons[i]);
2264    for (j = 0; j<leadexp(mons[i])[2];j++)
2265    {
2266      entry = entry * (theta-j);
2267    }
2268    mons[i] = entry;
2269  }//transforming the monomials as monomials in theta
2270  poly result = sum(mons);
2271  keepring KTheta;
2272}//proc fromZeroHomogToThetaPoly
2273
2274static proc fromHomogDistributionListToPoly(list l)
2275"returns the corresponding polynomial to a given output of homogDistribution"
2276{//proc fromHomogDistributionListToPoly
2277  poly result = 0;
2278  int i;
2279  for (i = 1; i<=size(l);i++)
2280  {
2281    result = result + l[i][2];
2282  }
2283  return(result);
2284}//proc fromHomogDistributionListToPoly
2285
2286static proc divides(poly p1, poly p2,map invo, list #)
2287  "Tests, whether p1 divides p2 either from left or from right. The involution invo is needed
2288for checking both sides. The optional argument is needed in order to also return the other factor.
2289RETURN: If no optional argument is given, it will just return 1 or 0.
2290 Otherwise a list with at least one element
2291 Case 1: p1 does not divide p2 from any side. Then the output will be the empty list.
2292 Case 2: p2 does divide p2 from one side at least.
2293  Then it returns a list with tuples p,q, such that p or q equals p1 and
2294  pq = p2.
2295ASSUMPTIONS: - The map invo is an involution on the basering."
2296{//proc divides
2297  list result = list();
2298  poly tempfactor;
2299  if (involution(reduce(involution(p2,invo),std(involution(ideal(p1),invo))),invo)==0)
2300  {//p1 is a left divisor
2301    if(size(#)==0){return(1);}
2302    tempfactor = involution(lift(involution(p1,invo),involution(p2,invo))[1,1],invo);
2303    if (leadcoef(p1)<0 && leadcoef(tempfactor)<0)
2304    {//both have a negative leading coefficient
2305      result = result +list(list(1,-p1, -tempfactor));
2306    }//both have a negative leading coefficient
2307    else
2308    {
2309      if (leadcoef(p1)*leadcoef(tempfactor)<0)
2310      {//One of them has a negative leading coefficient
2311        if (leadcoef(p1)<0)
2312        {
2313          result = result + list(list(-1,-p1,tempfactor));
2314        }
2315        else
2316        {
2317          result = result + list(list(-1,p1,-tempfactor));
2318        }
2319      }//One of them has a negative leading coefficient
2320      else
2321      {//no negative coefficient at all
2322        result = result + list(list(1,p1,tempfactor));
2323      }//no negative coefficient at all
2324    }
2325  }//p1 is a left divisor
2326
2327  if (reduce(p2,std(ideal(p1))) == 0)
2328  {//p1 is  a right divisor
2329    if(size(#)==0){return(1);}
2330    tempfactor = lift(p1, p2)[1,1];
2331    if (leadcoef(p1)<0 && leadcoef(tempfactor)<0)
2332    {//both have a negative leading coefficient
2333      result = result +list(list(1, -tempfactor,-p1));
2334    }//both have a negative leading coefficient
2335    else
2336    {
2337      if (leadcoef(p1)*leadcoef(tempfactor)<0)
2338      {//One of them has a negative leading coefficient
2339        if (leadcoef(p1)<0)
2340        {
2341          result = result + list(list(-1,tempfactor,-p1));
2342        }
2343        else
2344        {
2345          result = result + list(list(-1,-tempfactor,p1));
2346        }
2347      }//One of them has a negative leading coefficient
2348      else
2349      {//no negative coefficient at all
2350        result = result + list(list(1,tempfactor,p1));
2351      }//no negative coefficient at all
2352    }
2353  }//p1 is already a right divisor
2354  if (size(#)==0){return(0);}
2355  return(result);
2356}//proc divides
2357
2358static proc computeCombinationsMinMaxHomog(poly h)
2359"Input: A polynomial h in the first Weyl Algebra
2360Output: Combinations of the form (p_max + p_min)(q_max + q_min), such that p_max, p_min,
2361 q_max and q_min are homogeneous and p_max*q_max equals the maximal homogeneous
2362 part in h, and p_max * q_max equals the minimal homogeneous part in h.
2363 h is not homogeneous.
2364"{//proc computeCombinationsMinMaxHomog
2365  intvec ivm11 = intvec(-1,1);
2366  intvec iv11 = intvec(1,1);
2367  intvec iv10 = intvec(1,0);
2368  intvec iv01 = intvec(0,1);
2369  intvec iv1m1 = intvec(1,-1);
2370  poly maxh = jet(h,deg(h,ivm11),ivm11)-jet(h,deg(h,ivm11)-1,ivm11);
2371  poly minh = jet(h,deg(h,iv1m1),iv1m1)-jet(h,deg(h,iv1m1)-1,iv1m1);
2372  list f1 = homogfacFirstWeyl_all(maxh);
2373  list f2 = homogfacFirstWeyl_all(minh);
2374  list result = list();
2375  int i;
2376  int j;
2377  //TODOs: Nur die kombinieren, die auch wirklich eine brauchbare Kombination sind.
2378  //Also, wir duerfen nicht aus den Geraden herausfliegen
2379  list pqmax = list();
2380  list pqmin = list();
2381  list tempList = list();
2382
2383  //First, we are going to deal with our most hated guys: The Coefficients.
2384  //
2385  list coeffTuplesMax = getAllCoeffTuplesComb(factorizeInt(number(f1[1][1])));
2386  //We can assume without loss of generality, that p_max has a
2387  //nonnegative leading coefficient
2388  for (i = 1; i<=size(coeffTuplesMax);i++)
2389  {//Deleting all tuples with negative entries for p_max
2390    if (coeffTuplesMax[i][1]<0)
2391    {
2392      coeffTuplesMax = delete(coeffTuplesMax,i);
2393      continue;
2394    }
2395  }//Deleting all tuples with negative entries for p_max
2396  list coeffTuplesMin = getAllCoeffTuplesComb(factorizeInt(number(f2[1][1])));
2397
2398  //Now, we will be actally dealing with the Combinations.
2399  //Let's start with the pqmax
2400  if (size(f1[1]) == 1)
2401  {//the maximal homogeneous factor is a constant
2402    pqmax = coeffTuplesMax;
2403  }//the maximal homogeneous factor is a constant
2404  else
2405  {//the maximal homogeneous factor is not a constant
2406    for (i = 1; i<=size(f1); i++)
2407    {//We can forget about the first coefficient now. Therefore we will delete him from the list.
2408      f1[i] = delete(f1[i],1);
2409      if(size(f1[i])==1)
2410      {//trivial thing
2411        for (j = 1; j<=size(coeffTuplesMax); j++)
2412        {
2413          pqmax = pqmax + list(list(coeffTuplesMax[j][1],coeffTuplesMax[j][2]*f1[i][1]));
2414          pqmax = pqmax + list(list(coeffTuplesMax[j][1]*f1[i][1],coeffTuplesMax[j][2]));
2415        }
2416        f1 = delete(f1,i);
2417        continue;
2418      }//trivial thing
2419    }//We can forget about the first coefficient now. Therefore we will delete him from the list.
2420    if (size(f1)>0)
2421    {
2422      tempList = getAllCombOfHomogFact(f1);
2423      for (i = 1; i<=size(tempList); i++)
2424      {//Every combination combined with the coefficient possibilities
2425        for (j = 1; j<=size(coeffTuplesMax); j++)
2426        {//iterating through the possible coefficient choices
2427          pqmax = pqmax + list(list(coeffTuplesMax[j][1]*tempList[i][1],
2428                                    coeffTuplesMax[j][2]*tempList[i][2]));
2429        }//iterating through the possible coefficient choices
2430      }//Every combination combined with the coefficient possibilities
2431      for (i = 1; i<=size(coeffTuplesMax); i++)
2432      {
2433        pqmax = pqmax + list(list(coeffTuplesMax[i][1],maxh/coeffTuplesMax[i][1]));
2434        pqmax = pqmax + list(list(maxh/coeffTuplesMax[i][2],coeffTuplesMax[i][2]));
2435      }
2436    }
2437  }//the maximal homogeneous factor is not a constant
2438  //Now we go to pqmin
2439  if (size(f2[1]) == 1)
2440  {//the minimal homogeneous factor is a constant
2441    pqmin = coeffTuplesMin;
2442  }//the minimal homogeneous factor is a constant
2443  else
2444  {//the minimal homogeneous factor is not a constant
2445    for (i = 1; i<=size(f2); i++)
2446    {//We can forget about the first coefficient now. Therefore we will delete him from the list.
2447      f2[i] = delete(f2[i],1);
2448      if(size(f2[i])==1)
2449      {//trivial thing
2450        for (j = 1; j<=size(coeffTuplesMin); j++)
2451        {
2452          pqmin = pqmin + list(list(coeffTuplesMin[j][1],coeffTuplesMin[j][2]*f2[i][1]));
2453          pqmin = pqmin + list(list(coeffTuplesMin[j][1]*f2[i][1],coeffTuplesMin[j][2]));
2454        }
2455        f2 = delete(f2,i);
2456        continue;
2457      }
2458    }//We can forget about the first coefficient now. Therefore we will delete him from the list.
2459    if(size(f2)>0)
2460    {
2461      tempList = getAllCombOfHomogFact(f2);
2462      for (i = 1; i<=size(tempList); i++)
2463      {//Every combination combined with the coefficient possibilities
2464        for (j = 1; j<=size(coeffTuplesMin); j++)
2465        {//iterating through the possible coefficient choices
2466          pqmin = pqmin + list(list(coeffTuplesMin[j][1]*tempList[i][1],
2467                                    coeffTuplesMin[j][2]*tempList[i][2]));
2468        }//iterating through the possible coefficient choices
2469      }//Every combination combined with the coefficient possibilities
2470      for (i = 1; i<=size(coeffTuplesMin); i++)
2471      {
2472        pqmin = pqmin + list(list(coeffTuplesMin[i][1],minh/coeffTuplesMin[i][1]));
2473        pqmin = pqmin + list(list(minh/coeffTuplesMin[i][2],coeffTuplesMin[i][2]));
2474      }
2475    }
2476  }//the minimal homogeneous factor is not a constant
2477
2478  //and now we combine them together to obtain all possibilities.
2479  for (i = 1; i<=size(pqmax); i++)
2480  {//iterate over the maximal homogeneous combination possibilities
2481    for (j = 1; j<=size(pqmin); j++)
2482    {//iterate over the minimal homogeneous combiniation possibilities
2483      if (deg(pqmax[i][1], ivm11)>=deg(pqmin[j][1],ivm11) and deg(pqmax[i][2],
2484                                                                  ivm11)>=deg(pqmin[j][2],ivm11))
2485      {
2486        if (pqmax[i][1]+pqmin[j][1]!=0 and pqmax[i][2]+pqmin[j][2]!=0)
2487        {
2488
2489          if (deg(h,ivm11)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),ivm11))
2490          {
2491            j++;
2492            continue;
2493          }
2494          if (deg(h,iv1m1)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),iv1m1))
2495          {
2496            j++;
2497            continue;
2498          }
2499          result = result +list(list(pqmax[i][1]+pqmin[j][1],pqmax[i][2]+pqmin[j][2]));
2500        }
2501      }
2502    }//iterate over the minimal homogeneous combiniation possibilities
2503  }//iterate over the maximal homogeneous combination possibilities
2504  //Now deleting double entries
2505  result = delete_dublicates_noteval(result);
2506  return(result);
2507}//proc computeCombinationsMinMaxHomog
2508
2509static proc getAllCombOfHomogFact(list l)
2510"Gets called in computeCombinationsMinMaxHomog. It gets a list of different homogeneous
2511factorizations of
2512one homogeneous polynomial and returns the possibilities to combine them into two factors.
2513Assumptions:
2514- The list does not contain the first coefficient.
2515- The list contains at least one list with two elements."
2516{//proc getAllCombOfHomogFact
2517  list result;
2518  list leftAndRightHandSides;
2519  int i; int j;
2520  list tempset;
2521  if (size(l)==1 and size(l[1])==2)
2522  {
2523    result = result + list(list(l[1][1],l[1][2]));
2524    return(result);
2525  }
2526  leftAndRightHandSides = getPossibilitiesForRightSides(l);
2527  for (i = 1; i<=size(leftAndRightHandSides); i++)
2528  {
2529    result =result+list(list(leftAndRightHandSides[i][1],product(leftAndRightHandSides[i][2][1])));
2530    //tidy up the right hand sides, because, if it is just one irreducible factor, we are done
2531    for (j = 1; j<=size(leftAndRightHandSides[i][2]);j++)
2532    {//Tidy up right hand sides
2533      if (size(leftAndRightHandSides[i][2][j])<2)
2534      {//Element can be dismissed
2535        leftAndRightHandSides[i][2] = delete(leftAndRightHandSides[i][2],j);
2536        continue;
2537      }//Element can be dismissed
2538    }//Tidy up right hand sides
2539    if (size(leftAndRightHandSides[i][2])>0)
2540    {
2541      tempset = getAllCombOfHomogFact(leftAndRightHandSides[i][2]);
2542      for (j = 1; j<=size(tempset);j++)
2543      {//multiplying the first factor with the left hand side
2544        result = result + list(list(leftAndRightHandSides[i][1]*tempset[j][1],tempset[j][2]));
2545      }//multiplying the first factor with the left hand side
2546    }
2547  }
2548  return(result);
2549}//proc getAllCombOfHomogFact
2550
2551static proc getPossibilitiesForRightSides(list l)
2552"Given a list of different factorizations l, this function returns a list of the form
2553(a,{(a_2,...,a_n)| (a,a_2,...,a_n) in A})"
2554{//getPossibilitiesForRightSide
2555  list templ = l;
2556  list result;
2557  poly firstElement;
2558  list rightSides;
2559  list tempRightSide;
2560  int i; int j;
2561  while (size(templ)>0)
2562  {
2563    firstElement = templ[1][1];
2564    rightSides = list();
2565    for (i = 1; i<= size(templ); i++)
2566    {
2567      if (templ[i][1] == firstElement)
2568      {//save the right sides
2569        tempRightSide = list();
2570        for (j = 2; j<=size(templ[i]);j++)
2571        {
2572          tempRightSide = tempRightSide + list(templ[i][j]);
2573        }
2574        if (size(tempRightSide)!=0)
2575        {
2576          rightSides = rightSides + list(tempRightSide);
2577        }
2578        templ = delete(templ,i);
2579        continue;
2580      }//save the right sides
2581    }
2582    result = result + list(list(firstElement,rightSides));
2583  }
2584  return(result);
2585}//getPossibilitiesForRightSide
2586
2587static proc getAllCoeffTuplesComb(list l)"
2588Given the output of factorizeInt ((a_1,...,a_n),(i_1,...,i_n)) , it returns all possible tuples
2589of the set {(a,b) | There exists an real N!=emptyset subset of {1,...,n}, such that
2590a = prod_{i \in N}a_i, b=prod_{i \not\in N} a_i}
2591Assumption: The list is sorted from smallest integer to highest.
2592- it is not the factorization of 0.
2593"
2594{//proc getAllCoeffTuplesComb
2595  list result;
2596  if (l[1][1] == 0)
2597  {
2598    ERROR("getAllCoeffTuplesComb: Zero Coefficients as leading and Tail Coeffs?
2599That is not possible. Something went wrong.");
2600  }
2601  if (size(l[1]) == 1)
2602  {//Trivial Factorization, just 1
2603    if (l[1][1] == 1)
2604    {
2605      return(list(list(1,1),list(-1,-1)));
2606    }
2607    else
2608    {
2609      return(list(list(-1,1),list(1,-1)));
2610    }
2611  }//Trivial Factorization, just 1
2612  if (size(l[1]) == 2 and l[2][2]==1)
2613  {//Just a prime number
2614    if (l[1][1] == 1)
2615    {
2616      result = list(list(l[1][2],1),list(1,l[1][2]));
2617      result = result + list(list(-l[1][2],-1),list(-1,-l[1][2]));
2618      return(result);
2619    }
2620    else
2621    {
2622      result = list(list(l[1][2],-1),list(1,-l[1][2]));
2623      result = result + list(list(-l[1][2],1),list(-1,l[1][2]));
2624      return(result);
2625    }
2626  }//Just a prime number
2627  //Now comes the interesting case: a product of primes
2628  list tempPrimeFactors;
2629  list tempPowersOfThem;
2630  int i;
2631  for (i = 2; i<=size(l[1]);i++)
2632  {//Removing the starting 1 or -1 to get the N's
2633    tempPrimeFactors[i-1] = l[1][i];
2634    tempPowersOfThem[i-1] = l[2][i];
2635  }//Removing the starting 1 or -1 to get the N's
2636  list Ns = getAllSubsetsN(list(tempPrimeFactors,tempPowersOfThem));
2637  list tempTuples;
2638  number productOfl = multiplyFactIntOutput(l);
2639  if (productOfl<0){productOfl = -productOfl;}
2640  tempTuples = tempTuples + list(list(1,productOfl),list(productOfl,1));
2641  for (i = 1; i<=size(Ns); i++)
2642  {
2643    if (productOfl/Ns[i]>Ns[i])
2644    {//TODO: BEWEISEN, dass das die einzigen Combos sind
2645      tempTuples = tempTuples + list(list(Ns[i],productOfl/Ns[i]),list(productOfl/Ns[i],Ns[i]));
2646    }//TODO: BEWEISEN, dass das die einzigen Combos sind
2647    if (productOfl/Ns[i]==Ns[i])
2648    {
2649      tempTuples = tempTuples + list(list(Ns[i],Ns[i]));
2650    }
2651  }
2652  //And now, it just remains to get the -1s and 1-s correctly to the tuples
2653  list tempEntry;
2654  if (l[1][1] == 1)
2655  {
2656    for (i = 1; i<=size(tempTuples);i++)
2657    {//Adding everything to result
2658      tempEntry = tempTuples[i];
2659      result = result + list(tempEntry);
2660      result = result + list(list(-tempEntry[1], -tempEntry[2]));
2661    }//Adding everyThing to Result
2662  }
2663  else
2664  {
2665    for (i = 1; i<=size(tempTuples);i++)
2666    {//Adding everything to result
2667      tempEntry = tempTuples[i];
2668      result = result + list(list(tempEntry[1],-tempEntry[2]));
2669      result = result + list(list(-tempEntry[1], tempEntry[2]));
2670    }//Adding everyThing to Result
2671  }
2672  return(result);
2673}//proc getAllCoeffTuplesComb
2674
2675static proc contains(list l, int elem)
2676"Assumption: l is sorted"
2677{//Binary Search in list
2678  if (size(l)<=1)
2679  {
2680    if(size(l) == 0){return(0);}
2681    if (l[1]!=elem){return(0);}
2682    else{return(1);}
2683  }
2684  int imax = size(l);
2685  int imin = 1;
2686  int imid;
2687  while(imax >= imin)
2688  {
2689    imid = (imin + imax)/2;
2690    if (l[imid] == elem){return(1);}
2691    if (l[imid] <elem) {imin = imid +1;}
2692    else{imax = imid -1;}
2693  }
2694  return(0)
2695    }//Binary Search in list
2696
2697static proc getAllSubsetsN(list l)
2698"
2699Assumptions:
2700- The list is containing two lists. They can be assumed to be outputs of the function
2701factorizeInt. They have at least one entry. If it is exactly one entry, the second intvec should
2702contain a value at least 2.
2703  "
2704{
2705  list primeFactors=l[1];
2706  list powersOfThem = l[2];
2707  int i;int j;
2708  //Casting the entries to be numbers
2709  for (i=1; i<=size(primeFactors); i++)
2710  {
2711    primeFactors[i] = number(primeFactors[i]);
2712    powersOfThem[i] = number(powersOfThem[i]);
2713  }
2714
2715  //Done
2716  list result;
2717  list tempPrimeFactors;
2718  list tempPowersOfThem;
2719  list tempset;
2720  if (sum(powersOfThem) <=2)
2721  {//Easy Case
2722    return(list(primeFactors[1]));
2723  }//Easy Case
2724  if (size(primeFactors)==1)
2725  {//Also Easy Case
2726    for (j = 1; j<powersOfThem[1]; j++)
2727    {
2728      result = result + list(primeFactors[1]^j);
2729    }
2730    return(result);
2731  }//Also Easy Case
2732  for (i = 1; i<= size(primeFactors); i++)
2733  {//Going through every entry
2734    result = result + list(primeFactors[i]);
2735    if (i == size(primeFactors))
2736    {
2737      for (j = 1;j<powersOfThem[i];j++)
2738      {
2739        result = result + list (primeFactors[i]^j);
2740      }
2741      break;
2742    }
2743    if (powersOfThem[i]==1)
2744    {
2745      for (j = i+1;j<=size(primeFactors);j++)
2746      {
2747        tempPrimeFactors[j-i] = primeFactors[j];
2748        tempPowersOfThem[j-i] = powersOfThem[j];
2749      }
2750    }
2751    else
2752    {
2753      for (j = i; j<=size(primeFactors);j++)
2754      {
2755        tempPrimeFactors[j-i+1] = primeFactors[j];
2756        tempPowersOfThem[j-i+1] = powersOfThem[j];
2757        tempPowersOfThem[1] = tempPowersOfThem[1]-1;
2758      }
2759    }
2760    tempset = getAllSubsetsN(list(tempPrimeFactors,tempPowersOfThem));
2761    for (j = 1; j<=size(tempset); j++)
2762    {
2763      result = result +list((tempset[j])*(primeFactors[i]));
2764    }
2765  }//Going through every entry
2766  result = sort(result)[1];
2767  result = delete_dublicates_noteval(result);
2768  return(result);
2769}
2770
2771static proc multiplyFactIntOutput(list l)
2772"Given the output of factorizeInt, this method computes the product of it."
2773{//proc multiplyFactIntOutput
2774  int i;
2775  number result = 1;
2776  for (i = 1; i<=size(l[1]); i++)
2777  {
2778    result = result*(l[1][i])^(l[2][i]);
2779  }
2780  return(result);
2781}//proc multiplyFactIntOutput
2782
2783static proc fromListToIntvec(list l)
2784"Converter from List to intvec"
2785{
2786  intvec result; int i;
2787  for (i = 1; i<=size(l); i++)
2788  {
2789    result[i] = l[i];
2790  }
2791  return(result);
2792}
2793
2794static proc fromIntvecToList(intvec l)"
2795Converter from intvec to list"
2796{//proc fromIntvecToList
2797  list result = list();
2798  int i;
2799  for (i = size(l); i>=1; i--)
2800  {
2801    result = insert(result, l[i]);
2802  }
2803  return(result);
2804}//proc fromIntvecToList
2805
2806
2807static proc factorizeInt(number n)
2808"Given an integer n, factorizeInt computes its factorization. The output is a list
2809containing two intvecs. The first contains the prime factors, the second its powers.
2810ASSUMPTIONS:
2811- n is given as integer number
2812"{
2813  if (n==0)
2814  {return(list(list(0),list(1)));}
2815  int i;
2816  list temp = primefactors(n);
2817  if (n<0)
2818  {list result = list(list(-1),list(1));}
2819  else
2820  {list result = list(list(1),list(1));}
2821  result[1] = result[1] + temp[1];
2822  result[2] = result[2] + temp[2];
2823  return(result);
2824}
2825
2826
2827static proc homogDistribution(poly h)
2828"Input: A polynomial in the first Weyl Algebra.
2829  Output: A two-dimensional list of the following form. Every sublist contains exactly two entries.
2830   One for the Z-degree of the corresponding homogeneous part (integer), and the homogeneous
2831   polynomial itself, and those sublists are oredered by ascending degree.
2832   For example a call of homogDistribution(x+d+1) would have the output
2833   [1]:
2834     [1]:
2835       -1
2836     [2]:
2837       x
2838   [2]:
2839     [1]:
2840       0
2841     [2]:
2842       1
2843   [3]:
2844     [1]:
2845       1
2846     [2]:
2847       d
2848"{//homogDistribution
2849  if (h == 0)
2850  {//trivial case where input is 0
2851    return(list(list(0,0)));
2852  }//trivial case where input is 0
2853  if (!isWeyl())
2854  {//Our basering is not the Weyl algebra
2855    ERROR("Ring was not the first Weyl algebra");
2856    return(list());
2857  }//Our basering is not the Weyl algebra
2858  if(nvars(basering)!=2)
2859  {//Our basering is the Weyl algebra, but not the first
2860    ERROR("Ring is not the first Weyl algebra");
2861    return(list());
2862  }//Our basering is the Weyl algebra, but not the first
2863  intvec ivm11 = intvec(-1,1);
2864  intvec iv1m1 = intvec(1,-1);
2865  poly tempH = h;
2866  poly minh;
2867  list result = list();
2868  int nextExpectedDegree = -deg(tempH,iv1m1);
2869  while (tempH != 0)
2870  {
2871    minh = jet(tempH,deg(tempH,iv1m1),iv1m1)-jet(tempH,deg(tempH,iv1m1)-1,iv1m1);
2872    while (deg(minh,ivm11)>nextExpectedDegree)
2873    {//filling empty homogeneous spaces with 0
2874      result = result + list(list(nextExpectedDegree,0));
2875      nextExpectedDegree = nextExpectedDegree +1;
2876    }//filling empty homogeneous spaces with 0
2877    result = result + list(list(deg(minh,ivm11),minh));
2878    tempH = tempH - minh;
2879    nextExpectedDegree = nextExpectedDegree +1;
2880  }
2881  return(result);
2882}//homogDistribution
2883
2884static proc countHomogParts(poly h)
2885"Counts the homogeneous parts of a given polynomial h"
2886{
2887  int i;
2888  list outPutHD = homogDistribution(h);
2889  int result = 0;
2890  for (i = 1; i <=size(outPutHD); i++)
2891  {
2892    if (outPutHD[i][2] != 0){result++;}
2893  }
2894  return(result);
2895}
2896
2897//////////////////////////////////////////////////
2898/////BRANDNEW!!!!////////////////////
2899//////////////////////////////////////////////////
2900
2901//==================================================
2902/*Singular has no way implemented to test polynomials
2903  for homogenity with respect to a weight vector.
2904  The following procedure does exactly this*/
2905static proc homogwithorder(poly h, intvec weights)
2906{//proc homogwithorder
2907  if(size(weights) != nvars(basering))
2908  {//The user does not know how many variables the current ring has
2909    return(0);
2910  }//The user does not know how many variables the current ring has
2911  int i;
2912  int dofp = deg(h,weights); //degree of polynomial
2913  for (i = 1; i<=size(h);i++)
2914  {
2915    if (deg(h[i],weights)!=dofp)
2916    {
2917      return(0);
2918    }
2919  }
2920  return(1);
2921}//proc homogwithorder
2922
2923//==================================================
2924//Testfac: Given a list with different factorizations of
2925// one polynomial, the following procedure checks
2926// whether they all refer to the same polynomial.
2927// If they do, the output will be a list, that contains
2928// the product of each factorization. If not, the empty
2929// list will be returned.
2930// If the optional argument # is given (i.e. the polynomial
2931// which is factorized by the elements of the given list),
2932// then we look, if the entries are factorizations of p
2933// and if not, a list with the products subtracted by p
2934// will be returned
2935proc testNCfac(list l, list #)
2936"USAGE: testNCfac(l[,p,b]); l is a list, p is an optional poly, b is 1 or 0
2937RETURN: Case 1: No optional argument. In this case the output is 1, if the
2938  entries in the given list represent the same polynomial or 0
2939  otherwise.
2940 Case 2: One optional argument p is given. In this case it returns 1,
2941  if all the entries in l are factorizations of p, otherwise 0.
2942 Case 3: Second optional b is given. In this case a list is returned
2943  containing the difference between the product of each entry in
2944  l and p.
2945ASSUME: basering is the first Weyl algebra, the entries of l are polynomials
2946PURPOSE: Checks whether a list of factorizations contains factorizations of
2947  the same element in the first Weyl algebra
2948THEORY: @code{testNCfac} multiplies out each factorization and checks whether
2949 each factorization was a factorization of the same element.
2950@* - if there is only a list given, the output will be 0, if it
2951     does not contain factorizations of the same element. Otherwise the output
2952     will be 1.
2953@* - if there is a polynomial in the second argument, then the procedure checks
2954     whether the given list contains factorizations of this polynomial. If it
2955     does, then the output depends on the third argument. If it is not given,
2956     the procedure will check whether the factorizations in the list
2957     l are associated to this polynomial and return either 1 or 0, respectively.
2958     If the third argument is given, the output will be a list with
2959     the length of the given one and in each entry is the product of one
2960     entry in l subtracted by the polynomial.
2961EXAMPLE: example testNCfac; shows examples
2962SEE ALSO: facFirstWeyl, facSubWeyl, facFirstShift
2963"{//proc testfac
2964  int p = printlevel - voice + 2;
2965  int i;
2966  string dbprintWhitespace = "";
2967  for (i = 1; i<=voice;i++)
2968  {dbprintWhitespace = dbprintWhitespace + " ";}
2969  dbprint(p,dbprintWhitespace + " Checking the input");
2970  if (size(l)==0)
2971  {//The empty list is given
2972    dbprint(p,dbprintWhitespace + " Given list was empty");
2973    return(list());
2974  }//The empty list is given
2975  if (size(#)>2)
2976  {//We want max. two optional arguments
2977    dbprint(p,dbprintWhitespace + " More than two optional arguments");
2978    return(list());
2979  }//We want max. two optional arguments
2980  dbprint(p,dbprintWhitespace + " Done");
2981  list result;
2982  int j;
2983  if (size(#)==0)
2984  {//No optional argument is given
2985    dbprint(p,dbprintWhitespace + " No optional arguments");
2986    int valid = 1;
2987    for (i = size(l);i>=1;i--)
2988    {//iterate over the elements of the given list
2989      if (size(result)>0)
2990      {
2991        if (product(l[i])!=result[size(l)-i])
2992        {
2993          valid = 0;
2994          break;
2995        }
2996      }
2997      result = insert(result, product(l[i]));
2998    }//iterate over the elements of the given list
2999    return(valid);
3000  }//No optional argument is given
3001  else
3002  {
3003    dbprint(p,dbprintWhitespace + " Optional arguments are given.");
3004    int valid = 1;
3005    for (i = size(l);i>=1;i--)
3006    {//iterate over the elements of the given list
3007      if (product(l[i])!=#[1])
3008      {
3009        valid = 0;
3010      }
3011      result = insert(result, product(l[i])-#[1]);
3012    }//iterate over the elements of the given list
3013    if(size(#)==2)
3014    {
3015      dbprint(p,dbprintWhitespace + " A third argument is given. Output is a list now.");
3016      return(result);
3017    }
3018    return(valid);
3019  }
3020}//proc testfac
3021example
3022{
3023  "EXAMPLE:";echo=2;
3024  ring r = 0,(x,y),dp;
3025  def R = nc_algebra(1,1);
3026  setring R;
3027  poly h = (x^2*y^2+1)*(x^2);
3028  def t1 = facFirstWeyl(h);
3029  //fist a correct list
3030  testNCfac(t1);
3031  //now a correct list with the factorized polynomial
3032  testNCfac(t1,h);
3033  //now we put in an incorrect list without a polynomial
3034  t1[3][3] = y;
3035  testNCfac(t1);
3036  // take h as additional input
3037  testNCfac(t1,h);
3038  // take h as additional input and output list of differences
3039  testNCfac(t1,h,1);
3040}
3041//==================================================
3042//Procedure facSubWeyl:
3043//This procedure serves the purpose to compute a
3044//factorization of a given polynomial in a ring, whose subring
3045//is the first Weyl algebra. The polynomial must only contain
3046//the two arguments, which are also given by the user.
3047
3048proc facSubWeyl(poly h, X, D)
3049"USAGE:  facSubWeyl(h,x,y); h, X, D polynomials
3050RETURN: list
3051ASSUME: X and D are variables of a basering, which satisfy DX = XD +1.
3052@* That is,  they generate the copy of the first Weyl algebra in a basering.
3053@* Moreover, h is a polynomial in X and D only.
3054PURPOSE: compute factorizations of the polynomial, which depends on X and D.
3055EXAMPLE: example facSubWeyl; shows examples
3056SEE ALSO: facFirstWeyl, testNCfac, facFirstShift
3057"{
3058  int p = printlevel - voice + 2;
3059  dbprint(p," Start initial Checks of the input.");
3060  // basering can be anything having a Weyl algebra as subalgebra
3061  def @r = basering;
3062  //We begin to check the input for assumptions
3063  // which are: X,D are vars of the basering,
3064  if ( (isVar(X)!=1) || (isVar(D)!=1) || (size(X)>1) || (size(D)>1) ||
3065       (leadcoef(X) != number(1)) || (leadcoef(D) != number(1)) )
3066  {
3067    ERROR("expected pure variables as generators of a subalgebra");
3068  }
3069  // Weyl algebra:
3070  poly w = D*X-X*D-1; // [D,X]=1
3071  poly u = D*X-X*D+1; // [X,D]=1
3072  if (u*w!=0)
3073  {
3074    // that is no combination gives Weyl
3075    ERROR("2nd and 3rd argument do not generate a Weyl algebra");
3076  }
3077  // one of two is correct
3078  int isReverted = 0; // Reverted Weyl if dx=xd-1 holds
3079  if (u==0)
3080  {
3081    isReverted = 1;
3082  }
3083  // else: do nothing
3084  // DONE with assumptions, Input successfully checked
3085  dbprint(p," Successful");
3086  intvec lexpofX = leadexp(X);
3087  intvec lexpofD = leadexp(D);
3088  int varnumX=1;
3089  int varnumD=1;
3090  while(lexpofX[varnumX] != 1)
3091  {
3092    varnumX++;
3093  }
3094  while(lexpofD[varnumD] != 1)
3095  {
3096    varnumD++;
3097  }
3098  /* VL : to add printlevel stuff */
3099  dbprint(p," Change positions of the two variables in the list, if needed");
3100  if (isReverted)
3101  {
3102    ring firstweyl = 0,(var(varnumD),var(varnumX)),dp;
3103    def Firstweyl = nc_algebra(1,1);
3104    setring Firstweyl;
3105    ideal M = 0:nvars(@r);
3106    M[varnumX]=var(2);
3107    M[varnumD]=var(1);
3108    map Q = @r,M;
3109    poly h= Q(h);
3110  }
3111  else
3112  { // that is unReverted
3113    ring firstweyl = 0,(var(varnumX),var(varnumD)),dp;
3114    def Firstweyl = nc_algebra(1,1);
3115    setring Firstweyl;
3116    poly h= imap(@r,h);
3117  }
3118  dbprint(p," Done!");
3119  list result = facFirstWeyl(h);
3120  setring @r;
3121  list result;
3122  if (isReverted)
3123  {
3124    // map swap back
3125    ideal M; M[1] = var(varnumD); M[2] = var(varnumX);
3126    map S = Firstweyl, M;
3127    result = S(result);
3128  }
3129  else
3130  {
3131    // that is unReverted
3132    result = imap(Firstweyl,result);
3133  }
3134  return(result);
3135}//proc facSubWeyl
3136example
3137{
3138  "EXAMPLE:";echo=2;
3139  ring r = 0,(x,y,z),dp;
3140  matrix D[3][3]; D[1,3]=-1;
3141  def R = nc_algebra(1,D); // x,z generate Weyl subalgebra
3142  setring R;
3143  poly h = (x^2*z^2+x)*x;
3144  list fact1 = facSubWeyl(h,x,z);
3145  // compare with facFirstWeyl:
3146  ring s = 0,(z,x),dp;
3147  def S = nc_algebra(1,1); setring S;
3148  poly h = (x^2*z^2+x)*x;
3149  list fact2 = facFirstWeyl(h);
3150  map F = R,x,0,z;
3151  list fact1 = F(fact1); // it is identical to list fact2
3152  testNCfac(fact1); // check the correctness again
3153}
3154//==================================================
3155
3156//==================================================
3157//************From here: Shift-Algebra**************
3158//==================================================
3159//==================================================*
3160//one factorization of a homogeneous polynomial
3161//in the first Shift Algebra
3162static proc homogfacFirstShift(poly h)
3163{//proc homogfacFirstShift
3164  int p=printlevel-voice+2; //for dbprint
3165  def r = basering;
3166  poly hath;
3167  intvec iv01 = intvec(0,1);
3168  int i; int j;
3169  string dbprintWhitespace = "";
3170  for (i = 1; i<=voice;i++)
3171  {dbprintWhitespace = dbprintWhitespace + " ";}
3172  if (!homogwithorder(h,iv01))
3173  {//The given polynomial is not homogeneous
3174    ERROR("The given polynomial is not homogeneous.");
3175    return(list());
3176  }//The given polynomial is not homogeneous
3177  if (h==0)
3178  {
3179    return(list(0));
3180  }
3181  list result;
3182  int m = deg(h,iv01);
3183  dbprint(p,dbprintWhitespace+" exclude the homogeneous part of deg. 0");
3184  if (m>0)
3185  {//The degree is not zero
3186    hath = lift(var(2)^m,h)[1,1];
3187    for (i = 1; i<=m;i++)
3188    {
3189      result = result + list(var(2));
3190    }
3191  }//The degree is not zero
3192  else
3193  {//The degree is zero
3194    hath = h;
3195  }//The degree is zero
3196  ring tempRing = 0,(x),dp;
3197  setring tempRing;
3198  map thetamap = r,x,1;
3199  poly hath = thetamap(hath);
3200  dbprint(p,dbprintWhitespace+" Factorize it using commutative factorization.");
3201  list azeroresult = factorize(hath);
3202  list azeroresult_return_form;
3203  for (i = 1; i<=size(azeroresult[1]);i++)
3204  {//rewrite the result of the commutative factorization
3205    for (j = 1; j <= azeroresult[2][i];j++)
3206    {
3207      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
3208    }
3209  }//rewrite the result of the commutative factorization
3210  setring(r);
3211  map finalmap = tempRing,var(1);
3212  list tempresult = finalmap(azeroresult_return_form);
3213  result = tempresult+result;
3214  return(result);
3215}//proc homogfacFirstShift
3216
3217//==================================================
3218//Computes all possible homogeneous factorizations
3219static proc homogfacFirstShift_all(poly h)
3220{//proc HomogfacFirstShiftAll
3221  int p=printlevel-voice+2; //for dbprint
3222  intvec iv11 = intvec(1,1);
3223  if (deg(h,iv11) <= 0 )
3224  {//h is a constant
3225    return(list(list(h)));
3226  }//h is a constant
3227  def r = basering;
3228  list one_hom_fac; //stands for one homogeneous factorization
3229  int i; int j; int k;
3230  string dbprintWhitespace = "";
3231  for (i = 1; i<=voice;i++)
3232  {dbprintWhitespace = dbprintWhitespace + " ";}
3233  int shiftcounter;
3234  //Compute again a homogeneous factorization
3235  dbprint(p,dbprintWhitespace+" Computing one homog. factorization of the polynomial");
3236  one_hom_fac = homogfacFirstShift(h);
3237  one_hom_fac = delete(one_hom_fac,1);
3238  if (size(one_hom_fac) == 0)
3239  {//there is no homogeneous factorization or the polynomial was not homogeneous
3240    return(list());
3241  }//there is no homogeneous factorization or the polynomial was not homogeneous
3242  dbprint(p,dbprintWhitespace+" Permuting the 0-homogeneous part with the s");
3243  list result = permpp(one_hom_fac);
3244  for (i = 1; i<=size(result);i++)
3245  {
3246    shiftcounter = 0;
3247    for (j = 1; j<=size(result[i]); j++)
3248    {
3249      if (result[i][j]==var(2))
3250      {
3251        shiftcounter++;
3252      }
3253      else
3254      {
3255        result[i][j] = subst(result[i][j], var(1), var(1)-shiftcounter);
3256      }
3257    }
3258    result[i] = insert(result[i],1);
3259  }
3260  dbprint(p,dbprintWhitespace+" Deleting double entries in the resulting list");
3261  result = delete_dublicates_noteval(result);
3262  return(result);
3263}//proc HomogfacFirstShiftAll
3264
3265//==================================================
3266//factorization of the first Shift Algebra
3267proc facFirstShift(poly h)
3268"USAGE: facFirstShift(h); h a polynomial in the first shift algebra
3269RETURN: list
3270PURPOSE: compute all factorizations of a polynomial in the first shift algebra
3271THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle
3272ASSUME: basering is the first shift algebra
3273NOTE: Every entry of the output list is a list with factors for one possible factorization.
3274EXAMPLE: example facFirstShift; shows examples
3275SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl
3276"{//facFirstShift
3277  int p = printlevel - voice + 2;
3278  int i;
3279  string dbprintWhitespace = "";
3280  for (i = 1; i<=voice;i++)
3281  {dbprintWhitespace = dbprintWhitespace + " ";}
3282  dbprint(p,dbprintWhitespace +" Checking the input.");
3283  if(nvars(basering)!=2)
3284  {//Our basering is the Shift algebra, but not the first
3285    ERROR("Basering is not the first shift algebra");
3286    return(list());
3287  }//Our basering is the Shift algebra, but not the first
3288  def r = basering;
3289  setring r;
3290  list LR = ringlist(r);
3291  number @n = leadcoef(LR[5][1,2]);
3292  poly @p = LR[6][1,2];
3293  if  ( @n!=number(1) )
3294  {
3295    ERROR("Basering is not the first shift algebra");
3296    return(list());
3297  }
3298  dbprint(p,dbprintWhitespace +" Done");
3299  list result = list();
3300  int j; int k; int l; //counter
3301  // create a ring with the ordering which makes shift algebra
3302  // graded
3303  // def r = basering; // done before
3304  ring tempRing = LR[1][1],(x,s),(a(0,1),Dp);
3305  def tempRingnc = nc_algebra(1,s);
3306  setring r;
3307  // information on relations
3308  if (@p == -var(1)) // reverted shift algebra
3309  {
3310    dbprint(p,dbprintWhitespace +" Reverted shift algebra. Swaping variables in Ringlist");
3311    setring(tempRingnc);
3312    map transf = r, var(2), var(1);
3313    setring(r);
3314    map transfback = tempRingnc, var(2),var(1);
3315    //    result = transfback(resulttemp);
3316  }
3317  else
3318  {
3319    if ( @p == var(2)) // usual shift algebra
3320    {
3321      setring(tempRingnc);
3322      map transf = r, var(1), var(2);
3323      //    result = facshift(h);
3324      setring(r);
3325      map transfback = tempRingnc, var(1),var(2);
3326    }
3327    else
3328    {
3329      ERROR("Basering is not the first shift algebra");
3330      return(list());
3331    }
3332  }
3333  // main calls
3334  setring(tempRingnc);
3335  dbprint(p,dbprintWhitespace +" Factorize the given polynomial with the subroutine sFacShift");
3336  list resulttemp = sFacShift(transf(h));
3337  dbprint(p,dbprintWhitespace +" Successful");
3338  setring(r);
3339  result = transfback(resulttemp);
3340  return( delete_dublicates_noteval(result) );
3341}//facFirstShift
3342example
3343{
3344  "EXAMPLE:";echo=2;
3345  ring R = 0,(x,s),dp;
3346  def r = nc_algebra(1,s);
3347  setring(r);
3348  poly h = (s^2*x+x)*s;
3349  facFirstShift(h);
3350}
3351
3352static proc sFacShift(poly h)
3353"
3354USAGE: A static procedure to factorize a polynomial in the first Shift algebra, where all the
3355       validity checks were made in advance.
3356INPUT:  A polynomial h in the first Shift Algebra.
3357OUTPUT: A list of different factorizations of h, where the factors are irreducible
3358ASSUMPTIONS:
3359 - The basering is the first Shift algebra and has n as first, and s as second variable, i.e. we
3360   have var(2)*var(1) = var(1)*var(2)+1
3361THEORY: If the given polynomial h is [0,1]-homogeneous, the routines for homogeneous factorizations
3362 are called. Otherwise we map the polynomial into the first Weyl algebra (the first shift
3363 algebra is a subring of the first Weyl algebra), and use facFirstWeyl to factorize it. Later
3364 we map the factors back, if possible.
3365"
3366{//proc sFacShift
3367  int p = printlevel - voice + 2;
3368  int i; int j ;
3369  string dbprintWhitespace = "";
3370  number commonCoefficient = content(h);
3371  for (i = 1; i<=voice;i++)
3372  {dbprintWhitespace = dbprintWhitespace + " ";}
3373  //Checking if given polynomial is homogeneous
3374  if(homogwithorder(h,intvec(0,1)))
3375  {//The given polynomial is [0,1]-homogeneous
3376    dbprint(p,dbprintWhitespace+"The polynomial is [0,1]-homogeneous. Returning the
3377homogeneous factorization");
3378    return(homogfacFirstShift_all(h));
3379  }//The given polynomial is [0,1]-homogeneous
3380
3381  //---------- Start of interesting part ----------
3382
3383  dbprint(p,dbprintWhitespace+"Mapping the polynomial h into the first Weyl algebra.");
3384  poly temph = h/commonCoefficient;
3385  def ourBaseRing = basering;
3386  ring tempWeylAlgebraComm = 0,(x,d),dp;
3387  def tempWeylAlgebra = nc_algebra(1,1);
3388  setring(tempWeylAlgebra);
3389  map shiftMap = ourBaseRing, x*d, d;
3390  poly h = shiftMap(temph);
3391  dbprint(p,dbprintWhitespace+"Successful! The polynomial in the Weyl algebra is "+string(h));
3392  dbprint(p,dbprintWhitespace+"Factorizing the polynomial in the first Weyl algebra");
3393  list factorizationInWeyl = facFirstWeyl(h);
3394  dbprint(p,dbprintWhitespace+"Successful! The factorization is given by:");
3395  dbprint(p,factorizationInWeyl);
3396  list validCombinations;
3397
3398  dbprint(p,dbprintWhitespace+"Now we will map this back to the shift algebra and filter
3399valid results");
3400  //-Now we map the results back to the shift algebra. But first, we need to combine them properly.
3401  for (i = 1; i<=size(factorizationInWeyl); i++)
3402  {//Deleting the first Coefficient factor
3403    factorizationInWeyl[i] = delete(factorizationInWeyl[i],1);
3404    validCombinations = validCombinations + combineNonnegative(factorizationInWeyl[i]);
3405  }//Deleting the first Coefficient factor
3406  if (size(validCombinations) == 0)
3407  {//There are no valid combinations, therefore we can directly say, that h is irreducible
3408    setring(ourBaseRing);
3409    return(list(list(commonCoefficient, h/commonCoefficient)));
3410  }//There are no valid combinations, therefore we can directly say, that h is irreducible
3411  validCombinations = delete_dublicates_noteval(validCombinations);
3412  setring(ourBaseRing);
3413  map backFromWeyl = tempWeylAlgebra, var(1),var(2);
3414  list validCombinations = backFromWeyl(validCombinations);
3415  for (i = 1; i<=size(validCombinations); i++)
3416  {
3417    for (j = 1; j<=size(validCombinations[i]);j++)
3418    {
3419      setring(tempWeylAlgebra);
3420      fromWeylToShiftPoly(validCombinations[i][j],ourBaseRing);
3421      validCombinations[i][j] = result;
3422      kill result;
3423      kill tempResult;
3424      kill zeroPoly;
3425      kill fromWeyl;
3426    }
3427  }
3428  for (i = 1; i<=size(validCombinations); i++)
3429  {//Adding the common factor in the first position of the list
3430    validCombinations[i] = insert(validCombinations[i],commonCoefficient);
3431  }//Adding the common factor in the first position of the list
3432  dbprint(dbprintWhitespace+"Done.");
3433  //mapping
3434  return(validCombinations);
3435}//proc sFacShift
3436
3437static proc combineNonnegative(list l)
3438"
3439USAGE: In sFacShift, when we want to map back the results of the factorization of the polynomial in
3440       the first Weyl algebra to the shift algebra. We need to recombine the factors such that
3441       we can map it back to the shift algebra without any problems.
3442INPUT: A list l containing one factorization of a polynomial in the first Weyl algebra. For example
3443       for the polynomial (1+x)*(1+x+d) we would have the list [1,x+1,x+d+1].
3444OUTPUT:If we can map every factor without a problem back to the shift algebra (i.e. if the smallest
3445 homogeneous summand of every factor is of nonnegative degree), a list containing the same
3446 list as given in the input is returned.
3447 If otherwise some factors cause problems, we consider every possible combination (i.e.
3448 products of the factors) and extract those where all factors have a smallest homogeneous
3449 summand of nonnegative degree.
3450ASSUMPTIONS:
3451 - Weyl algebra is given, and we have var(2)*var(1)=var(1)*var(2) +1
3452"
3453{//combineNonnegative
3454  int p = printlevel - voice + 2;
3455  int i;
3456  string dbprintWhitespace = "";
3457  for (i = 1; i<=voice;i++)
3458  {dbprintWhitespace = dbprintWhitespace + " ";}
3459  //First the easy case: all of the factors fulfill the condition of mapping to shift:
3460  dbprint(p,dbprintWhitespace+"Checking, if the given factors
3461can already be mapped without a problem.");
3462  int isValid = 1;
3463  for (i = 1; i<=size(l);i++)
3464  {//Checking for every entry if the condition is fulfilled.
3465    if (deg(l[i],intvec(1,-1))>0)
3466    {//Found one, where it is not fulfilled
3467      isValid = 0;
3468      break;
3469    }//Found one, where it is not fulfilled
3470  }//Checking for every entry if the condition is fulfilled.
3471  dbprint(p,dbprintWhitespace+"Done.");
3472  if (isValid)
3473  {//We can map every factor to the shift algebra and do not need to combine anything
3474    dbprint(p,dbprintWhitespace+"They can be mapped. Therefore we return them directly.");
3475    return(list(l));
3476  }//We can map every factor to the shift algebra and do not need to combine anything
3477  dbprint(p,dbprintWhitespace+"They cannot be mapped. Looking for valid combinations.");
3478  //Starting with the case, where l only consists of 1 or two elements.
3479  if(size(l)<=2)
3480  {//The case where we won't call the function a second time
3481    if (deg(product(l),intvec(1,-1))>0)
3482    {//No way of a valid combination
3483      return(list());
3484    }//No way of a valid combination
3485    else
3486    {//The product is the only possible and valid combination
3487      return(list(list(product(l))));
3488    }//The product is the only possible and valid combination
3489  }//The case where we won't call the function a second time
3490  //---------- Easy pre-stuff done. now we combine the factors.----------
3491  int pos;
3492  int j; int k;
3493  dbprint(p,dbprintWhitespace+"Making combinations of two.");
3494  list combinationsOfTwo = combinekfinlf(l,2);
3495  dbprint(p,dbprintWhitespace+"Done. Now checking, if there are valid ones in between.");
3496  list result;
3497  list validLHS;
3498  list validRHS;
3499  for (i = 1; i<=size(combinationsOfTwo); i++)
3500  {//go through all combinations and detect the valid ones
3501    if(deg(combinationsOfTwo[i][1],intvec(1,-1))>0 or deg(combinationsOfTwo[i][2],intvec(1,-1))>0)
3502    {//No chance, so no further treatment needed
3503      i++;
3504      continue;
3505    }//No chance, so no further treatment needed
3506    for (pos = 1; pos<=size(l);pos++)
3507    {//find the position where the combination splits
3508      if (product(l[1..pos]) == combinationsOfTwo[i][1])
3509      {//Found the position
3510        break;
3511      }//Found the position
3512    }//find the position where the combination splits
3513    dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " +
3514            string(list(l[1..pos])));
3515    validLHS = combineNonnegative(list(l[1..pos]));
3516    dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " +
3517            string(list(l[pos+1..size(l)])));
3518    validRHS = combineNonnegative(list(l[pos+1..size(l)]));
3519    for (j = 1; j<=size(validLHS); j++)
3520    {//Combining the left hand side valid combnations...
3521      for (k = 1; k<=size(validRHS); k++)
3522      {//... with the right hand side valid combinations
3523        result = insert(result, validLHS[j]+validRHS[k]);
3524      }//... with the right hand side valid combinations
3525    }//Combining the left hand side valid combnations...
3526  }//go through all combinations and detect the valid ones
3527  result = delete_dublicates_noteval(result);
3528  dbprint(p,dbprintWhitespace+"Done.");
3529  return(result);
3530}//combineNonnegative
3531
3532static proc fromWeylToShiftPoly(poly h, sAlgebra)
3533"
3534USAGE: Given a polynomial in the first Weyl algebra, this method returns it -- if possible --
3535       as an element in the first shift algebra, which is given in the method header.
3536INPUT: A polynomial h, and the first shift algebra as a ring
3537OUTPUT: The correct mapping in the shift Algebra
3538ASSUMPTIONS:
3539 - The lowest [-1,1]-homogeneous summand of h is of nonnegative degree
3540 - The shift algebra is given in the way that var(2)*var(1) = (var(1)+1)*var(2)
3541"
3542{//fromWeylToShiftPoly
3543  int p = printlevel - voice + 2;
3544  int i;
3545  string dbprintWhitespace = "";
3546  for (i = 1; i<=voice;i++)
3547  {dbprintWhitespace = dbprintWhitespace + " ";}
3548  if (deg(h,intvec(1,-1))>0)
3549  {//Wrong input polynomial
3550    ERROR("The lowest [-1,1] homogeneous summand of "+string(h)+" is of negative degree.");
3551  }//Wrong input polynomial
3552  def ourHomeBase = basering;
3553  list hDist = homogDistribution(h);
3554  setring(sAlgebra);
3555  poly result = 0;
3556  poly tempResult;
3557  poly zeroPoly;
3558  map fromWeyl = ourHomeBase, var(1), var(2);
3559  setring(ourHomeBase);
3560  poly zeroPoly;
3561  poly tempZeroPoly;
3562  int j; int k;
3563  int derDeg;
3564  for (i = 1; i<=size(hDist);i++)
3565  {
3566    derDeg = hDist[i][1];
3567    setring(sAlgebra);
3568    tempResult = 1;
3569    setring(ourHomeBase);
3570    zeroPoly = lift(d^derDeg, hDist[i][2])[1,1];
3571    for (j = 1; j<=size(zeroPoly); j++)
3572    {
3573      tempZeroPoly = zeroPoly[j];
3574      setring(sAlgebra);
3575      zeroPoly = fromWeyl(tempZeroPoly);
3576      tempResult = tempResult * leadcoef(zeroPoly);
3577      setring(ourHomeBase);
3578      for (k = 1; k<=deg(zeroPoly[j],intvec(0,1));k++)
3579      {
3580        setring(sAlgebra);
3581        tempResult = tempResult*(var(1)-(k-1));
3582        setring(ourHomeBase);
3583      }
3584      setring(sAlgebra);
3585      result = result + tempResult*var(2)^derDeg;
3586      tempResult = 1;
3587      setring(ourHomeBase);
3588    }
3589  }
3590  setring(sAlgebra);
3591  keepring(sAlgebra);
3592}//fromWeylToShiftPoly
3593
3594static proc refineFactList(list L)
3595{
3596  // assume: list L is an output of factorization proc
3597  // doing: remove doubled entries
3598  int s = size(L); int sm;
3599  int i,j,k,cnt;
3600  list M, U, A, B;
3601  A = L;
3602  k = 0;
3603  cnt  = 1;
3604  for (i=1; i<=s; i++)
3605  {
3606    if (size(A[i]) != 0)
3607    {
3608      M = A[i];
3609      //      "probing with"; M; i;
3610      B[cnt] = M; cnt++;
3611      for (j=i+1; j<=s; j++)
3612      {
3613        if ( isEqualList(M,A[j]) )
3614        {
3615          k++;
3616          // U consists of intvecs with equal pairs
3617          U[k] = intvec(i,j);
3618          A[j] = 0;
3619        }
3620      }
3621    }
3622  }
3623  kill A,U,M;
3624  return(B);
3625}
3626example
3627{
3628  "EXAMPLE:";echo=2;
3629  ring R = 0,(x,s),dp;
3630  def r = nc_algebra(1,1);
3631  setring(r);
3632  list l,m;
3633  l = list(1,s2+1,x,s,x+s);
3634  m = l,list(1,s,x,s,x),l;
3635  refineFactList(m);
3636}
3637
3638static proc isEqualList(list L, list M)
3639{
3640  // int boolean: 1=yes, 0 =no : test whether two lists are identical
3641  int s = size(L);
3642  if (size(M)!=s) { return(0); }
3643  int j=1;
3644  while ( (L[j]==M[j]) && (j<s) )
3645  {
3646    j++;
3647  }
3648  if (L[j]==M[j])
3649  {
3650    return(1);
3651  }
3652  return(0);
3653}
3654example
3655{
3656  "EXAMPLE:";echo=2;
3657  ring R = 0,(x,s),dp;
3658  def r = nc_algebra(1,1);
3659  setring(r);
3660  list l,m;
3661  l = list(1,s2+1,x,s,x+s);
3662  m = l;
3663  isEqualList(m,l);
3664}
3665
3666
3667//////////////////////////////////////////////////
3668// Q-WEYL-SECTION
3669//////////////////////////////////////////////////
3670
3671//==================================================
3672//A function to get the i'th triangular number
3673static proc triangNum(int n)
3674{
3675  if (n == 0)
3676  {
3677    return(0);
3678  }
3679  return (n*(n+1) div 2);
3680}
3681
3682//==================================================*
3683//one factorization of a homogeneous polynomial
3684//in the first Q Weyl Algebra
3685proc homogfacFirstQWeyl(poly h)
3686"USAGE: homogfacFirstQWeyl(h); h is a homogeneous polynomial in the
3687 first q-Weyl algebra with respect to the weight vector [-1,1]
3688RETURN: list
3689PURPOSE: Computes a factorization of a homogeneous polynomial h with
3690  respect to the weight vector [-1,1] in the first q-Weyl algebra
3691THEORY: @code{homogfacFirstQWeyl} returns a list with a factorization of the given,
3692 [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with
3693 k positive, the last k entries in the output list are the second
3694 variable. If k is positive, the last k entries will be x. The other
3695 entries will be irreducible polynomials of degree zero or 1 resp. -1.
3696SEE ALSO: homogfacFirstQWeyl_all
3697"{//proc homogfacFirstQWeyl
3698  int p = printlevel-voice+2;//for dbprint
3699  def r = basering;
3700  poly hath;
3701  int i; int j;
3702  string dbprintWhitespace = "";
3703  for (i = 1; i<=voice;i++)
3704  {dbprintWhitespace = dbprintWhitespace + " ";}
3705  intvec ivm11 = intvec(-1,1);
3706  if (!homogwithorder(h,ivm11))
3707  {//The given polynomial is not homogeneous
3708    ERROR("Given polynomial was not [-1,1]-homogeneous");
3709    return(list());
3710  }//The given polynomial is not homogeneous
3711  if (h==0)
3712  {
3713    return(list(0));
3714  }
3715  list result;
3716  int m = deg(h,ivm11);
3717  dbprint(p,dbprintWhitespace+" Splitting the polynomial in A_0 and A_k-Part");
3718  if (m!=0)
3719  {//The degree is not zero
3720    if (m <0)
3721    {//There are more x than y
3722      hath = lift(var(1)^(-m),h)[1,1];
3723      for (i = 1; i<=-m; i++)
3724      {
3725        result = result + list(var(1));
3726      }
3727    }//There are more x than y
3728    else
3729    {//There are more y than x
3730      hath = lift(var(2)^m,h)[1,1];
3731      for (i = 1; i<=m;i++)
3732      {
3733        result = result + list(var(2));
3734      }
3735    }//There are more y than x
3736  }//The degree is not zero
3737  else
3738  {//The degree is zero
3739    hath = h;
3740  }//The degree is zero
3741  dbprint(p," Done");
3742  //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1)
3743  list mons;
3744  dbprint(p," Putting the monomials in the A_0-part in a list.");
3745  for(i = 1; i<=size(hath);i++)
3746  {//Putting the monomials in a list
3747    mons = mons+list(hath[i]);
3748  }//Putting the monomials in a list
3749  dbprint(p," Done");
3750  dbprint(p," Mapping this monomials to K(q)[theta]");
3751  //Now, map to the commutative ring with theta:
3752  list tempRingList = ringlist(r);
3753  tempRingList[2] = insert(tempRingList[2],"theta",2); //New variable theta = x*d
3754  tempRingList = delete(tempRingList,5);
3755  tempRingList = delete(tempRingList,5); //The ring should now be commutative
3756  def tempRing = ring(tempRingList);
3757  setring tempRing;
3758  map thetamap = r,var(1),var(2);
3759  list mons = thetamap(mons);
3760  poly entry;
3761  poly tempSummand;
3762  for (i = 1; i<=size(mons);i++)
3763  {//transforming the monomials as monomials in theta
3764    entry = 1;//leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
3765    for (j = 0; j<leadexp(mons[i])[2];j++)
3766    {
3767      tempSummand = (par(1)^j-1)/(par(1)-1);
3768      entry = entry * theta-tempSummand*entry;
3769    }
3770    //entry;
3771    //leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
3772    mons[i] = entry*leadcoef(mons[i]) * par(1)^(-triangNum(leadexp(mons[i])[2]-1));
3773  }//transforming the monomials as monomials in theta
3774  dbprint(p," Done");
3775  dbprint(p," Factorize the A_0-Part in K[theta]");
3776  list azeroresult = factorize(sum(mons));
3777  dbprint(p," Successful");
3778  list azeroresult_return_form;
3779  for (i = 1; i<=size(azeroresult[1]);i++)
3780  {//rewrite the result of the commutative factorization
3781    for (j = 1; j <= azeroresult[2][i];j++)
3782    {
3783      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
3784    }
3785  }//rewrite the result of the commutative factorization
3786  dbprint(p," Mapping back to A_0.");
3787  setring(r);
3788  map finalmap = tempRing,var(1),var(2),var(1)*var(2);
3789  list tempresult = finalmap(azeroresult_return_form);
3790  dbprint(p,"Successful.");
3791  for (i = 1; i<=size(tempresult);i++)
3792  {//factorizations of theta resp. theta +1
3793    if(tempresult[i]==var(1)*var(2))
3794    {
3795      tempresult = insert(tempresult,var(1),i-1);
3796      i++;
3797      tempresult[i]=var(2);
3798    }
3799    if(tempresult[i]==var(2)*var(1))
3800    {
3801      tempresult = insert(tempresult,var(2),i-1);
3802      i++;
3803      tempresult[i]=var(1);
3804    }
3805  }//factorizations of theta resp. theta +1
3806  result = tempresult+result;
3807  //Correction of the result in the special q-Case:
3808  for (j = 2 ; j<= size(result);j++)
3809  {//Div the whole Term by the leading coefficient and multiply it to the first entry in result[i]
3810    result[1] = result[1] * leadcoef(result[j]);
3811    result[j] = 1/leadcoef(result[j]) * result[j];
3812  }//Div the whole Term by the leading coefficient and multiply it to the first entry in result[i]
3813  return(result);
3814}//proc homogfacFirstQWeyl
3815example
3816{
3817  "EXAMPLE:";echo=2;
3818  ring R = (0,q),(x,d),dp;
3819  def r = nc_algebra (q,1);
3820  setring(r);
3821  poly h = q^25*x^10*d^10+q^16*(q^4+q^3+q^2+q+1)^2*x^9*d^9+
3822    q^9*(q^13+3*q^12+7*q^11+13*q^10+20*q^9+26*q^8+30*q^7+
3823    31*q^6+26*q^5+20*q^4+13*q^3+7*q^2+3*q+1)*x^8*d^8+
3824    q^4*(q^9+2*q^8+4*q^7+6*q^6+7*q^5+8*q^4+6*q^3+
3825     4*q^2+2q+1)*(q^4+q^3+q^2+q+1)*(q^2+q+1)*x^7*d^7+
3826    q*(q^2+q+1)*(q^5+2*q^4+2*q^3+3*q^2+2*q+1)*(q^4+q^3+q^2+q+1)*(q^2+1)*(q+1)*x^6*d^6+
3827    (q^10+5*q^9+12*q^8+21*q^7+29*q^6+33*q^5+31*q^4+24*q^3+15*q^2+7*q+12)*x^5*d^5+
3828    6*x^3*d^3+24;
3829  homogfacFirstQWeyl(h);
3830}
3831
3832//==================================================
3833//Computes all possible homogeneous factorizations for an element in the first Q-Weyl Algebra
3834proc homogfacFirstQWeyl_all(poly h)
3835"USAGE: homogfacFirstQWeyl_all(h); h is a homogeneous polynomial in the first q-Weyl algebra
3836 with respect to the weight vector [-1,1]
3837RETURN: list
3838PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect
3839  to the weight vector [-1,1] in the first q-Weyl algebra
3840THEORY: @code{homogfacFirstQWeyl} returns a list with all factorization of the given,
3841 homogeneous polynomial. It uses the output of homogfacFirstQWeyl and permutes
3842 its entries with respect to the commutation rule. Furthermore, if a
3843 factor of degree zero is irreducible in K[  heta], but reducible in
3844 the first q-Weyl algebra, the permutations of this element with the other
3845 entries will also be computed.
3846SEE ALSO: homogfacFirstQWeyl
3847"{//proc HomogfacFirstQWeylAll
3848  int p=printlevel-voice+2;//for dbprint
3849  intvec iv11= intvec(1,1);
3850  if (deg(h,iv11) <= 0 )
3851  {//h is a constant
3852    dbprint(p,"Given polynomial was not homogeneous");
3853    return(list(list(h)));
3854  }//h is a constant
3855  def r = basering;
3856  list one_hom_fac; //stands for one homogeneous factorization
3857  int i; int j; int k;
3858  intvec ivm11 = intvec(-1,1);
3859  dbprint(p," Calculate one homogeneous factorization using homogfacFirstQWeyl");
3860  //Compute again a homogeneous factorization
3861  one_hom_fac = homogfacFirstQWeyl(h);
3862  dbprint(p,"Successful");
3863  if (size(one_hom_fac) == 0)
3864  {//there is no homogeneous factorization or the polynomial was not homogeneous
3865    return(list());
3866  }//there is no homogeneous factorization or the polynomial was not homogeneous
3867  //divide list in A0-Part and a list of x's resp. y's
3868  list list_not_azero = list();
3869  list list_azero;
3870  list k_factor;
3871  int is_list_not_azero_empty = 1;
3872  int is_list_azero_empty = 1;
3873  k_factor = list(one_hom_fac[1]);
3874  if (absValue(deg(h,ivm11))<size(one_hom_fac)-1)
3875  {//There is a nontrivial A_0-part
3876    list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,ivm11)))];
3877    is_list_azero_empty = 0;
3878  }//There is a nontrivial A_0 part
3879  dbprint(p," Combine x,y to xy in the factorization again.");
3880  for (i = 1; i<=size(list_azero)-1;i++)
3881  {//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
3882    if (list_azero[i] == var(1))
3883    {
3884      if (list_azero[i+1]==var(2))
3885      {
3886        list_azero[i] = var(1)*var(2);
3887        list_azero = delete(list_azero,i+1);
3888      }
3889    }
3890    if (list_azero[i] == var(2))
3891    {
3892      if (list_azero[i+1]==var(1))
3893      {
3894        list_azero[i] = var(2)*var(1);
3895        list_azero = delete(list_azero,i+1);
3896      }
3897    }
3898  }//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
3899  dbprint(p," Done");
3900  if(deg(h,ivm11)!=0)
3901  {//list_not_azero is not empty
3902    list_not_azero =
3903      one_hom_fac[(size(one_hom_fac)-absValue(deg(h,ivm11))+1)..size(one_hom_fac)];
3904    is_list_not_azero_empty = 0;
3905  }//list_not_azero is not empty
3906  //Map list_azero in K[theta]
3907  dbprint(p," Map list_azero to K[theta]");
3908  //Now, map to the commutative ring with theta:
3909  list tempRingList = ringlist(r);
3910  tempRingList[2] = insert(tempRingList[2],"theta",2); //New variable theta = x*d
3911  tempRingList = delete(tempRingList,5);
3912  tempRingList = delete(tempRingList,5); //The ring should now be commutative
3913  def tempRing = ring(tempRingList);
3914  setring(tempRing);
3915  poly entry;
3916  map thetamap = r,var(1),var(2);
3917  if(!is_list_not_azero_empty)
3918  {//Mapping in Singular is only possible, if the list before
3919    //contained at least one element of the other ring
3920    list list_not_azero = thetamap(list_not_azero);
3921  }//Mapping in Singular is only possible, if the list before
3922  //contained at least one element of the other ring
3923  if(!is_list_azero_empty)
3924  {//Mapping in Singular is only possible, if the list before
3925    //contained at least one element of the other ring
3926    list list_azero= thetamap(list_azero);
3927  }//Mapping in Singular is only possible, if the list before
3928  //contained at least one element of the other ring
3929  list k_factor = thetamap(k_factor);
3930  list tempmons;
3931  dbprint(p," Done");
3932  for(i = 1; i<=size(list_azero);i++)
3933  {//rewrite the polynomials in A1 as polynomials in K[theta]
3934    tempmons = list();
3935    for (j = 1; j<=size(list_azero[i]);j++)
3936    {
3937      tempmons = tempmons + list(list_azero[i][j]);
3938    }
3939    for (j = 1 ; j<=size(tempmons);j++)
3940    {
3941      //entry = leadcoef(tempmons[j]);
3942      entry = leadcoef(tempmons[j]) * par(1)^(-triangNum(leadexp(tempmons[j])[2]-1));
3943      for (k = 0; k < leadexp(tempmons[j])[2];k++)
3944      {
3945        entry = entry*(theta-(par(1)^k-1)/(par(1)-1));
3946      }
3947      tempmons[j] = entry;
3948    }
3949    list_azero[i] = sum(tempmons);
3950  }//rewrite the polynomials in A1 as polynomials in K[theta]
3951  //Compute all permutations of the A0-part
3952  dbprint(p," Compute all permutations of the A_0-part with the first resp. the snd. variable");
3953  list result;
3954  int shift_sign;
3955  int shift;
3956  poly shiftvar;
3957  if (size(list_not_azero)!=0)
3958  {//Compute all possibilities to permute the x's resp. the y's in the list
3959    if (list_not_azero[1] == var(1))
3960    {//h had a negative weighted degree
3961      shift_sign = 1;
3962      shiftvar = var(1);
3963    }//h had a negative weighted degree
3964    else
3965    {//h had a positive weighted degree
3966      shift_sign = -1;
3967      shiftvar = var(2);
3968    }//h had a positive weighted degree
3969    result = permpp(list_azero + list_not_azero);
3970    for (i = 1; i<= size(result); i++)
3971    {//adjust the a_0-parts
3972      shift = 0;
3973      for (j=1; j<=size(result[i]);j++)
3974      {
3975        if (result[i][j]==shiftvar)
3976        {
3977          shift = shift + shift_sign;
3978        }
3979        else
3980        {
3981          if (shift < 0)
3982          {//We have two distict formulas for x and y. In this case use formula for y
3983            if (shift == -1)
3984            {
3985              result[i][j] = subst(result[i][j],theta,1/par(1)*(theta - 1));
3986            }
3987            else
3988            {
3989              result[i][j] =
3990                subst(result[i][j],
3991                      theta,
3992                      1/par(1)*((theta - 1)/par(1)^(absValue(shift)-1)
3993                                - (par(1)^(shift +2)-par(1))/(1-par(1))));
3994            }
3995          }//We have two distict formulas for x and y. In this case use formula for y
3996          if (shift > 0)
3997          {//We have two distict formulas for x and y. In this case use formula for x
3998            if (shift == 1)
3999            {
4000              result[i][j] = subst(result[i][j],theta,par(1)*theta + 1);
4001            }
4002            else
4003            {
4004              result[i][j] =
4005                subst(result[i][j],
4006                      theta,par(1)^shift*theta+(par(1)^shift-1)/(par(1)-1));
4007            }
4008          }//We have two distict formulas for x and y. In this case use formula for x
4009        }
4010      }
4011    }//adjust the a_0-parts
4012  }//Compute all possibilities to permute the x's resp. the y's in the list
4013  else
4014  {//The result is just all the permutations of the a_0-part
4015    result = permpp(list_azero);
4016  }//The result is just all the permutations of the a_0 part
4017  if (size(result)==0)
4018  {
4019    return(result);
4020  }
4021  dbprint(p," Done");
4022  dbprint(p," Searching for theta resp. theta + 1 in the list and factorize them");
4023  //Now we are going deeper and search for theta resp. theta + 1, substitute
4024  //them by xy resp. yx and go on permuting
4025  int found_theta;
4026  int thetapos;
4027  list leftpart;
4028  list rightpart;
4029  list lparts;
4030  list rparts;
4031  list tempadd;
4032  for (i = 1; i<=size(result) ; i++)
4033  {//checking every entry of result for theta or theta +1
4034    found_theta = 0;
4035    for(j=1;j<=size(result[i]);j++)
4036    {
4037      if (result[i][j]==theta)
4038      {//the jth entry is theta and can be written as x*y
4039        thetapos = j;
4040        result[i]= insert(result[i],var(1),j-1);
4041        j++;
4042        result[i][j] = var(2);
4043        found_theta = 1;
4044        break;
4045      }//the jth entry is theta and can be written as x*y
4046      if(result[i][j] == par(1)*theta +1)
4047      {
4048        thetapos = j;
4049        result[i] = insert(result[i],var(2),j-1);
4050        j++;
4051        result[i][j] = var(1);
4052        found_theta = 1;
4053        break;
4054      }
4055    }
4056    if (found_theta)
4057    {//One entry was theta resp. theta +1
4058      leftpart = result[i];
4059      leftpart = leftpart[1..thetapos];
4060      rightpart = result[i];
4061      rightpart = rightpart[(thetapos+1)..size(rightpart)];
4062      lparts = list(leftpart);
4063      rparts = list(rightpart);
4064      //first deal with the left part
4065      if (leftpart[thetapos] == var(1))
4066      {
4067        shift_sign = 1;
4068        shiftvar = var(1);
4069      }
4070      else
4071      {
4072        shift_sign = -1;
4073        shiftvar = var(2);
4074      }
4075      for (j = size(leftpart); j>1;j--)
4076      {//drip x resp. y
4077        if (leftpart[j-1]==shiftvar)
4078        {//commutative
4079          j--;
4080          continue;
4081        }//commutative
4082        if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
4083        {//stop here
4084          break;
4085        }//stop here
4086        //Here, we can only have a a0- part
4087        if (shift_sign<0)
4088        {
4089          leftpart[j] = subst(leftpart[j-1],theta, 1/par(1)*(theta +shift_sign));
4090        }
4091        if (shift_sign>0)
4092        {
4093          leftpart[j] = subst(leftpart[j-1],theta, par(1)*theta + shift_sign);
4094        }
4095        leftpart[j-1] = shiftvar;
4096        lparts = lparts + list(leftpart);
4097      }//drip x resp. y
4098      //and now deal with the right part
4099      if (rightpart[1] == var(1))
4100      {
4101        shift_sign = 1;
4102        shiftvar = var(1);
4103      }
4104      else
4105      {
4106        shift_sign = -1;
4107        shiftvar = var(2);
4108      }
4109      for (j = 1 ; j < size(rightpart); j++)
4110      {
4111        if (rightpart[j+1] == shiftvar)
4112        {
4113          j++;
4114          continue;
4115        }
4116        if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
4117        {
4118          break;
4119        }
4120        if (shift_sign<0)
4121        {
4122          rightpart[j] = subst(rightpart[j+1], theta, par(1)*theta - shift_sign);
4123        }
4124        if (shift_sign>0)
4125        {
4126          rightpart[j] = subst(rightpart[j+1], theta, 1/par(1)*(theta - shift_sign));
4127        }
4128        rightpart[j+1] = shiftvar;
4129        rparts = rparts + list(rightpart);
4130      }
4131      //And now, we put all possibilities together
4132      tempadd = list();
4133      for (j = 1; j<=size(lparts); j++)
4134      {
4135        for (k = 1; k<=size(rparts);k++)
4136        {
4137          tempadd = tempadd + list(lparts[j]+rparts[k]);
4138        }
4139      }
4140      tempadd = delete(tempadd,1); // The first entry is already in the list
4141      result = result + tempadd;
4142      continue; //We can may be not be done already with the ith entry
4143    }//One entry was theta resp. theta +1
4144  }//checking every entry of result for theta or theta +1
4145  dbprint(p," Done");
4146  //map back to the basering
4147  dbprint(p," Mapping back everything to the basering");
4148  setring(r);
4149  map finalmap = tempRing, var(1), var(2),var(1)*var(2);
4150  list result = finalmap(result);
4151  for (i=1; i<=size(result);i++)
4152  {//adding the K factor
4153    result[i] = k_factor + result[i];
4154  }//adding the k-factor
4155  dbprint(p," Done");
4156  dbprint(p," Delete double entries in the list.");
4157  result = delete_dublicates_noteval(result);
4158  dbprint(p," Done");
4159  return(result);
4160}//proc HomogfacFirstQWeylAll
4161example
4162{
4163  "EXAMPLE:";echo=2;
4164  ring R = (0,q),(x,d),dp;
4165  def r = nc_algebra (q,1);
4166  setring(r);
4167  poly h = q^25*x^10*d^10+q^16*(q^4+q^3+q^2+q+1)^2*x^9*d^9+
4168    q^9*(q^13+3*q^12+7*q^11+13*q^10+20*q^9+26*q^8+30*q^7+
4169    31*q^6+26*q^5+20*q^4+13*q^3+7*q^2+3*q+1)*x^8*d^8+
4170    q^4*(q^9+2*q^8+4*q^7+6*q^6+7*q^5+8*q^4+6*q^3+
4171    4*q^2+2q+1)*(q^4+q^3+q^2+q+1)*(q^2+q+1)*x^7*d^7+
4172    q*(q^2+q+1)*(q^5+2*q^4+2*q^3+3*q^2+2*q+1)*(q^4+q^3+q^2+q+1)*(q^2+1)*(q+1)*x^6*d^6+
4173    (q^10+5*q^9+12*q^8+21*q^7+29*q^6+33*q^5+31*q^4+24*q^3+15*q^2+7*q+12)*x^5*d^5+
4174    6*x^3*d^3+24;
4175  homogfacFirstQWeyl_all(h);
4176}
4177
4178//TODO: FirstQWeyl check the parameters...
4179
4180//==================================================
4181// EASY EXAMPLES FOR WEYL ALGEBRA
4182//==================================================
4183/*
4184  Easy and fast example polynomials where one can find factorizations: K<x,d |dx=xd+1>
4185  (x^2+d)*(x^2+d);
4186  (x^2+x)*(x^2+d);
4187  (x^3+x+1)*(x^4+d*x+2);
4188  (x^2*d+d)*(d+x*d);
4189  d^3+x*d^3+2*d^2+2*(x+1)*d^2+d+(x+2)*d; //Example 5 Grigoriev-Schwarz.
4190  (d+1)*(d+1)*(d+x*d); //Landau Example projected to the first dimension.
4191*/
4192
4193//==================================================
4194//Some Bugs(fixed)/hard examples from Martin Lee:
4195//==================================================
4196// ex1, ex2
4197/*
4198ring s = 0,(x,d),Ws(-1,1);
4199def S = nc_algebra(1,1); setring S;
4200poly a = 10x5d4+26x4d5+47x5d2-97x4d3; //Not so hard any more... Done in around 4 minutes
4201def l= facFirstWeyl (a); l;
4202kill l;
4203poly b = -5328x8d5-5328x7d6+720x9d2+720x8d3-16976x7d4-38880x6d5
4204-5184x7d3-5184x6d4-3774x5d5+2080x8d+5760x7d2-6144x6d3-59616x5d4
4205+3108x3d6-4098x6d2-25704x5d3-21186x4d4+8640x6d-17916x4d3+22680x2d5
4206+2040x5d-4848x4d2-9792x3d3+3024x2d4-10704x3d2-3519x2d3+34776xd4
4207+12096xd3+2898d4-5040x2d+8064d3+6048d2; //Still very hard... But it seems to be only because of the
4208//combinatorial explosion
4209def l= facFirstWeyl (b); l;
4210
4211// ex3: there was difference in answers => fixed
4212LIB "ncfactor.lib";
4213ring r = 0,(x,y,z),dp;
4214matrix D[3][3]; D[1,3]=-1;
4215def R = nc_algebra(1,D);
4216setring R;
4217poly g= 7*z4*x+62*z3+26*z;
4218def l1= facSubWeyl (g, x, z);
4219l1;
4220//---- other ring
4221ring s = 0,(x,z),dp;
4222def S = nc_algebra(1,-1); setring S;
4223poly g= 7*z4*x+62*z3+26*z;
4224def l2= facFirstWeyl (g);
4225l2;
4226map F = R,x,0,z;
4227list l1 = F(l1);
4228l1;
4229//---- so the answers look different, check them!
4230testNCfac(l2); // ok
4231testNCfac(l1); // was not ok, but now it's been fixed!!!
4232
4233// selbst D und X so vertauschen dass sie erfuellt ist : ist gemacht
4234
4235*/
4236
4237/*
4238// bug from M Lee
4239LIB "ncfactor.lib";
4240ring s = 0,(z,x),dp;
4241def S = nc_algebra(1,1); setring S;
4242poly f= -60z4x2-54z4-56zx3-59z2x-64;
4243def l= facFirstWeyl (f);
4244l; // before: empty list; after fix: 1 entry, f is irreducible
4245poly g = 75z3x2+92z3+24;
4246def l= facFirstWeyl (g);
4247l; //before: empty list, now: correct
4248*/
4249
4250/* more things from Martin Lee; fixed
4251ring R = 0,(x,s),dp;
4252def r = nc_algebra(1,s);
4253setring(r);
4254poly h = (s2*x+x)*s;
4255h= h* (x+s);
4256def l= facFirstShift(h);
4257l; // contained doubled entries: not anymore, fixed!
4258
4259ring R = 0,(x,s),dp;
4260def r = nc_algebra(1,-1);
4261setring(r);
4262poly h = (s2*x+x)*s;
4263h= h* (x+s);
4264def l= facFirstWeyl(h);
4265l; // contained doubled entries: not anymore, fixed!
4266
4267*/
4268
4269//======================================================================
4270//Examples from TestSuite that are terminating in a reasonable time.
4271//======================================================================
4272
4273//Counter example for old Algorithm, but now working:
4274/*
4275ring R = 0,(x,d),dp;
4276def r = nc_algebra(1,1);
4277setring(r);
4278LIB "ncfactor.lib";
4279poly h = (1+x^2*d)^4;
4280list lsng = facFirstWeyl(h);
4281print(lsng);
4282*/
4283
4284//Example 2.7. from Master thesis
4285/*
4286ring R = 0,(x,d),dp;
4287def r = nc_algebra(1,1);
4288setring(r);
4289LIB "ncfactor.lib";
4290poly h = (xdd + xd+1+ (xd+5)*x)*(((x*d)^2+1)*d + xd+3+ (xd+7)*x);
4291list lsng = facFirstWeyl(h);
4292print(lsng);
4293 */
4294
4295//Example with high combinatorial income
4296/*
4297ring R = 0,(x,d),dp;
4298def r = nc_algebra(1,1);
4299setring(r);
4300LIB "ncfactor.lib";
4301poly h = (xdddd + (xd+1)*d*d+ (xd+5)*x*d*d)*(((x*d)^2+1)*d*x*x + (xd+3)*x*x+ (xd+7)*x*x*x);
4302list lsng = facFirstWeyl(h);
4303print(lsng);
4304 */
4305
4306//Once a bug, now working
4307/*
4308ring R = 0,(x,d),dp;
4309def r = nc_algebra(1,1);
4310setring(r);
4311LIB "ncfactor.lib";
4312poly h = (x^2*d^2+x)*(x+1);
4313list lsng = facFirstWeyl(h);
4314print(lsng);
4315*/
4316
4317//Another one of that kind
4318/*
4319ring R = 0,(x,d),dp;
4320def r = nc_algebra(1,1);
4321setring(r);
4322LIB "ncfactor.lib";
4323poly h = (x*d*d + (x*d)^5 +x)*((x*d+1)*d-(x*d-1)^5+x);
4324list lsng = facFirstWeyl(h);
4325print(lsng);
4326 */
4327
4328//Example of Victor for Shift Algebra
4329/*
4330ring s = 0,(n,Sn),dp;
4331def S = nc_algebra(1,Sn); setring S;
4332LIB "ncfactor.lib";
4333list lsng = facFirstShift(n^2*Sn^2+3*n*Sn^2-n^2+2*Sn^2-3*n-2);
4334print(lsng);
4335 */
4336
4337//Interesting example, as there are actually also some complex solutions to it:
4338/*
4339ring R = 0,(x,d),dp;
4340def r = nc_algebra(1,1);
4341setring(r);
4342LIB "/Users/albertheinle/Studium/forschung/ncfactor/versionen/ncfactor.lib";
4343poly h =(x^3+x+1)*(x^4+d*x+2);//Example for finitely many, but more than one solution in between.
4344list lsng = facFirstWeyl(h);
4345print(lsng);
4346 */
4347
4348//Another one of that kind:
4349/*
4350ring R = 0,(x,d),dp;
4351def r = nc_algebra(1,1);
4352setring(r);
4353LIB "ncfactor.lib";
4354poly h =(x^2+d)*(x^2+d);//Example for finitely many, but more than one solution in between.
4355list lsng = facFirstWeyl(h);
4356print(lsng);
4357*/
4358
4359//Example by W. Koepf:
4360/*
4361ring R = 0,(x,d),dp;
4362def r = nc_algebra(1,1);
4363setring(r);
4364LIB "ncfactor.lib";
4365poly h = (x^4-1)*x*d^2+(1+7*x^4)*d+8*x^3;
4366list lsng = facFirstWeyl(h);
4367print(lsng);
4368 */
4369
4370//Shift Example from W. Koepf
4371/*
4372ring R = 0,(n,s),dp;
4373def r = nc_algebra(1,s);
4374setring(r);
4375LIB "ncfactor.lib";
4376poly h = n*(n+1)*s^2-2*n*(n+100)*s+(n+99)*(n+100);
4377list lsng = facFirstShift(h);
4378print(lsng);
4379 */
4380
4381//Tsai Example... Once hard, now easy...
4382/*
4383ring R = 0,(x,d),dp;
4384def r = nc_algebra(1,1);
4385setring(r);
4386LIB "ncfactor.lib";
4387poly h = (x^6+2*x^4-3*x^2)*d^2-(4*x^5-4*x^4-12*x^2-12*x)*d + (6*x^4-12*x^3-6*x^2-24*x-12);
4388list lsng =facFirstWeyl(h);
4389print(lsng);
4390 */
4391
4392//======================================================================
4393// Hard examples not yet calculatable in feasible amount of time
4394//======================================================================
4395
4396//Also a counterexample for REDUCE. Very long Groebner basis computation in between.
4397/*
4398ring R = 0,(x,d),dp;
4399def r = nc_algebra(1,1);
4400setring(r);
4401LIB "ncfactor.lib";
4402poly h = (d^4+x^2+dx+x)*(d^2+x^4+xd+d);
4403list lsng = facFirstWeyl(h);
4404print(lsng);
4405*/
4406
4407//Example from the Mainz-Group
4408/*
4409ring R = 0,(x,d),dp;
4410def r = nc_algebra(1,1);
4411setring(r);
4412poly dop6 = 1/35*x^4*(27-70*x+35*x^2)+ 1/35*x*(32+152*x+100*x^2-59*x^3+210*x^4+105*x^5)*d+
4413(-10368/35-67056/35*x-35512/7*x^2-50328/7*x^3-40240/7*x^4-2400*x^5-400*x^6)*d^2+
4414(-144/35*(x+1)*(1225*x^5+11025*x^4+37485*x^3+61335*x^2+50138*x+16584)-6912/35*(x+2)*
4415(x+1)*(105*x^4+1155*x^3+4456*x^2+7150*x+4212) -27648/35*(x+3)*(x+1)*(35*x^2+350*x+867)*
4416(x+2)^2)*d^3;
4417LIB "ncfactor.lib";
4418printlevel = 5;
4419facFirstWeyl(dop6);
4420$;*/
4421
4422//Another Mainz Example:
4423/*
4424LIB "ncfactor.lib";
4425ring R = 0,(x,d),dp;
4426def r = nc_algebra(1,1);
4427setring(r);
4428poly dopp = 82547*x^4*d^4+60237*x^3*d^3+26772*x^5*d^5+2231*x^6*d^6+x*(1140138*
4429x^2*d^2-55872*x*d-3959658*x^3*d^3-8381805*x^4*d^4-3089576*x^5*d^5-274786*
4430x^6*d^6)+x^2*(-16658622*x*d-83427714*x^2*d^2-19715033*x^3*d^3+78915395*x^4
4431*d^4+35337930*x^5*d^5+3354194*x^6*d^6)+x^3*(-99752472-1164881352*x*d+
44324408536996*x^2*d^2+11774185985*x^3*d^3+5262196786*x^4*d^4+1046030561/2*x^5*
4433d^5-10564451/2*x^6*d^6)+x^4*(-1925782272+21995375398*x*d+123415803356*x^2*
4434d^2+302465300831/2*x^3*d^3+34140803907/2*x^4*d^4-15535653409*x^5*d^5-\
44352277687768*x^6*d^6)+x^5*(71273525520+691398212366*x*d+901772633569*x^2*d^2+
44362281275427069*x^3*d^3+2944352819911/2*x^4*d^4+836872370039/4*x^5*d^5+
44379066399237/4*x^6*d^6)+x^6*(2365174430376+9596715855542*x*d+29459572469704*x^
44382*d^2+92502197003786*x^3*d^3+65712473180525*x^4*d^4+13829360193674*x^5*d^5
4439+3231449477251/4*x^6*d^6)+x^7*(26771079436836+117709870166226*x*d+
4440821686455179082*x^2*d^2+1803972139232179*x^3*d^3+1083654460691481*x^4*d^4+
4441858903621851785/4*x^5*d^5+50096565802957/4*x^6*d^6)+x^8*(179341727601960+
44422144653944040630*x*d+13123246960284302*x^2*d^2+41138357917778169/2*x^3*d^3+
444320605819587976401/2*x^4*d^4+3677396642905423/2*x^5*d^5+402688260229369/4*x^6
4444*d^6)+x^9*(2579190935961288+43587063726809764*x*d+157045086382352387*x^2*d^
44452+172175668477370223*x^3*d^3+138636285385875407/2*x^4*d^4+10707836398626232*
4446x^5*d^5+529435530567584*x^6*d^6)+x^10*(41501953525903392+558336731465626084*
4447x*d+1407267553543222268*x^2*d^2+1153046693323226808*x^3*d^3+
4448372331468563656085*x^4*d^4+48654019090240214*x^5*d^5+2114661191282167*x^6*d
4449^6)+x^11*(364526077273381884+4158060401095928464*x*d+8646807662899324262*x^2*
4450d^2+5914675753405705400*x^3*d^3+1631934058875116005*x^4*d^4+
4451187371894330537204*x^5*d^5+7366806367019734*x^6*d^6)+x^12*(
44521759850321214603648+18265471270535733520*x*d+34201910114871110912*x^2*d^2+
445321265221434709398152*x^3*d^3+5437363546219595036*x^4*d^4+594029113431041060*
4454x^5*d^5+22881659624561644*x^6*d^6)+x^13*(4648382639403200688+
445545699084277107816096*x*d+81049061578449009384*x^2*d^2+48858488665016574368*x
4456^3*d^3+12515362110098721444*x^4*d^4+1412152747420021048*x^5*d^5+
445757196947123984972*x^6*d^6)+x^14*(5459369397960020544+55837825300341621824*x*
4458d+105671876924055409696*x^2*d^2+71551727420848766624*x^3*d^3+
445921094786205096577808*x^4*d^4+2695663190297032192*x^5*d^5+118791751565613264*
4460x^6*d^6)+x^15*(1023333653580043776+47171127937488813824*x*d+
4461157258351906685700352*x^2*d^2+145765192195300531840*x^3*d^3+
446249876215785510342176*x^4*d^4+6647374188802036864*x^5*d^5+287310278455067312*
4463x^6*d^6)+x^16*(11960091747366236160+250326608568269289472*x*d+
4464677587171115580981248*x^2*d^2+538246374825683603456*x^3*d^3+
4465161380433451548754048*x^4*d^4+19149099315354950144*x^5*d^5+
4466746433247985092544*x^6*d^6)+x^17*(42246252365448668160+657220532737851248640*
4467x*d+1531751689216283911680*x^2*d^2+1090829514212206064640*x^3*d^3+
4468299280728709430851840*x^4*d^4+32932767387222323200*x^5*d^5+
44691202281367574179840*x^6*d^6)+x^18*(6239106101942784000+320638742839606579200*
4470x*d+873857213570556364800*x^2*d^2+645649080101933721600*x^3*d^3+
4471177008238160627276800*x^4*d^4+19165088507111475200*x^5*d^5+
4472683600826675660800*x^6*d^6)+x^19*(-60440251454613504000-476055211197689856000
4473*x*d-733497382597635072000*x^2*d^2-386038662982742016000*x^3*d^3-\
447483361486778142976000*x^4*d^4-7524999543181824000*x^5*d^5-232189492987008000*
4475x^6*d^6)+x^20*(1578562930483200000+12628503443865600000*x*d+
447619732036631040000000*x^2*d^2+10523752869888000000*x^3*d^3+
44772302070940288000000*x^4*d^4+210475057397760000*x^5*d^5+6577345543680000*x^6*
4478d^6);
4479printlevel = 3;
4480facFirstWeyl(dopp);
4481*/
4482
4483
4484
4485//Hard Example by Viktor:
4486/*
4487  ring r = 0,(x,d), (dp);
4488def R = nc_algebra(1,1);
4489setring R;
4490LIB "/Users/albertheinle/Studium/forschung/ncfactor/versionen/ncfactor.lib";
4491poly t = x; poly D =d;
4492poly p = 2*t^2*D^8-6*t*D^8+2*t^2*D^7+8*t*D^7+12*D^7-2*t^4*D^6+6*t^3*D^6+12*t*D^6-20*D^6
4493-2*t^4*D^5-8*t^3*D^5-4*t^2*D^5+12*t*D^5-28*D^5-12*t^3*D^4-4*t^2*D^4-4*t*D^4-24*D^4+4*t^4*D^3
4494-12*t^3*D^3+2*t^2*D^3-18*t*D^3+16*D^3+6*t^4*D^2-2*t^3*D^2+2*t^2*D^2-2*t*D^2+44*D^2+2*t^4*D
4495+12*t^3*D+2*t*D+4*t^3-8;
4496list lsng = facFirstWeyl(p);
4497print(lsng);
4498*/
Note: See TracBrowser for help on using the repository browser.