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

spielwiese
Last change on this file since d91c9d was d91c9d, checked in by Hans Schoenemann <hannes@…>, 11 years ago
update by albert.heinle@rwth-aachen.de Conflicts: Singular/LIB/ncfactor.lib
  • Property mode set to 100644
File size: 152.7 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 +" recursively check factors for irreducibility");
1173  list recursivetemp;
1174  int changedSomething;
1175  for(i = 1; i<=size(result);i++)
1176  {//recursively factorize factors
1177    if(size(result[i])>2)
1178    {//Nontrivial factorization
1179      for (j=2;j<=size(result[i]);j++)
1180      {//Factorize every factor
1181        recursivetemp = facFirstWeyl(result[i][j]);
1182        //if(size(recursivetemp)>1)
1183        //{//we have a nontrivial factorization
1184        changedSomething = 0;
1185        for(k=1; k<=size(recursivetemp);k++)
1186        {//insert factorized factors
1187          if(size(recursivetemp[k])>2)
1188          {//nontrivial
1189            changedSomething = 1;
1190            result = insert(result,result[i],i);
1191            for(l = size(recursivetemp[k]);l>=2;l--)
1192            {
1193              result[i+1] = insert(result[i+1],recursivetemp[k][l],j);
1194            }
1195            result[i+1] = delete(result[i+1],j);
1196          }//nontrivial
1197        }//insert factorized factors
1198        if (changedSomething)
1199        {
1200          result = delete(result,i);
1201        }
1202        //}//we have a nontrivial factorization
1203      }//Factorize every factor
1204    }//Nontrivial factorization
1205  }//recursively factorize factors
1206  dbprint(p,dbprintWhitespace +" Done");
1207  if (size(result)==0)
1208  {//only the trivial factorization could be found
1209    result = list(list(1,h));
1210  }//only the trivial factorization could be found
1211  list resultWithInterchanges;
1212  dbprint(p,dbprintWhitespace+ "And the result without interchanges with homogeneous factors is:");
1213  dbprint(p,result);
1214  for (i = 1; i <= size(result) ; i++)
1215  {//applying the interchanges to result
1216    resultWithInterchanges = resultWithInterchanges +
1217                             checkForHomogInhomogInterchangability(result[i],2,size(result[i]));
1218  }//applying the interchanges to result
1219  //now, refine the possible redundant list
1220  return( delete_dublicates_noteval(resultWithInterchanges) );
1221}//proc facFirstWeyl
1222example
1223{
1224  "EXAMPLE:";echo=2;
1225  ring R = 0,(x,y),dp;
1226  def r = nc_algebra(1,1);
1227  setring(r);
1228  poly h = (x^2*y^2+x)*(x+1);
1229  facFirstWeyl(h);
1230}
1231
1232//////////////////////////////////////////////////
1233/////BRANDNEW!!!!////////////////////
1234//////////////////////////////////////////////////
1235
1236static proc checkForHomogInhomogInterchangability(list factors, posLeft, posRight)
1237"
1238INPUT:  A list consisting of factors of a certain polynomial in the first Weyl
1239        algebra, factors, and a position from the left and the right, where the last swap was done.
1240OUTPUT: A list containing lists consisting of factors of a certain polynomial in the first Weyl
1241        algebra.
1242The purpose of this function is to check whether we can interchange certain inhomogeneous factors
1243with homogeneous ones. If it is possible, this function returns a list of lists
1244of possible interchanges.
1245
1246The idea came because of an example, where we need an extra swap in the end, otherwise we would
1247not capture all factorizations. The example was
1248h = x4d7+11x3d6+x2d7+x2d6+x3d4+29x2d5+xd6+8xd5+d6+5x2d3+14xd4+13d4+5xd2+d3+d;
1249
1250ASSUMPTIONS:
1251
1252- All factors are irreducible
1253"
1254{//checkForHomogInhomogInterchangability
1255  int p = printlevel-voice+2;
1256  string dbprintWhitespace = "";
1257  int i; int j;
1258  for (i = 1; i<=voice;i++)
1259  {dbprintWhitespace = dbprintWhitespace + " ";}
1260  if (size(factors) <= 2 || posLeft == posRight - 1)
1261  {//easiest case: There is nothing to swap
1262    return (list(factors));
1263  }//easiest case: There is nothing to swap
1264  list result = list(factors);
1265  list tempResultEntries;
1266  list tempSwaps;
1267  list tempSwapsTempEntry;
1268  list attemptToSwap;
1269  intvec ivm11 = intvec(-1,1);
1270  dbprint(p, dbprintWhitespace+"We try to swap elements in the following list:");
1271  dbprint(p, factors);
1272  for (i = posLeft; i < posRight; i++)
1273  {//checking within the window posLeft <--> posRight, if there are interchanges possible
1274    if (homogwithorder(factors[i],ivm11) && !homogwithorder(factors[i+1],ivm11))
1275    {//position i is homogeneous, position i+1 is not ==> trying to swap
1276      attemptToSwap = extractHomogeneousDivisorsRight(factors[i]*factors[i+1]);
1277      if (size(attemptToSwap[1])>1)
1278      {//Bingo, we were able to swap this one element
1279        dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1));
1280        dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:");
1281        dbprint(p,attemptToSwap);
1282        tempSwapsTempEntry = list();
1283        for (j = size(factors); j >=1; j--)
1284        {//creating a new entry for the resulting list, replacing the swap in factors
1285          if (j==i+1)
1286          {tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][2]);}
1287          else
1288          {
1289            if (j == i)
1290            {tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][1]);}
1291            else
1292            {tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]);}
1293          }
1294        }//creating a new entry for the resulting list, replacing the swap in factors
1295        tempSwaps = insert(tempSwaps,list(list(i+1,posRight),tempSwapsTempEntry));
1296      }//Bingo, we were able to swap this one element
1297    }//position i is homogeneous, position i+1 is not ==> trying to swap
1298    else
1299    {
1300      if(!homogwithorder(factors[i],ivm11) && homogwithorder(factors[i+1],ivm11))
1301      {//position i+1 is homogeneous, position i is not ==> trying to swap
1302        attemptToSwap = extractHomogeneousDivisorsLeft(factors[i]*factors[i+1]);
1303        if (size(attemptToSwap[1])>1)
1304        {//Bingo, we were able to swap this one element
1305          dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1));
1306          dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:");
1307          dbprint(p,attemptToSwap);
1308          tempSwapsTempEntry = list();
1309          for (j = size(factors); j >=1; j--)
1310          {//creating a new entry for the resulting list, replacing the swap in factors
1311            if (j==i+1)
1312            {tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][2]);}
1313            else
1314            {
1315              if (j == i)
1316              {tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][1]);}
1317              else
1318              {tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]);}
1319            }
1320          }//creating a new entry for the resulting list, replacing the swap in factors
1321          tempSwaps = insert(tempSwaps,list(list(posLeft,i),tempSwapsTempEntry));
1322        }//Bingo, we were able to swap this one element
1323      }//position i+1 is homogeneous, position i is not ==> trying to swap
1324    }
1325  }//checking within the window posLeft <--> posRight, if there are interchanges possible
1326
1327  //Now we will recursively call the function for all swapped entries.
1328  dbprint(p,dbprintWhitespace+ "Our list of different factorizations is now:");
1329  dbprint(p,tempSwaps);
1330  for (i = 1; i<=size(tempSwaps);i++)
1331  {//recursive call to all formerly attempted swaps.
1332    tempResultEntries=checkForHomogInhomogInterchangability(tempSwaps[i][2],
1333                                                            tempSwaps[i][1][1],tempSwaps[i][1][2]);
1334    result = result + tempResultEntries;
1335  }//recursive call to all formerly attempted swaps.
1336  result = delete_dublicates_noteval(result);
1337  return(result);
1338}//checkForHomogInhomogInterchangability
1339
1340static proc sfacwa(poly h)
1341"INPUT: A polynomial h in the first Weyl algebra
1342OUTPUT: A list of factorizations, where the factors might still be reducible.
1343ASSUMPTIONS:
1344- Our basering is the first Weyl algebra; the x is the first variable,
1345  the differential operator the second.
1346"
1347{//proc sfacwa
1348  int i; int j; int k;
1349  int p = printlevel-voice+2;
1350  string dbprintWhitespace = "";
1351  number commonCoefficient = content(h);
1352  for (i = 1; i<=voice;i++)
1353  {dbprintWhitespace = dbprintWhitespace + " ";}
1354  dbprint(p,dbprintWhitespace + " Extracting homogeneous left and right factors");
1355  if(homogwithorder(h,intvec(-1,1)))
1356  {//we are already dealing with a -1,1 homogeneous poly
1357    dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous. Start homog.
1358fac. and ret. its result");
1359    return(homogfacFirstWeyl_all(h));
1360  }//we are already dealing with a -1,1 homogeneous poly
1361  list resulttemp = extractHomogeneousDivisors(h/commonCoefficient);
1362  //resulttemp = resulttemp + list(list(h/commonCoefficient));
1363  list inhomogeneousFactorsToFactorize;
1364  int isAlreadyInInhomogList;
1365  dbprint(p,dbprintWhitespace +" Done");
1366  dbprint(p,dbprintWhitespace +" Making Set of inhomogeneous polynomials we have to factorize.");
1367  for (i = 1; i<=size(resulttemp); i++)
1368  {//Going through all different kinds of factorizations where we extracted homogeneous factors
1369    for (j = 1;j<=size(resulttemp[i]);j++)
1370    {//searching for the inhomogeneous factor
1371      if (!homogwithorder(resulttemp[i][j],intvec(-1,1)))
1372      {//We have found our candidate
1373        isAlreadyInInhomogList = 0;
1374        for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++)
1375        {//Checking if our candidate is already in our tofactorize-list
1376          if (inhomogeneousFactorsToFactorize[k]==resulttemp[i][j])
1377          {//The candidate was already in the list
1378            isAlreadyInInhomogList = 1;
1379            break;
1380          }//The candidate was already in the list
1381        }//Checking if our candidate is already in our tofactorize-list
1382        if (!isAlreadyInInhomogList)
1383        {
1384          inhomogeneousFactorsToFactorize=inhomogeneousFactorsToFactorize + list(resulttemp[i][j]);
1385        }
1386      }//We have found our candidate
1387    }//searching for the inhomogeneous factor
1388  }//Going through all different kinds of factorizations where we extracted homogeneous factors
1389  dbprint(p,dbprintWhitespace +" Done");
1390  dbprint(p,dbprintWhitespace + "The set is:");
1391  dbprint(p,inhomogeneousFactorsToFactorize);
1392  dbprint(p,dbprintWhitespace+ "Factorizing the different occuring inhomogeneous factors");
1393  for (i = 1; i<= size(inhomogeneousFactorsToFactorize); i++)
1394  {//Factorizing all kinds of inhomogeneous factors
1395    inhomogeneousFactorsToFactorize[i] = sfacwa2(inhomogeneousFactorsToFactorize[i]);
1396    for (j = 1; j<=size(inhomogeneousFactorsToFactorize[i]);j++)
1397    {//Deleting the leading coefficient since we don't need him
1398      if (deg(inhomogeneousFactorsToFactorize[i][j][1],intvec(1,1))==0)
1399      {
1400        inhomogeneousFactorsToFactorize[i][j] = delete(inhomogeneousFactorsToFactorize[i][j],1);
1401      }
1402    }//Deleting the leading coefficient since we don't need him
1403  }//Factorizing all kinds of inhomogeneous factors
1404  dbprint(p,dbprintWhitespace +" Done");
1405  dbprint(p,dbprintWhitespace +" Putting the factorizations in the lists");
1406  list result;
1407  int posInhomogPoly;
1408  int posInhomogFac;
1409  for (i = 1; i<=size(resulttemp); i++)
1410  {//going through all by now calculated factorizations
1411    for (j = 1;j<=size(resulttemp[i]); j++)
1412    {//Finding the inhomogeneous factor
1413      if (!homogwithorder(resulttemp[i][j],intvec(-1,1)))
1414      {//Found it
1415        posInhomogPoly = j;
1416        break;
1417      }//Found it
1418    }//Finding the inhomogeneous factor
1419    for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++)
1420    {//Finding the matching inhomogeneous factorization we already determined
1421      if(product(inhomogeneousFactorsToFactorize[k][1]) == resulttemp[i][j])
1422      {//found it
1423        posInhomogFac = k;
1424        break;
1425      }//Found it
1426    }//Finding the matching inhomogeneous factorization we already determined
1427    for (j = 1; j <= size(inhomogeneousFactorsToFactorize[posInhomogFac]); j++)
1428    {
1429      result = insert(result, resulttemp[i]);
1430      result[1] = delete(result[1],posInhomogPoly);
1431      for (k =size(inhomogeneousFactorsToFactorize[posInhomogFac][j]);k>=1; k--)
1432      {//Inserting factorizations
1433        result[1] = insert(result[1],inhomogeneousFactorsToFactorize[posInhomogFac][j][k],
1434                           posInhomogPoly-1);
1435      }//Inserting factorizations
1436    }
1437  }//going through all by now calculated factorizations
1438  dbprint(p,dbprintWhitespace +" Done");
1439  result = delete_dublicates_noteval(result);
1440  for (i = 1; i<=size(result);i++)
1441  {//Putting the content everywhere
1442    result[i] = insert(result[i],commonCoefficient);
1443  }//Putting the content everywhere
1444  return(result);
1445}//proc sfacwa
1446
1447static proc sfacwa2(poly h)
1448"
1449Subprocedure of sfacwa
1450Assumptions:
1451- h is not in K[x] or in K[d], or even in K. These cases are caught by the input
1452- The coefficients are integer values and the gcd of the coefficients is 1
1453"
1454{//proc sfacwa2
1455  int p=printlevel-voice+2; // for dbprint
1456  int i;
1457  string dbprintWhitespace = "";
1458  for (i = 1; i<=voice;i++)
1459  {dbprintWhitespace = dbprintWhitespace + " ";}
1460  intvec ivm11 = intvec(-1,1);
1461  intvec iv11 = intvec(1,1);
1462  intvec iv10 = intvec(1,0);
1463  intvec iv01 = intvec(0,1);
1464  intvec iv1m1 = intvec(1,-1);
1465  poly p_max; poly p_min; poly q_max; poly q_min;
1466  map invo = basering,-var(1),var(2);
1467  list calculatedRightFactors;
1468  if(homogwithorder(h,ivm11))
1469  {//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here
1470    dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous.
1471Start homog. fac. and ret. its result");
1472    return(homogfacFirstWeyl_all(h));
1473  }//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here
1474  list result = list();
1475  int j; int k; int l;
1476  dbprint(p,dbprintWhitespace+" Computing the degree-limits of the factorization");
1477  //end finding the limits
1478  dbprint(p,dbprintWhitespace+" Computing the maximal and the minimal
1479homogeneous part of the given polynomial");
1480  list M = computeCombinationsMinMaxHomog(h);
1481  dbprint(p,dbprintWhitespace+" Done.");
1482  dbprint(p,dbprintWhitespace+" Filtering invalid combinations in M.");
1483  for (i = 1 ; i<= size(M); i++)
1484  {//filter valid combinations
1485    if (product(M[i]) == h)
1486    {//We have one factorization
1487      result = result + divides(M[i][1],h,invo,1);
1488      dbprint(p,dbprintWhitespace+"Result list updated:");
1489      dbprint(p,dbprintWhitespace+string(result));
1490      M = delete(M,i);
1491      continue;
1492    }//We have one factorization
1493  }//filter valid combinations
1494  dbprint(p,dbprintWhitespace+"Done.");
1495  dbprint(p,dbprintWhitespace+"The size of M is "+string(size(M)));
1496  for (i = 1; i<=size(M); i++)
1497  {//Iterate over all first combinations (p_max + p_min)(q_max + q_min)
1498    dbprint(p,dbprintWhitespace+" Combination No. "+string(i)+" in M:" );
1499    p_max = jet(M[i][1],deg(M[i][1],ivm11),ivm11)-jet(M[i][1],deg(M[i][1],ivm11)-1,ivm11);
1500    p_min = jet(M[i][1],deg(M[i][1],iv1m1),iv1m1)-jet(M[i][1],deg(M[i][1],iv1m1)-1,iv1m1);
1501    q_max = jet(M[i][2],deg(M[i][2],ivm11),ivm11)-jet(M[i][2],deg(M[i][2],ivm11)-1,ivm11);
1502    q_min = jet(M[i][2],deg(M[i][2],iv1m1),iv1m1)-jet(M[i][2],deg(M[i][2],iv1m1)-1,iv1m1);
1503    dbprint(p,dbprintWhitespace+" pmax = "+string(p_max));
1504    dbprint(p,dbprintWhitespace+" pmin = "+string(p_min));
1505    dbprint(p,dbprintWhitespace+" qmax = "+string(q_max));
1506    dbprint(p,dbprintWhitespace+" qmin = "+string(q_min));
1507    //Check, whether p_max + p_min or q_max and q_min are already left or right divisors.
1508    if (divides(p_min + p_max,h,invo))
1509    {
1510      dbprint(p,dbprintWhitespace+" Got one result.");
1511      result = result + divides(p_min + p_max,h,invo,1);
1512    }
1513    else
1514    {
1515      if (divides(q_min + q_max,h,invo))
1516      {
1517        dbprint(p,dbprintWhitespace+" Got one result.");
1518        result = result + divides(q_min + q_max, h , invo, 1);
1519      }
1520    }
1521    //Now the check, if deg(p_max) = deg(p_min)+1 (and the same with q_max and q_min)
1522
1523    if (deg(p_max, ivm11) == deg(p_min, ivm11) +1 or  deg(q_max, ivm11) == deg(q_min, ivm11) +1 )
1524    {//Therefore, p_max + p_min must be a left factor or we can dismiss the combination
1525      dbprint(p,dbprintWhitespace+" There are no homogeneous parts we can put between
1526pmax and pmin resp. qmax and qmin.");
1527      //TODO: Prove, that then also a valid right factor is not possible
1528      M = delete(M,i);
1529      continue;
1530    }//Therefore, p_max + p_min must be a left factor or we can dismiss the combination
1531
1532    //Done with the Check
1533
1534    //If we come here, there are still homogeneous parts to be added to p_max + p_min
1535    //AND to q_max and q_min in
1536    //order to obtain a real factor
1537    //We use the procedure determineRestOfHomogParts to find our q.
1538    dbprint(p,dbprintWhitespace+" Solving for the other homogeneous parts in q");
1539    calculatedRightFactors = determineRestOfHomogParts(p_max,p_min,q_max,q_min,h);
1540    dbprint(p,dbprintWhitespace+" Done with it. Found "+string(size(calculatedRightFactors))
1541            +" solutions.");
1542    for (j = 1; j<=size(calculatedRightFactors);j++)
1543    {//Check out whether we really have right factors of h in calculatedRightFactors
1544      if (divides(calculatedRightFactors[j],h,invo))
1545      {
1546        result = result + divides(calculatedRightFactors[j],h,invo,1);
1547      }
1548      else
1549      {
1550        dbprint(p,"Solution for max and min homog found, but not a divisor of h");
1551        //TODO: Proof, why this can happen.
1552      }
1553    }//Check out whether we really have right factors of h in calculatedRightFactors
1554  }//Iterate over all first combinations (p_max + p_min)(q_max + q_min)
1555
1556
1557  result = delete_dublicates_noteval(result);
1558  //print(M);
1559  if (size(result) == 0)
1560  {//no factorization found
1561    result = list(list(h));
1562  }//no factorization found
1563  return(result);
1564}//proc sfacwa2
1565
1566static proc determineRestOfHomogParts(poly pmax, poly pmin, poly qmax, poly qmin, poly h)
1567"INPUT: Polynomials p_max, p_min, q_max, q_min and h. The maximum homogeneous part h_max of h is
1568 given by p_max*pmin, the minimum homogeneous part h_min of h is given by p_min*q_min.
1569OUTPUT: A list of right factors q of h that have q_max and q_min as their maximum respectively
1570 minimum homogeneous part. Empty list, if those elements are not existent
1571ASSUMPTIONS:
1572 - deg(p_max,intvec(-1,1))>deg(p_min,intvec(-1,1)) +1
1573 - deg(q_max,intvec(-1,1))>deg(q_min,intvec(-1,1)) +1
1574 - p_max*q_max = h_max
1575 - p_min*q_min = h_min
1576 - The basering is the first Weyl algebra
1577"
1578{//proc determineRestOfHomogParts
1579  int p=printlevel-voice+2; // for dbprint
1580  string dbprintWhitespace = "";
1581  int i;
1582  for (i = 1; i<=voice;i++)
1583  {dbprintWhitespace = dbprintWhitespace + " ";}
1584  int kappa = Min(intvec(deg(h,intvec(1,0)), deg(h,intvec(0,1))));
1585  def R = basering;
1586  int n1 = deg(pmax,intvec(-1,1));
1587  int nk = -deg(pmin,intvec(1,-1));
1588  int m1 = deg(qmax,intvec(-1,1));
1589  int ml = -deg(qmin,intvec(1,-1));
1590  int j; int k;
1591  ideal mons;
1592
1593  dbprint(p,dbprintWhitespace+" Extracting zero homog. parts of pmax, qmax, pmin, qmin and h.");
1594  //Extracting the zero homogeneous part of the given polynomials
1595  ideal pandqZero = pmax,pmin,qmax,qmin;
1596  if (n1 > 0){pandqZero[1] = lift(var(2)^n1,pmax)[1,1];}
1597  else{if (n1 < 0){pandqZero[1] = lift(var(1)^(-n1),pmax)[1,1];}
1598    else{pandqZero[1] = pmax;}}
1599  if (nk > 0){pandqZero[2] = lift(var(2)^nk,pmin)[1,1];}
1600  else{if (nk < 0){pandqZero[2] = lift(var(1)^(-nk),pmin)[1,1];}
1601    else{pandqZero[2] = pmin;}}
1602  if (m1 > 0){pandqZero[3] = lift(var(2)^m1,qmax)[1,1];}
1603  else{if (m1 < 0){pandqZero[3] = lift(var(1)^(-m1),qmax)[1,1];}
1604    else{pandqZero[3] = qmax;}}
1605  if (ml > 0){pandqZero[4] = lift(var(2)^ml,qmin)[1,1];}
1606  else{if (ml < 0){pandqZero[4] = lift(var(1)^(-ml),qmin)[1,1];}
1607    else{pandqZero[4] = qmin;}}
1608  list hZeroinR = homogDistribution(h);
1609  for (i = 1; i<=size(hZeroinR);i++)
1610  {//Extracting the zero homogeneous parts of the homogeneous summands of h
1611    if (hZeroinR[i][1] > 0){hZeroinR[i][2] = lift(var(2)^hZeroinR[i][1],hZeroinR[i][2])[1,1];}
1612    if (hZeroinR[i][1] < 0){hZeroinR[i][2] = lift(var(1)^(-hZeroinR[i][1]),hZeroinR[i][2])[1,1];}
1613  }//Extracting the zero homogeneous parts of the homogeneous summands of h
1614  dbprint(p,dbprintWhitespace+" Done!");
1615  //Moving everything into the ring K[theta]
1616  dbprint(p,dbprintWhitespace+" Moving everything into the ring K[theta]");
1617  ring KTheta = 0,(x,d,theta),dp;
1618  map thetamap = R, x, d;
1619  poly entry;
1620  ideal mons;
1621  ideal pandqZero;
1622  list hZeroinKTheta;
1623  setring(R);
1624
1625  //Starting with p and q
1626  for (k=1; k<=4; k++)
1627  {//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys
1628    mons = ideal();
1629    for(i = 1; i<=size(pandqZero[k]);i++)
1630    {//Putting the monomials in a list
1631      mons[size(mons)+1] = pandqZero[k][i];
1632    }//Putting the monomials in a list
1633    setring(KTheta);
1634    mons = thetamap(mons);
1635    for (i = 1; i<=size(mons);i++)
1636    {//transforming the monomials as monomials in theta
1637      entry = leadcoef(mons[i]);
1638      for (j = 0; j<leadexp(mons[i])[2];j++)
1639      {
1640        entry = entry * (theta-j);
1641      }
1642      mons[i] = entry;
1643    }//transforming the monomials as monomials in theta
1644    pandqZero[size(pandqZero)+1] = sum(mons);
1645    setring(R);
1646  }//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys
1647
1648  //Now hZero
1649  for (k = size(hZeroinR); k>= 1;k--)
1650  {//Transforming the different homogeneous parts of h into polys in K[theta]
1651    mons = ideal();
1652    for(i = 1; i<=size(hZeroinR[k][2]);i++)
1653    {//Putting the monomials in a list
1654      mons[size(mons)+1] = hZeroinR[k][2][i];
1655    }//Putting the monomials in a list
1656    setring(KTheta);
1657    mons = thetamap(mons);
1658    for (i = 1; i<=size(mons);i++)
1659    {//transforming the monomials as monomials in theta
1660      entry = leadcoef(mons[i]);
1661      for (j = 0; j<leadexp(mons[i])[2];j++)
1662      {
1663        entry = entry * (theta-j);
1664      }
1665      mons[i] = entry;
1666    }//transforming the monomials as monomials in theta
1667    hZeroinKTheta = hZeroinKTheta + list(sum(mons));
1668    setring(R);
1669  }//Transforming the different homogeneous parts of h into polys in K[theta]
1670  dbprint(p,dbprintWhitespace+" Done!");
1671  //Making the solutionRing
1672  ring solutionRing = 0,(theta,q(0..(kappa+1)*(m1-ml-1)-1)),lp;
1673  dbprint(p,dbprintWhitespace+" Our solution ring is given by "+ string(solutionRing));
1674  //mapping the different ps and qs and HZeros
1675  dbprint(p,dbprintWhitespace+" Setting up our solution system.");
1676  list ps;
1677  ideal pandqZero = imap(KTheta,pandqZero);
1678  ps[1] = list(n1,pandqZero[1]);
1679  ps[n1-nk+1] = list(nk,pandqZero[2]);
1680  for (i = 2; i<=n1-nk; i++)
1681  {
1682    ps[i] = list(n1-i+1,0);
1683  }
1684  list qs;
1685  qs[1] = list(m1,pandqZero[3]);
1686  qs[m1-ml+1] = list(ml,pandqZero[4]);
1687  for (i = 2; i<=m1-ml; i++)
1688  {
1689    qs[i] = list(m1-i+1,0);
1690    for (j = 0; j<=kappa; j++)
1691    {
1692      qs[i][2] = qs[i][2] + q((i-2)*(kappa+1)+j)*theta^j;
1693    }
1694  }
1695  list hZero = imap(KTheta,hZeroinKTheta);
1696  for (i = 1; i<=size(hZero); i++)
1697  {
1698    hZero[i] = list(n1+m1-i+1,hZero[i]);
1699  }
1700
1701  //writing and solving the system
1702  list lhs;
1703  lhs[1] = list(ps[1][2],1);
1704  list rhs;
1705  rhs[n1-nk+1] = list(ps[size(ps)][2],1);
1706  for (i = 2; i<= n1-nk+1;i++)
1707  {
1708    lhs[i] = list(hZero[i][2],subst(qs[1][2],theta,theta+n1-i+1)*gammaForTheta(ps[i][1],qs[1][1]));
1709    rhs[n1-nk-i+2] = list(hZero[size(hZero)-i+1][2],subst(qs[m1-ml+1][2],theta,
1710                          theta+nk+i-1)*gammaForTheta(ps[n1-nk-i+2][1],qs[m1-ml+1][1]));
1711    for (j = 1; j<i; j++)
1712    {
1713      for (k = 1; k<=size(qs);k++)
1714      {
1715        if(ps[j][1]+qs[k][1] == hZero[i][1])
1716        {
1717          lhs[i][1] = lhs[i][1]*lhs[j][2];
1718          lhs[i][2] = lhs[i][2]*lhs[j][2];
1719          lhs[i][1] = lhs[i][1]-lhs[j][1]*subst(qs[k][2],theta, theta + n1- j +1)
1720            *gammaForTheta(ps[j][1],qs[k][1]);
1721        }
1722        if(ps[n1-nk+2-j][1] + qs[m1-ml+2-k][1] ==hZero[size(hZero)-i+1][1])
1723        {
1724          rhs[n1-nk-i+2][1] = rhs[n1-nk-i+2][1]*rhs[n1-nk+2-j][2];
1725          rhs[n1-nk-i+2][2] = rhs[n1-nk-i+2][2]*rhs[n1-nk+2-j][2];
1726          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,
1727                             theta + nk -j+1)*gammaForTheta(ps[n1-nk+2-j][1],qs[m1-ml+2-k][1]);
1728        }
1729      }
1730    }
1731  }
1732  list eqs;
1733  poly tempgcd;
1734  poly templhscoeff;
1735  poly temprhscoeff;
1736  for (i = 2; i<=n1-nk;i++)
1737  {
1738    if (gcd(rhs[i][2],lhs[i][2]) == 1)
1739    {
1740      eqs = eqs + list(lhs[i][1]*rhs[i][2] -rhs[i][1]*lhs[i][2]);
1741    }
1742    else
1743    {
1744      tempgcd = gcd(rhs[i][2],lhs[i][2]);
1745      templhscoeff = quotient(rhs[i][2],tempgcd)[1];
1746      temprhscoeff = quotient(lhs[i][2],tempgcd)[1];
1747      eqs = eqs + list(lhs[i][1]*templhscoeff -rhs[i][1]*temprhscoeff);
1748    }
1749  }
1750  matrix tempCoefMatrix;
1751  ideal solutionSystemforqs;
1752  for (i= 1; i<=size(eqs); i++)
1753  {
1754    tempCoefMatrix = coef(eqs[i],theta);
1755    solutionSystemforqs = solutionSystemforqs + ideal(submat(tempCoefMatrix,intvec(2),
1756                                                             intvec(1..ncols(tempCoefMatrix))));
1757  }
1758  dbprint(p,dbprintWhitespace+" Solution system for the coefficients of q is given by:");
1759  dbprint(p,solutionSystemforqs);
1760  option(redSB);
1761  dbprint(p,dbprintWhitespace+" Calculating reduced Groebner Basis of that system.");
1762  solutionSystemforqs = slimgb(solutionSystemforqs);
1763  dbprint(p,dbprintWhitespace+" Done!, the solution for the system is:");
1764  dbprint(p,dbprintWhitespace+string(solutionSystemforqs));
1765  if(vdim(slimgb(solutionSystemforqs+theta))==0)
1766  {//No solution in this case. Return the empty list
1767    dbprint(p,dbprintWhitespace+"The Groebner Basis of the solution system was <1>.");
1768    setring(R);
1769    return(list());
1770  }//No solution in this case. Return the empty list
1771  if(vdim(slimgb(solutionSystemforqs+theta))==-1)
1772  {//My conjecture is that this would never happen
1773    //ERROR("This is an counterexample to your conjecture. We have infinitely many solutions");
1774    //TODO: See, what we would do here
1775    dbprint(p,dbprintWhitespace+"There are infinitely many solution to this system.
1776We will return the empty list.");
1777    setring(R);
1778    return(list());
1779  }//My conjecture is that this would never happen
1780  else
1781  {//We have finitely many solutions
1782    if(vdim(slimgb(solutionSystemforqs+theta))==1)
1783    {//exactly one solution
1784      for (i = 2; i<= size(qs)-1;i++)
1785      {
1786        qs[i][2] = NF(qs[i][2],solutionSystemforqs);
1787      }
1788      setring(R);
1789      map backFromSolutionRing = solutionRing,var(1)*var(2);
1790      list qs = backFromSolutionRing(qs);
1791      list result = list(0);
1792      for (i = 1; i<=size(qs); i++)
1793      {
1794        if (qs[i][1]>0){qs[i][2] = qs[i][2]*var(2)^qs[i][1];}
1795        if (qs[i][1]<0){qs[i][2] = qs[i][2]*var(1)^(-qs[i][1]);}
1796        result[1] = result[1] + qs[i][2];
1797      }
1798      dbprint(p,dbprintWhitespace+"Found one unique solution. Returning the result.");
1799      return(result);
1800    }//exactly one solution
1801    else
1802    {//We have more than one solution, but finitely many
1803      dbprint(p,dbprintWhitespace+"Finitely many, but more than one solution.
1804Right now the conjecture is that this cannot happen.");
1805      dbprint(p,dbprintWhitespace+"Yay to counterexample :-) We can have more than one solution.");
1806      //TODO: Some theoretical work why this can never happan
1807      return(list());
1808    }//We have more than one solution, but finitely many
1809  }//We have finitely many solutions
1810}//proc determineRestOfHomogParts
1811
1812static proc gammaForTheta(int j1,int j2)
1813"
1814INPUT: Two integers j1 and j2
1815OUTPUT: A polynomial in the first variable of the given ring. It calculates the following function:
1816   / 1,                                                if j1,j2>0 or j1,j2 <= 0
1817   | prod_{kappa = 0}^{|j1|-1}(var(1)-kappa),          if j1<0, j2>0, |j1|<=|j2|
1818gamma_{j1,j2}:= <  prod_{kappa = 0}^{j2-1}(var(1)-kappa-|j1|+|j2|),  if j1<0, j2>0, |j1|>|j2|
1819   | prod_{kappa = 1}^{j1}(var(1)+kappa),              if j1>0, j2<0, |j1|<=|j2|
1820   \ prod_{kappa = 1}^{|j2|}(\var(1)+kappa+|j1|-|j2|), if j1>0, j2<0, |j1|>|j2|
1821ASSUMPTION:
1822- Ring has at least one variable
1823"
1824{//gammaForTheta
1825  if (j1<=0 && j2 <=0) {return(1);}
1826  if (j1>=0 && j2 >=0){return(1);}
1827  poly result;
1828  int i;
1829  if (j1<0 && j2>0)
1830  {//case 2 or 3 from description above
1831    if (absValue(j1)<=absValue(j2))
1832    {//Case 2 holds here
1833      result = 1;
1834      for (i = 0;i<absValue(j1);i++)
1835      {
1836        result = result*(var(1)-i);
1837      }
1838      return(result);
1839    }//Case 2 holds here
1840    else
1841    {//Case 3 holds here
1842      result = 1;
1843      for (i = 0; i<j2; i++)
1844      {
1845        result = result*(var(1)-i-absValue(j1)+absValue(j2));
1846      }
1847      return(result);
1848    }//Case 3 holds here
1849  }//case 2 or 3 from description above
1850  else
1851  {//Case 4 or 5 from description above hold
1852    if (absValue(j1)<=absValue(j2))
1853    {//Case 4 holds
1854      result = 1;
1855      for (i = 1; i<=j1; i++)
1856      {
1857        result = result*(var(1)+i);
1858      }
1859      return(result);
1860    }//Case 4 holds
1861    else
1862    {//Case 5 holds
1863      result = 1;
1864      for (i = 1; i<=absValue(j2); i++)
1865      {
1866        result = result*(var(1)+i+absValue(j1)-absValue(j2));
1867      }
1868      return(result);
1869    }//Case 5 holds
1870  }//Case 4 or 5 from description above hold
1871}//gammaForTheta
1872
1873static proc extractHomogeneousDivisors(poly h)
1874"INPUT: A polynomial h in the first Weyl algebra
1875OUTPUT: If h is homogeneous, then all factorizations of h are returned.
1876 If h is inhomogeneous, then a list l is returned whose entries
1877 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h and there
1878 exists an i in {1,...,n}, such that k_i is inhomogeneous and there
1879 is no homogeneous polynomial that divides this k_i neither from the
1880 left nor from the right. All the other entries in k are homogeneous
1881 polynomials.
1882"
1883{//extractHomogeneousDivisors
1884  int p=printlevel-voice+2; // for dbprint
1885  string dbprintWhitespace = "";
1886  int i; int j; int k; int l;
1887  list result;
1888  for (i = 1; i<=voice;i++)
1889  {dbprintWhitespace = dbprintWhitespace + " ";}
1890  if (homogwithorder(h,intvec(-1,1)))
1891  {//given polynomial was homogeneous already
1892    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
1893    result = homogfacFirstWeyl_all(h);
1894    for (i = 1; i<=size(result);i++)
1895    {//removing the first entry (coefficient) from the list result
1896      result[i] = delete(result[i],1);
1897    }//removing the first entry (coefficient) from the list result
1898    return(result);
1899  }//given polynomial was homogeneous already
1900  dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous left divisors extracted");
1901  list leftDivisionPossibilities = extractHomogeneousDivisorsLeft(h);
1902  dbprint(p,dbprintWhitespace+"Done. The result is:");
1903  dbprint(p,leftDivisionPossibilities);
1904  dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous Right divisors extracted");
1905  list rightDivisionPossibilities = extractHomogeneousDivisorsRight(h);
1906  dbprint(p,dbprintWhitespace+"Done. The result is:");
1907  dbprint(p,rightDivisionPossibilities);
1908  list tempList;
1909  dbprint(p,dbprintWhitespace+"Calculating remaining right and left homogeneous divisors");
1910  for (i = 1; i<=size(leftDivisionPossibilities); i++)
1911  {//iterating through the list with extracted left divisors
1912    tempList = extractHomogeneousDivisorsRight
1913      (leftDivisionPossibilities[i][size(leftDivisionPossibilities[i])]);
1914    leftDivisionPossibilities[i] = delete(leftDivisionPossibilities[i],
1915                                          size(leftDivisionPossibilities[i]));
1916    for (j=1;j<=size(tempList);j++)
1917    {//Updating the list for Result
1918      tempList[j] = leftDivisionPossibilities[i] + tempList[j];
1919    }//Updating the list for Result
1920    result = result + tempList;
1921  }//iterating through the list with extracted left divisors
1922  for (i = 1; i<=size(rightDivisionPossibilities); i++)
1923  {//iterating through the list with extracted left divisors
1924    tempList = extractHomogeneousDivisorsLeft(rightDivisionPossibilities[i][1]);
1925    rightDivisionPossibilities[i] = delete(rightDivisionPossibilities[i],1);
1926    for (j=1;j<=size(tempList);j++)
1927    {//Updating the list for Result
1928      tempList[j] =  tempList[j]+rightDivisionPossibilities[i];
1929    }//Updating the list for Result
1930    result = result + tempList;
1931  }//iterating through the list with extracted left divisors
1932  dbprint(p,dbprintWhitespace+"Done");
1933  int posInhomog;
1934  poly hath = 1;
1935  list tempResult;
1936  dbprint(p,dbprintWhitespace+"Checking if we can swap left resp. right
1937divisors and updating result.");
1938  for (i = 1; i<= size(result); i++)
1939  {//Checking if we can swap left resp. right divisors
1940    for (j = 1; j<=size(result[i]);j++)
1941    {//finding the position of the inhomogeneous element in the list
1942      if(!homogwithorder(result[i][j],intvec(-1,1)))
1943      {
1944        posInhomog = j;
1945        break;
1946      }
1947    }//finding the position of the inhomogeneous element in the list
1948    hath = result[i][posInhomog];
1949    for(j=posInhomog-1;j>=1;j--)
1950    {
1951      hath = result[i][j]*hath;
1952      tempList = extractHomogeneousDivisorsRight(hath);
1953      if(size(tempList[1])==1)
1954      {//We could not swap this element to the right
1955        break;
1956      }//We could not swap this element to the right
1957      dbprint(p,dbprintWhitespace+"A swapping (left) of an element was possible");
1958      for(k = 1; k<=size(tempList);k++)
1959      {
1960        tempResult = insert(tempResult,result[i]);
1961        for (l = j;l<=posInhomog;l++)
1962        {
1963          tempResult[1] = delete(tempResult[1],j);
1964        }
1965        for (l = size(tempList[k]);l>=1;l--)
1966        {
1967          tempResult[1] = insert(tempResult[1],tempList[k][l],j-1);
1968        }
1969      }
1970    }
1971    hath = result[i][posInhomog];
1972    for(j=posInhomog+1;j<=size(result[i]);j++)
1973    {
1974      hath = hath*result[i][j];
1975      tempList = extractHomogeneousDivisorsLeft(hath);
1976      if(size(tempList[1])==1)
1977      {//We could not swap this element to the right
1978        break;
1979      }//We could not swap this element to the right
1980      dbprint(p,dbprintWhitespace+"A swapping (right) of an element was possible");
1981      for(k = 1; k<=size(tempList);k++)
1982      {
1983        tempResult = insert(tempResult,result[i]);
1984        for (l=posInhomog; l<=j;l++)
1985        {
1986          tempResult[1] = delete(tempResult[1],posInhomog);
1987        }
1988        for (l = size(tempList[k]);l>=1;l--)
1989        {
1990          tempResult[1] = insert(tempResult[1],tempList[k][l],posInhomog-1);
1991        }
1992      }
1993    }
1994  }//Checking if we can swap left resp. right divisors
1995  result = result + tempResult;
1996  result = delete_dublicates_noteval(result);
1997  return(result);
1998}//extractHomogeneousDivisors
1999
2000static proc extractHomogeneousDivisorsLeft(poly h)
2001"INPUT: A polynomial h in the first Weyl algebra
2002OUTPUT: If h is homogeneous, then all factorizations of h are returned.
2003 If h is inhomogeneous, then a list l is returned whose entries
2004 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h.
2005 The entry k_n is inhomogeneous and has no other homogeneous
2006 left divisors any more.
2007 All the other entries in k are homogeneous
2008 polynomials.
2009"
2010{//extractHomogeneousDivisorsLeft
2011  int p=printlevel-voice+2; // for dbprint
2012  string dbprintWhitespace = "";
2013  int i;int j; int k;
2014  list result;
2015  poly hath;
2016  list recResult;
2017  map invo = basering,-var(1),var(2);
2018  for (i = 1; i<=voice;i++)
2019  {dbprintWhitespace = dbprintWhitespace + " ";}
2020  if (homogwithorder(h,intvec(-1,1)))
2021  {//given polynomial was homogeneous already
2022    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
2023    result = homogfacFirstWeyl_all(h);
2024    for (i = 1; i<=size(result);i++)
2025    {//removing the first entry (coefficient) from the list result
2026      result[i] = delete(result[i],1);
2027    }//removing the first entry (coefficient) from the list result
2028    return(result);
2029  }//given polynomial was homogeneous already
2030  list hlist = homogDistribution(h);
2031  dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands.");
2032  for (i = 1; i<= size(hlist); i++)
2033  {
2034    hlist[i] = homogfacFirstWeyl_all(hlist[i][2]);
2035    if (size(hlist[i][1])==1)
2036    {//One homogeneous part just has a trivial factorization
2037      if(hlist[i][1][1] == 0)
2038      {
2039        hlist = delete(hlist,i);
2040        continue;
2041      }
2042      else
2043      {
2044        return(list(list(h)));
2045      }
2046    }//One homogeneous part just has a trivial factorization
2047  }
2048  dbprint(p,dbprintWhitespace+ " Done.");
2049  dbprint(p,dbprintWhitespace+ " Trying to find Left divisors");
2050  list alreadyConsideredCandidates;
2051  poly candidate;
2052  int isCandidate;
2053  for (i = 1; i<=size(hlist[1]);i++)
2054  {//Finding candidates for homogeneous left divisors of h
2055    candidate = hlist[1][i][2];
2056    isCandidate = 0;
2057    for (j=1;j<=size(alreadyConsideredCandidates);j++)
2058    {
2059      if(alreadyConsideredCandidates[j] == candidate)
2060      {
2061        isCandidate =1;
2062        break;
2063      }
2064    }
2065    if(isCandidate)
2066    {
2067      i++;
2068      continue;
2069    }
2070    else
2071    {
2072      alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate);
2073    }
2074    dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous left divisor");
2075    for (j = 2; j<=size(hlist);j++)
2076    {//Iterating through the other homogeneous parts
2077      isCandidate = 0;
2078      for(k=1; k<=size(hlist[j]);k++)
2079      {
2080        if(hlist[j][k][2]==candidate)
2081        {
2082          isCandidate = 1;
2083          break;
2084        }
2085      }
2086      if(!isCandidate)
2087      {
2088        break;
2089      }
2090    }//Iterating through the other homogeneous parts
2091    if(isCandidate)
2092    {//candidate was really a left divisor
2093      dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous left divisor");
2094      hath = involution(lift(involution(candidate,invo),involution(h,invo))[1,1],invo);
2095      recResult = extractHomogeneousDivisorsLeft(hath);
2096      for (j = 1; j<=size(recResult); j++)
2097      {
2098        recResult[j] = insert(recResult[j],candidate);
2099      }
2100      result = result + recResult;
2101    }//Candidate was really a left divisor
2102  }//Finding candidates for homogeneous left divisors of h
2103  if (size(result)==0)
2104  {
2105    return(list(list(h)));
2106  }
2107  return(result);
2108}//extractHomogeneousDivisorsLeft
2109
2110static proc extractHomogeneousDivisorsRight(poly h)
2111"INPUT: A polynomial h in the first Weyl algebra
2112OUTPUT: If h is homogeneous, then all factorizations of h are returned.
2113 If h is inhomogeneous, then a list l is returned whose entries
2114 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h.
2115 The entry k_1 is inhomogeneous and has no other homogeneous
2116 right divisors any more.
2117 All the other entries in k are homogeneous
2118 polynomials.
2119"
2120{//extractHomogeneousDivisorsRight
2121  int p=printlevel-voice+2; // for dbprint
2122  string dbprintWhitespace = "";
2123  int i;int j; int k;
2124  list result;
2125  poly hath;
2126  list recResult;
2127  for (i = 1; i<=voice;i++)
2128  {dbprintWhitespace = dbprintWhitespace + " ";}
2129  if (homogwithorder(h,intvec(-1,1)))
2130  {//given polynomial was homogeneous already
2131    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
2132    result = homogfacFirstWeyl_all(h);
2133    for (i = 1; i<=size(result);i++)
2134    {//removing the first entry (coefficient) from the list result
2135      result[i] = delete(result[i],1);
2136    }//removing the first entry (coefficient) from the list result
2137    return(result);
2138  }//given polynomial was homogeneous already
2139  list hlist = homogDistribution(h);
2140  dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands.");
2141  for (i = 1; i<= size(hlist); i++)
2142  {
2143    hlist[i] = homogfacFirstWeyl_all(hlist[i][2]);
2144    if (size(hlist[i][1])==1)
2145    {//One homogeneous part just has a trivial factorization
2146      if(hlist[i][1][1] == 0)
2147      {
2148        hlist = delete(hlist,i);
2149        continue;
2150      }
2151      else
2152      {
2153        return(list(list(h)));
2154      }
2155    }//One homogeneous part just has a trivial factorization
2156  }
2157  dbprint(p,dbprintWhitespace+ " Done.");
2158  dbprint(p,dbprintWhitespace+ " Trying to find right divisors");
2159  list alreadyConsideredCandidates;
2160  poly candidate;
2161  int isCandidate;
2162  for (i = 1; i<=size(hlist[1]);i++)
2163  {//Finding candidates for homogeneous left divisors of h
2164    candidate = hlist[1][i][size(hlist[1][i])];
2165    isCandidate = 0;
2166    for (j=1;j<=size(alreadyConsideredCandidates);j++)
2167    {
2168      if(alreadyConsideredCandidates[j] == candidate)
2169      {
2170        isCandidate =1;
2171        break;
2172      }
2173    }
2174    if(isCandidate)
2175    {
2176      i++;
2177      continue;
2178    }
2179    else
2180    {
2181      alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate);
2182    }
2183    dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous r-divisor");
2184    for (j = 2; j<=size(hlist);j++)
2185    {//Iterating through the other homogeneous parts
2186      isCandidate = 0;
2187      for(k=1; k<=size(hlist[j]);k++)
2188      {
2189        if(hlist[j][k][size(hlist[j][k])]==candidate)
2190        {
2191          isCandidate = 1;
2192          break;
2193        }
2194      }
2195      if(!isCandidate)
2196      {
2197        break;
2198      }
2199    }//Iterating through the other homogeneous parts
2200    if(isCandidate)
2201    {//candidate was really a left divisor
2202      dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous right divisor");
2203      hath = lift(candidate,h)[1,1];
2204      recResult = extractHomogeneousDivisorsRight(hath);
2205      for (j = 1; j<=size(recResult); j++)
2206      {
2207        recResult[j] = insert(recResult[j],candidate,size(recResult[j]));
2208      }
2209      result = result + recResult;
2210    }//Candidate was really a left divisor
2211  }//Finding candidates for homogeneous left divisors of h
2212  if (size(result)==0)
2213  {
2214    result = list(list(h));
2215  }
2216  return(result);
2217}//extractHomogeneousDivisorsRight
2218
2219static proc fromZeroHomogToThetaPoly(poly h)
2220"
2221//DEPRECATED
2222INPUT: A polynomial h in the first Weyl algebra, homogeneous of degree 0
2223 OUTPUT: The ring Ktheta with a polynomial result representing h as polynomial in
2224  theta
2225 ASSUMPTIONS:
2226- h is homogeneous of degree 0 with respect to the [-1,1] weight vector
2227- The basering is the first Weyl algebra
2228"
2229{//proc fromZeroHomogToThetaPoly
2230  int i; int j;
2231  list mons;
2232  if(!homogwithorder(h,intvec(-1,1)))
2233  {//Input not a homogeneous polynomial causes an error
2234    ERROR("The input was not a homogeneous polynomial");
2235  }//Input not a homogeneous polynomial causes an error
2236  if(deg(h,intvec(-1,1))!=0)
2237  {//Input does not have degree 0
2238    ERROR("The input did not have degree 0");
2239  }//Input does not have degree 0
2240  for(i = 1; i<=size(h);i++)
2241  {//Putting the monomials in a list
2242    mons = mons+list(h[i]);
2243  }//Putting the monomials in a list
2244  ring KTheta = 0,(x,y,theta),dp;
2245  setring KTheta;
2246  map thetamap = r,x,y;
2247  list mons = thetamap(mons);
2248  poly entry;
2249  for (i = 1; i<=size(mons);i++)
2250  {//transforming the monomials as monomials in theta
2251    entry = leadcoef(mons[i]);
2252    for (j = 0; j<leadexp(mons[i])[2];j++)
2253    {
2254      entry = entry * (theta-j);
2255    }
2256    mons[i] = entry;
2257  }//transforming the monomials as monomials in theta
2258  poly result = sum(mons);
2259  keepring KTheta;
2260}//proc fromZeroHomogToThetaPoly
2261
2262static proc fromHomogDistributionListToPoly(list l)
2263"returns the corresponding polynomial to a given output of homogDistribution"
2264{//proc fromHomogDistributionListToPoly
2265  poly result = 0;
2266  int i;
2267  for (i = 1; i<=size(l);i++)
2268  {
2269    result = result + l[i][2];
2270  }
2271  return(result);
2272}//proc fromHomogDistributionListToPoly
2273
2274static proc divides(poly p1, poly p2,map invo, list #)
2275  "Tests, whether p1 divides p2 either from left or from right. The involution invo is needed
2276for checking both sides. The optional argument is needed in order to also return the other factor.
2277RETURN: If no optional argument is given, it will just return 1 or 0.
2278 Otherwise a list with at least one element
2279 Case 1: p1 does not divide p2 from any side. Then the output will be the empty list.
2280 Case 2: p2 does divide p2 from one side at least.
2281  Then it returns a list with tuples p,q, such that p or q equals p1 and
2282  pq = p2.
2283ASSUMPTIONS: - The map invo is an involution on the basering."
2284{//proc divides
2285  list result = list();
2286  poly tempfactor;
2287  if (involution(reduce(involution(p2,invo),std(involution(ideal(p1),invo))),invo)==0)
2288  {//p1 is a left divisor
2289    if(size(#)==0){return(1);}
2290    tempfactor = involution(lift(involution(p1,invo),involution(p2,invo))[1,1],invo);
2291    if (leadcoef(p1)<0 && leadcoef(tempfactor)<0)
2292    {//both have a negative leading coefficient
2293      result = result +list(list(1,-p1, -tempfactor));
2294    }//both have a negative leading coefficient
2295    else
2296    {
2297      if (leadcoef(p1)*leadcoef(tempfactor)<0)
2298      {//One of them has a negative leading coefficient
2299        if (leadcoef(p1)<0)
2300        {
2301          result = result + list(list(-1,-p1,tempfactor));
2302        }
2303        else
2304        {
2305          result = result + list(list(-1,p1,-tempfactor));
2306        }
2307      }//One of them has a negative leading coefficient
2308      else
2309      {//no negative coefficient at all
2310        result = result + list(list(1,p1,tempfactor));
2311      }//no negative coefficient at all
2312    }
2313  }//p1 is a left divisor
2314
2315  if (reduce(p2,std(ideal(p1))) == 0)
2316  {//p1 is  a right divisor
2317    if(size(#)==0){return(1);}
2318    tempfactor = lift(p1, p2)[1,1];
2319    if (leadcoef(p1)<0 && leadcoef(tempfactor)<0)
2320    {//both have a negative leading coefficient
2321      result = result +list(list(1, -tempfactor,-p1));
2322    }//both have a negative leading coefficient
2323    else
2324    {
2325      if (leadcoef(p1)*leadcoef(tempfactor)<0)
2326      {//One of them has a negative leading coefficient
2327        if (leadcoef(p1)<0)
2328        {
2329          result = result + list(list(-1,tempfactor,-p1));
2330        }
2331        else
2332        {
2333          result = result + list(list(-1,-tempfactor,p1));
2334        }
2335      }//One of them has a negative leading coefficient
2336      else
2337      {//no negative coefficient at all
2338        result = result + list(list(1,tempfactor,p1));
2339      }//no negative coefficient at all
2340    }
2341  }//p1 is already a right divisor
2342  if (size(#)==0){return(0);}
2343  return(result);
2344}//proc divides
2345
2346static proc computeCombinationsMinMaxHomog(poly h)
2347"Input: A polynomial h in the first Weyl Algebra
2348Output: Combinations of the form (p_max + p_min)(q_max + q_min), such that p_max, p_min,
2349 q_max and q_min are homogeneous and p_max*q_max equals the maximal homogeneous
2350 part in h, and p_max * q_max equals the minimal homogeneous part in h.
2351 h is not homogeneous.
2352"{//proc computeCombinationsMinMaxHomog
2353  intvec ivm11 = intvec(-1,1);
2354  intvec iv11 = intvec(1,1);
2355  intvec iv10 = intvec(1,0);
2356  intvec iv01 = intvec(0,1);
2357  intvec iv1m1 = intvec(1,-1);
2358  poly maxh = jet(h,deg(h,ivm11),ivm11)-jet(h,deg(h,ivm11)-1,ivm11);
2359  poly minh = jet(h,deg(h,iv1m1),iv1m1)-jet(h,deg(h,iv1m1)-1,iv1m1);
2360  list f1 = homogfacFirstWeyl_all(maxh);
2361  list f2 = homogfacFirstWeyl_all(minh);
2362  list result = list();
2363  int i;
2364  int j;
2365  //TODOs: Nur die kombinieren, die auch wirklich eine brauchbare Kombination sind.
2366  //Also, wir duerfen nicht aus den Geraden herausfliegen
2367  list pqmax = list();
2368  list pqmin = list();
2369  list tempList = list();
2370
2371  //First, we are going to deal with our most hated guys: The Coefficients.
2372  //
2373  list coeffTuplesMax = getAllCoeffTuplesComb(factorizeInt(number(f1[1][1])));
2374  //We can assume without loss of generality, that p_max has a
2375  //nonnegative leading coefficient
2376  for (i = 1; i<=size(coeffTuplesMax);i++)
2377  {//Deleting all tuples with negative entries for p_max
2378    if (coeffTuplesMax[i][1]<0)
2379    {
2380      coeffTuplesMax = delete(coeffTuplesMax,i);
2381      continue;
2382    }
2383  }//Deleting all tuples with negative entries for p_max
2384  list coeffTuplesMin = getAllCoeffTuplesComb(factorizeInt(number(f2[1][1])));
2385
2386  //Now, we will be actally dealing with the Combinations.
2387  //Let's start with the pqmax
2388  if (size(f1[1]) == 1)
2389  {//the maximal homogeneous factor is a constant
2390    pqmax = coeffTuplesMax;
2391  }//the maximal homogeneous factor is a constant
2392  else
2393  {//the maximal homogeneous factor is not a constant
2394    for (i = 1; i<=size(f1); i++)
2395    {//We can forget about the first coefficient now. Therefore we will delete him from the list.
2396      f1[i] = delete(f1[i],1);
2397      if(size(f1[i])==1)
2398      {//trivial thing
2399        for (j = 1; j<=size(coeffTuplesMax); j++)
2400        {
2401          pqmax = pqmax + list(list(coeffTuplesMax[j][1],coeffTuplesMax[j][2]*f1[i][1]));
2402          pqmax = pqmax + list(list(coeffTuplesMax[j][1]*f1[i][1],coeffTuplesMax[j][2]));
2403        }
2404        f1 = delete(f1,i);
2405        continue;
2406      }//trivial thing
2407    }//We can forget about the first coefficient now. Therefore we will delete him from the list.
2408    if (size(f1)>0)
2409    {
2410      tempList = getAllCombOfHomogFact(f1);
2411      for (i = 1; i<=size(tempList); i++)
2412      {//Every combination combined with the coefficient possibilities
2413        for (j = 1; j<=size(coeffTuplesMax); j++)
2414        {//iterating through the possible coefficient choices
2415          pqmax = pqmax + list(list(coeffTuplesMax[j][1]*tempList[i][1],
2416                                    coeffTuplesMax[j][2]*tempList[i][2]));
2417        }//iterating through the possible coefficient choices
2418      }//Every combination combined with the coefficient possibilities
2419      for (i = 1; i<=size(coeffTuplesMax); i++)
2420      {
2421        pqmax = pqmax + list(list(coeffTuplesMax[i][1],maxh/coeffTuplesMax[i][1]));
2422        pqmax = pqmax + list(list(maxh/coeffTuplesMax[i][2],coeffTuplesMax[i][2]));
2423      }
2424    }
2425  }//the maximal homogeneous factor is not a constant
2426  //Now we go to pqmin
2427  if (size(f2[1]) == 1)
2428  {//the minimal homogeneous factor is a constant
2429    pqmin = coeffTuplesMin;
2430  }//the minimal homogeneous factor is a constant
2431  else
2432  {//the minimal homogeneous factor is not a constant
2433    for (i = 1; i<=size(f2); i++)
2434    {//We can forget about the first coefficient now. Therefore we will delete him from the list.
2435      f2[i] = delete(f2[i],1);
2436      if(size(f2[i])==1)
2437      {//trivial thing
2438        for (j = 1; j<=size(coeffTuplesMin); j++)
2439        {
2440          pqmin = pqmin + list(list(coeffTuplesMin[j][1],coeffTuplesMin[j][2]*f2[i][1]));
2441          pqmin = pqmin + list(list(coeffTuplesMin[j][1]*f2[i][1],coeffTuplesMin[j][2]));
2442        }
2443        f2 = delete(f2,i);
2444        continue;
2445      }
2446    }//We can forget about the first coefficient now. Therefore we will delete him from the list.
2447    if(size(f2)>0)
2448    {
2449      tempList = getAllCombOfHomogFact(f2);
2450      for (i = 1; i<=size(tempList); i++)
2451      {//Every combination combined with the coefficient possibilities
2452        for (j = 1; j<=size(coeffTuplesMin); j++)
2453        {//iterating through the possible coefficient choices
2454          pqmin = pqmin + list(list(coeffTuplesMin[j][1]*tempList[i][1],
2455                                    coeffTuplesMin[j][2]*tempList[i][2]));
2456        }//iterating through the possible coefficient choices
2457      }//Every combination combined with the coefficient possibilities
2458      for (i = 1; i<=size(coeffTuplesMin); i++)
2459      {
2460        pqmin = pqmin + list(list(coeffTuplesMin[i][1],minh/coeffTuplesMin[i][1]));
2461        pqmin = pqmin + list(list(minh/coeffTuplesMin[i][2],coeffTuplesMin[i][2]));
2462      }
2463    }
2464  }//the minimal homogeneous factor is not a constant
2465
2466  //and now we combine them together to obtain all possibilities.
2467  for (i = 1; i<=size(pqmax); i++)
2468  {//iterate over the maximal homogeneous combination possibilities
2469    for (j = 1; j<=size(pqmin); j++)
2470    {//iterate over the minimal homogeneous combiniation possibilities
2471      if (deg(pqmax[i][1], ivm11)>=deg(pqmin[j][1],ivm11) and deg(pqmax[i][2],
2472                                                                  ivm11)>=deg(pqmin[j][2],ivm11))
2473      {
2474        if (pqmax[i][1]+pqmin[j][1]!=0 and pqmax[i][2]+pqmin[j][2]!=0)
2475        {
2476
2477          if (deg(h,ivm11)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),ivm11))
2478          {
2479            j++;
2480            continue;
2481          }
2482          if (deg(h,iv1m1)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),iv1m1))
2483          {
2484            j++;
2485            continue;
2486          }
2487          result = result +list(list(pqmax[i][1]+pqmin[j][1],pqmax[i][2]+pqmin[j][2]));
2488        }
2489      }
2490    }//iterate over the minimal homogeneous combiniation possibilities
2491  }//iterate over the maximal homogeneous combination possibilities
2492  //Now deleting double entries
2493  result = delete_dublicates_noteval(result);
2494  return(result);
2495}//proc computeCombinationsMinMaxHomog
2496
2497static proc getAllCombOfHomogFact(list l)
2498"Gets called in computeCombinationsMinMaxHomog. It gets a list of different homogeneous
2499factorizations of
2500one homogeneous polynomial and returns the possibilities to combine them into two factors.
2501Assumptions:
2502- The list does not contain the first coefficient.
2503- The list contains at least one list with two elements."
2504{//proc getAllCombOfHomogFact
2505  list result;
2506  list leftAndRightHandSides;
2507  int i; int j;
2508  list tempset;
2509  if (size(l)==1 and size(l[1])==2)
2510  {
2511    result = result + list(list(l[1][1],l[1][2]));
2512    return(result);
2513  }
2514  leftAndRightHandSides = getPossibilitiesForRightSides(l);
2515  for (i = 1; i<=size(leftAndRightHandSides); i++)
2516  {
2517    result =result+list(list(leftAndRightHandSides[i][1],product(leftAndRightHandSides[i][2][1])));
2518    //tidy up the right hand sides, because, if it is just one irreducible factor, we are done
2519    for (j = 1; j<=size(leftAndRightHandSides[i][2]);j++)
2520    {//Tidy up right hand sides
2521      if (size(leftAndRightHandSides[i][2][j])<2)
2522      {//Element can be dismissed
2523        leftAndRightHandSides[i][2] = delete(leftAndRightHandSides[i][2],j);
2524        continue;
2525      }//Element can be dismissed
2526    }//Tidy up right hand sides
2527    if (size(leftAndRightHandSides[i][2])>0)
2528    {
2529      tempset = getAllCombOfHomogFact(leftAndRightHandSides[i][2]);
2530      for (j = 1; j<=size(tempset);j++)
2531      {//multiplying the first factor with the left hand side
2532        result = result + list(list(leftAndRightHandSides[i][1]*tempset[j][1],tempset[j][2]));
2533      }//multiplying the first factor with the left hand side
2534    }
2535  }
2536  return(result);
2537}//proc getAllCombOfHomogFact
2538
2539static proc getPossibilitiesForRightSides(list l)
2540"Given a list of different factorizations l, this function returns a list of the form
2541(a,{(a_2,...,a_n)| (a,a_2,...,a_n) in A})"
2542{//getPossibilitiesForRightSide
2543  list templ = l;
2544  list result;
2545  poly firstElement;
2546  list rightSides;
2547  list tempRightSide;
2548  int i; int j;
2549  while (size(templ)>0)
2550  {
2551    firstElement = templ[1][1];
2552    rightSides = list();
2553    for (i = 1; i<= size(templ); i++)
2554    {
2555      if (templ[i][1] == firstElement)
2556      {//save the right sides
2557        tempRightSide = list();
2558        for (j = 2; j<=size(templ[i]);j++)
2559        {
2560          tempRightSide = tempRightSide + list(templ[i][j]);
2561        }
2562        if (size(tempRightSide)!=0)
2563        {
2564          rightSides = rightSides + list(tempRightSide);
2565        }
2566        templ = delete(templ,i);
2567        continue;
2568      }//save the right sides
2569    }
2570    result = result + list(list(firstElement,rightSides));
2571  }
2572  return(result);
2573}//getPossibilitiesForRightSide
2574
2575static proc getAllCoeffTuplesComb(list l)"
2576Given the output of factorizeInt ((a_1,...,a_n),(i_1,...,i_n)) , it returns all possible tuples
2577of the set {(a,b) | There exists an real N!=emptyset subset of {1,...,n}, such that
2578a = prod_{i \in N}a_i, b=prod_{i \not\in N} a_i}
2579Assumption: The list is sorted from smallest integer to highest.
2580- it is not the factorization of 0.
2581"
2582{//proc getAllCoeffTuplesComb
2583  list result;
2584  if (l[1][1] == 0)
2585  {
2586    ERROR("getAllCoeffTuplesComb: Zero Coefficients as leading and Tail Coeffs?
2587That is not possible. Something went wrong.");
2588  }
2589  if (size(l[1]) == 1)
2590  {//Trivial Factorization, just 1
2591    if (l[1][1] == 1)
2592    {
2593      return(list(list(1,1),list(-1,-1)));
2594    }
2595    else
2596    {
2597      return(list(list(-1,1),list(1,-1)));
2598    }
2599  }//Trivial Factorization, just 1
2600  if (size(l[1]) == 2 and l[2][2]==1)
2601  {//Just a prime number
2602    if (l[1][1] == 1)
2603    {
2604      result = list(list(l[1][2],1),list(1,l[1][2]));
2605      result = result + list(list(-l[1][2],-1),list(-1,-l[1][2]));
2606      return(result);
2607    }
2608    else
2609    {
2610      result = list(list(l[1][2],-1),list(1,-l[1][2]));
2611      result = result + list(list(-l[1][2],1),list(-1,l[1][2]));
2612      return(result);
2613    }
2614  }//Just a prime number
2615  //Now comes the interesting case: a product of primes
2616  list tempPrimeFactors;
2617  list tempPowersOfThem;
2618  int i;
2619  for (i = 2; i<=size(l[1]);i++)
2620  {//Removing the starting 1 or -1 to get the N's
2621    tempPrimeFactors[i-1] = l[1][i];
2622    tempPowersOfThem[i-1] = l[2][i];
2623  }//Removing the starting 1 or -1 to get the N's
2624  list Ns = getAllSubsetsN(list(tempPrimeFactors,tempPowersOfThem));
2625  list tempTuples;
2626  number productOfl = multiplyFactIntOutput(l);
2627  if (productOfl<0){productOfl = -productOfl;}
2628  tempTuples = tempTuples + list(list(1,productOfl),list(productOfl,1));
2629  for (i = 1; i<=size(Ns); i++)
2630  {
2631    if (productOfl/Ns[i]>Ns[i])
2632    {//TODO: BEWEISEN, dass das die einzigen Combos sind
2633      tempTuples = tempTuples + list(list(Ns[i],productOfl/Ns[i]),list(productOfl/Ns[i],Ns[i]));
2634    }//TODO: BEWEISEN, dass das die einzigen Combos sind
2635    if (productOfl/Ns[i]==Ns[i])
2636    {
2637      tempTuples = tempTuples + list(list(Ns[i],Ns[i]));
2638    }
2639  }
2640  //And now, it just remains to get the -1s and 1-s correctly to the tuples
2641  list tempEntry;
2642  if (l[1][1] == 1)
2643  {
2644    for (i = 1; i<=size(tempTuples);i++)
2645    {//Adding everything to result
2646      tempEntry = tempTuples[i];
2647      result = result + list(tempEntry);
2648      result = result + list(list(-tempEntry[1], -tempEntry[2]));
2649    }//Adding everyThing to Result
2650  }
2651  else
2652  {
2653    for (i = 1; i<=size(tempTuples);i++)
2654    {//Adding everything to result
2655      tempEntry = tempTuples[i];
2656      result = result + list(list(tempEntry[1],-tempEntry[2]));
2657      result = result + list(list(-tempEntry[1], tempEntry[2]));
2658    }//Adding everyThing to Result
2659  }
2660  return(result);
2661}//proc getAllCoeffTuplesComb
2662
2663static proc contains(list l, int elem)
2664"Assumption: l is sorted"
2665{//Binary Search in list
2666  if (size(l)<=1)
2667  {
2668    if(size(l) == 0){return(0);}
2669    if (l[1]!=elem){return(0);}
2670    else{return(1);}
2671  }
2672  int imax = size(l);
2673  int imin = 1;
2674  int imid;
2675  while(imax >= imin)
2676  {
2677    imid = (imin + imax)/2;
2678    if (l[imid] == elem){return(1);}
2679    if (l[imid] <elem) {imin = imid +1;}
2680    else{imax = imid -1;}
2681  }
2682  return(0)
2683    }//Binary Search in list
2684
2685static proc getAllSubsetsN(list l)
2686"
2687Assumptions:
2688- The list is containing two lists. They can be assumed to be outputs of the function
2689factorizeInt. They have at least one entry. If it is exactly one entry, the second intvec should
2690contain a value at least 2.
2691  "
2692{
2693  list primeFactors=l[1];
2694  list powersOfThem = l[2];
2695  int i;int j;
2696  //Casting the entries to be numbers
2697  for (i=1; i<=size(primeFactors); i++)
2698  {
2699    primeFactors[i] = number(primeFactors[i]);
2700    powersOfThem[i] = number(powersOfThem[i]);
2701  }
2702
2703  //Done
2704  list result;
2705  list tempPrimeFactors;
2706  list tempPowersOfThem;
2707  list tempset;
2708  if (sum(powersOfThem) <=2)
2709  {//Easy Case
2710    return(list(primeFactors[1]));
2711  }//Easy Case
2712  if (size(primeFactors)==1)
2713  {//Also Easy Case
2714    for (j = 1; j<powersOfThem[1]; j++)
2715    {
2716      result = result + list(primeFactors[1]^j);
2717    }
2718    return(result);
2719  }//Also Easy Case
2720  for (i = 1; i<= size(primeFactors); i++)
2721  {//Going through every entry
2722    result = result + list(primeFactors[i]);
2723    if (i == size(primeFactors))
2724    {
2725      for (j = 1;j<powersOfThem[i];j++)
2726      {
2727        result = result + list (primeFactors[i]^j);
2728      }
2729      break;
2730    }
2731    if (powersOfThem[i]==1)
2732    {
2733      for (j = i+1;j<=size(primeFactors);j++)
2734      {
2735        tempPrimeFactors[j-i] = primeFactors[j];
2736        tempPowersOfThem[j-i] = powersOfThem[j];
2737      }
2738    }
2739    else
2740    {
2741      for (j = i; j<=size(primeFactors);j++)
2742      {
2743        tempPrimeFactors[j-i+1] = primeFactors[j];
2744        tempPowersOfThem[j-i+1] = powersOfThem[j];
2745        tempPowersOfThem[1] = tempPowersOfThem[1]-1;
2746      }
2747    }
2748    tempset = getAllSubsetsN(list(tempPrimeFactors,tempPowersOfThem));
2749    for (j = 1; j<=size(tempset); j++)
2750    {
2751      result = result +list((tempset[j])*(primeFactors[i]));
2752    }
2753  }//Going through every entry
2754  result = sort(result)[1];
2755  result = delete_dublicates_noteval(result);
2756  return(result);
2757}
2758
2759static proc multiplyFactIntOutput(list l)
2760"Given the output of factorizeInt, this method computes the product of it."
2761{//proc multiplyFactIntOutput
2762  int i;
2763  number result = 1;
2764  for (i = 1; i<=size(l[1]); i++)
2765  {
2766    result = result*(l[1][i])^(l[2][i]);
2767  }
2768  return(result);
2769}//proc multiplyFactIntOutput
2770
2771static proc fromListToIntvec(list l)
2772"Converter from List to intvec"
2773{
2774  intvec result; int i;
2775  for (i = 1; i<=size(l); i++)
2776  {
2777    result[i] = l[i];
2778  }
2779  return(result);
2780}
2781
2782static proc fromIntvecToList(intvec l)"
2783Converter from intvec to list"
2784{//proc fromIntvecToList
2785  list result = list();
2786  int i;
2787  for (i = size(l); i>=1; i--)
2788  {
2789    result = insert(result, l[i]);
2790  }
2791  return(result);
2792}//proc fromIntvecToList
2793
2794
2795static proc factorizeInt(number n)
2796"Given an integer n, factorizeInt computes its factorization. The output is a list
2797containing two intvecs. The first contains the prime factors, the second its powers.
2798ASSUMPTIONS:
2799- n is given as integer number
2800"{
2801  if (n==0)
2802  {return(list(list(0),list(1)));}
2803  int i;
2804  list temp = primefactors(n);
2805  if (n<0)
2806  {list result = list(list(-1),list(1));}
2807  else
2808  {list result = list(list(1),list(1));}
2809  result[1] = result[1] + temp[1];
2810  result[2] = result[2] + temp[2];
2811  return(result);
2812}
2813
2814
2815static proc homogDistribution(poly h)
2816"Input: A polynomial in the first Weyl Algebra.
2817  Output: A two-dimensional list of the following form. Every sublist contains exactly two entries.
2818   One for the Z-degree of the corresponding homogeneous part (integer), and the homogeneous
2819   polynomial itself, and those sublists are oredered by ascending degree.
2820   For example a call of homogDistribution(x+d+1) would have the output
2821   [1]:
2822     [1]:
2823       -1
2824     [2]:
2825       x
2826   [2]:
2827     [1]:
2828       0
2829     [2]:
2830       1
2831   [3]:
2832     [1]:
2833       1
2834     [2]:
2835       d
2836"{//homogDistribution
2837  if (h == 0)
2838  {//trivial case where input is 0
2839    return(list(list(0,0)));
2840  }//trivial case where input is 0
2841  if (!isWeyl())
2842  {//Our basering is not the Weyl algebra
2843    ERROR("Ring was not the first Weyl algebra");
2844    return(list());
2845  }//Our basering is not the Weyl algebra
2846  if(nvars(basering)!=2)
2847  {//Our basering is the Weyl algebra, but not the first
2848    ERROR("Ring is not the first Weyl algebra");
2849    return(list());
2850  }//Our basering is the Weyl algebra, but not the first
2851  intvec ivm11 = intvec(-1,1);
2852  intvec iv1m1 = intvec(1,-1);
2853  poly tempH = h;
2854  poly minh;
2855  list result = list();
2856  int nextExpectedDegree = -deg(tempH,iv1m1);
2857  while (tempH != 0)
2858  {
2859    minh = jet(tempH,deg(tempH,iv1m1),iv1m1)-jet(tempH,deg(tempH,iv1m1)-1,iv1m1);
2860    while (deg(minh,ivm11)>nextExpectedDegree)
2861    {//filling empty homogeneous spaces with 0
2862      result = result + list(list(nextExpectedDegree,0));
2863      nextExpectedDegree = nextExpectedDegree +1;
2864    }//filling empty homogeneous spaces with 0
2865    result = result + list(list(deg(minh,ivm11),minh));
2866    tempH = tempH - minh;
2867    nextExpectedDegree = nextExpectedDegree +1;
2868  }
2869  return(result);
2870}//homogDistribution
2871
2872static proc countHomogParts(poly h)
2873"Counts the homogeneous parts of a given polynomial h"
2874{
2875  int i;
2876  list outPutHD = homogDistribution(h);
2877  int result = 0;
2878  for (i = 1; i <=size(outPutHD); i++)
2879  {
2880    if (outPutHD[i][2] != 0){result++;}
2881  }
2882  return(result);
2883}
2884
2885//////////////////////////////////////////////////
2886/////BRANDNEW!!!!////////////////////
2887//////////////////////////////////////////////////
2888
2889//==================================================
2890/*Singular has no way implemented to test polynomials
2891  for homogenity with respect to a weight vector.
2892  The following procedure does exactly this*/
2893static proc homogwithorder(poly h, intvec weights)
2894{//proc homogwithorder
2895  if(size(weights) != nvars(basering))
2896  {//The user does not know how many variables the current ring has
2897    return(0);
2898  }//The user does not know how many variables the current ring has
2899  int i;
2900  int dofp = deg(h,weights); //degree of polynomial
2901  for (i = 1; i<=size(h);i++)
2902  {
2903    if (deg(h[i],weights)!=dofp)
2904    {
2905      return(0);
2906    }
2907  }
2908  return(1);
2909}//proc homogwithorder
2910
2911//==================================================
2912//Testfac: Given a list with different factorizations of
2913// one polynomial, the following procedure checks
2914// whether they all refer to the same polynomial.
2915// If they do, the output will be a list, that contains
2916// the product of each factorization. If not, the empty
2917// list will be returned.
2918// If the optional argument # is given (i.e. the polynomial
2919// which is factorized by the elements of the given list),
2920// then we look, if the entries are factorizations of p
2921// and if not, a list with the products subtracted by p
2922// will be returned
2923proc testNCfac(list l, list #)
2924"USAGE: testNCfac(l[,p,b]); l is a list, p is an optional poly, b is 1 or 0
2925RETURN: Case 1: No optional argument. In this case the output is 1, if the
2926  entries in the given list represent the same polynomial or 0
2927  otherwise.
2928 Case 2: One optional argument p is given. In this case it returns 1,
2929  if all the entries in l are factorizations of p, otherwise 0.
2930 Case 3: Second optional b is given. In this case a list is returned
2931  containing the difference between the product of each entry in
2932  l and p.
2933ASSUME: basering is the first Weyl algebra, the entries of l are polynomials
2934PURPOSE: Checks whether a list of factorizations contains factorizations of
2935  the same element in the first Weyl algebra
2936THEORY: @code{testNCfac} multiplies out each factorization and checks whether
2937 each factorization was a factorization of the same element.
2938@* - if there is only a list given, the output will be 0, if it
2939     does not contain factorizations of the same element. Otherwise the output
2940     will be 1.
2941@* - if there is a polynomial in the second argument, then the procedure checks
2942     whether the given list contains factorizations of this polynomial. If it
2943     does, then the output depends on the third argument. If it is not given,
2944     the procedure will check whether the factorizations in the list
2945     l are associated to this polynomial and return either 1 or 0, respectively.
2946     If the third argument is given, the output will be a list with
2947     the length of the given one and in each entry is the product of one
2948     entry in l subtracted by the polynomial.
2949EXAMPLE: example testNCfac; shows examples
2950SEE ALSO: facFirstWeyl, facSubWeyl, facFirstShift
2951"{//proc testfac
2952  int p = printlevel - voice + 2;
2953  int i;
2954  string dbprintWhitespace = "";
2955  for (i = 1; i<=voice;i++)
2956  {dbprintWhitespace = dbprintWhitespace + " ";}
2957  dbprint(p,dbprintWhitespace + " Checking the input");
2958  if (size(l)==0)
2959  {//The empty list is given
2960    dbprint(p,dbprintWhitespace + " Given list was empty");
2961    return(list());
2962  }//The empty list is given
2963  if (size(#)>2)
2964  {//We want max. two optional arguments
2965    dbprint(p,dbprintWhitespace + " More than two optional arguments");
2966    return(list());
2967  }//We want max. two optional arguments
2968  dbprint(p,dbprintWhitespace + " Done");
2969  list result;
2970  int j;
2971  if (size(#)==0)
2972  {//No optional argument is given
2973    dbprint(p,dbprintWhitespace + " No optional arguments");
2974    int valid = 1;
2975    for (i = size(l);i>=1;i--)
2976    {//iterate over the elements of the given list
2977      if (size(result)>0)
2978      {
2979        if (product(l[i])!=result[size(l)-i])
2980        {
2981          valid = 0;
2982          break;
2983        }
2984      }
2985      result = insert(result, product(l[i]));
2986    }//iterate over the elements of the given list
2987    return(valid);
2988  }//No optional argument is given
2989  else
2990  {
2991    dbprint(p,dbprintWhitespace + " Optional arguments are given.");
2992    int valid = 1;
2993    for (i = size(l);i>=1;i--)
2994    {//iterate over the elements of the given list
2995      if (product(l[i])!=#[1])
2996      {
2997        valid = 0;
2998      }
2999      result = insert(result, product(l[i])-#[1]);
3000    }//iterate over the elements of the given list
3001    if(size(#)==2)
3002    {
3003      dbprint(p,dbprintWhitespace + " A third argument is given. Output is a list now.");
3004      return(result);
3005    }
3006    return(valid);
3007  }
3008}//proc testfac
3009example
3010{
3011  "EXAMPLE:";echo=2;
3012  ring r = 0,(x,y),dp;
3013  def R = nc_algebra(1,1);
3014  setring R;
3015  poly h = (x^2*y^2+1)*(x^2);
3016  def t1 = facFirstWeyl(h);
3017  //fist a correct list
3018  testNCfac(t1);
3019  //now a correct list with the factorized polynomial
3020  testNCfac(t1,h);
3021  //now we put in an incorrect list without a polynomial
3022  t1[3][3] = y;
3023  testNCfac(t1);
3024  // take h as additional input
3025  testNCfac(t1,h);
3026  // take h as additional input and output list of differences
3027  testNCfac(t1,h,1);
3028}
3029//==================================================
3030//Procedure facSubWeyl:
3031//This procedure serves the purpose to compute a
3032//factorization of a given polynomial in a ring, whose subring
3033//is the first Weyl algebra. The polynomial must only contain
3034//the two arguments, which are also given by the user.
3035
3036proc facSubWeyl(poly h, X, D)
3037"USAGE:  facSubWeyl(h,x,y); h, X, D polynomials
3038RETURN: list
3039ASSUME: X and D are variables of a basering, which satisfy DX = XD +1.
3040@* That is,  they generate the copy of the first Weyl algebra in a basering.
3041@* Moreover, h is a polynomial in X and D only.
3042PURPOSE: compute factorizations of the polynomial, which depends on X and D.
3043EXAMPLE: example facSubWeyl; shows examples
3044SEE ALSO: facFirstWeyl, testNCfac, facFirstShift
3045"{
3046  int p = printlevel - voice + 2;
3047  dbprint(p," Start initial Checks of the input.");
3048  // basering can be anything having a Weyl algebra as subalgebra
3049  def @r = basering;
3050  //We begin to check the input for assumptions
3051  // which are: X,D are vars of the basering,
3052  if ( (isVar(X)!=1) || (isVar(D)!=1) || (size(X)>1) || (size(D)>1) ||
3053       (leadcoef(X) != number(1)) || (leadcoef(D) != number(1)) )
3054  {
3055    ERROR("expected pure variables as generators of a subalgebra");
3056  }
3057  // Weyl algebra:
3058  poly w = D*X-X*D-1; // [D,X]=1
3059  poly u = D*X-X*D+1; // [X,D]=1
3060  if (u*w!=0)
3061  {
3062    // that is no combination gives Weyl
3063    ERROR("2nd and 3rd argument do not generate a Weyl algebra");
3064  }
3065  // one of two is correct
3066  int isReverted = 0; // Reverted Weyl if dx=xd-1 holds
3067  if (u==0)
3068  {
3069    isReverted = 1;
3070  }
3071  // else: do nothing
3072  // DONE with assumptions, Input successfully checked
3073  dbprint(p," Successful");
3074  intvec lexpofX = leadexp(X);
3075  intvec lexpofD = leadexp(D);
3076  int varnumX=1;
3077  int varnumD=1;
3078  while(lexpofX[varnumX] != 1)
3079  {
3080    varnumX++;
3081  }
3082  while(lexpofD[varnumD] != 1)
3083  {
3084    varnumD++;
3085  }
3086  /* VL : to add printlevel stuff */
3087  dbprint(p," Change positions of the two variables in the list, if needed");
3088  if (isReverted)
3089  {
3090    ring firstweyl = 0,(var(varnumD),var(varnumX)),dp;
3091    def Firstweyl = nc_algebra(1,1);
3092    setring Firstweyl;
3093    ideal M = 0:nvars(@r);
3094    M[varnumX]=var(2);
3095    M[varnumD]=var(1);
3096    map Q = @r,M;
3097    poly h= Q(h);
3098  }
3099  else
3100  { // that is unReverted
3101    ring firstweyl = 0,(var(varnumX),var(varnumD)),dp;
3102    def Firstweyl = nc_algebra(1,1);
3103    setring Firstweyl;
3104    poly h= imap(@r,h);
3105  }
3106  dbprint(p," Done!");
3107  list result = facFirstWeyl(h);
3108  setring @r;
3109  list result;
3110  if (isReverted)
3111  {
3112    // map swap back
3113    ideal M; M[1] = var(varnumD); M[2] = var(varnumX);
3114    map S = Firstweyl, M;
3115    result = S(result);
3116  }
3117  else
3118  {
3119    // that is unReverted
3120    result = imap(Firstweyl,result);
3121  }
3122  return(result);
3123}//proc facSubWeyl
3124example
3125{
3126  "EXAMPLE:";echo=2;
3127  ring r = 0,(x,y,z),dp;
3128  matrix D[3][3]; D[1,3]=-1;
3129  def R = nc_algebra(1,D); // x,z generate Weyl subalgebra
3130  setring R;
3131  poly h = (x^2*z^2+x)*x;
3132  list fact1 = facSubWeyl(h,x,z);
3133  // compare with facFirstWeyl:
3134  ring s = 0,(z,x),dp;
3135  def S = nc_algebra(1,1); setring S;
3136  poly h = (x^2*z^2+x)*x;
3137  list fact2 = facFirstWeyl(h);
3138  map F = R,x,0,z;
3139  list fact1 = F(fact1); // it is identical to list fact2
3140  testNCfac(fact1); // check the correctness again
3141}
3142//==================================================
3143
3144//==================================================
3145//************From here: Shift-Algebra**************
3146//==================================================
3147//==================================================*
3148//one factorization of a homogeneous polynomial
3149//in the first Shift Algebra
3150static proc homogfacFirstShift(poly h)
3151{//proc homogfacFirstShift
3152  int p=printlevel-voice+2; //for dbprint
3153  def r = basering;
3154  poly hath;
3155  intvec iv01 = intvec(0,1);
3156  int i; int j;
3157  string dbprintWhitespace = "";
3158  for (i = 1; i<=voice;i++)
3159  {dbprintWhitespace = dbprintWhitespace + " ";}
3160  if (!homogwithorder(h,iv01))
3161  {//The given polynomial is not homogeneous
3162    ERROR("The given polynomial is not homogeneous.");
3163    return(list());
3164  }//The given polynomial is not homogeneous
3165  if (h==0)
3166  {
3167    return(list(0));
3168  }
3169  list result;
3170  int m = deg(h,iv01);
3171  dbprint(p,dbprintWhitespace+" exclude the homogeneous part of deg. 0");
3172  if (m>0)
3173  {//The degree is not zero
3174    hath = lift(var(2)^m,h)[1,1];
3175    for (i = 1; i<=m;i++)
3176    {
3177      result = result + list(var(2));
3178    }
3179  }//The degree is not zero
3180  else
3181  {//The degree is zero
3182    hath = h;
3183  }//The degree is zero
3184  ring tempRing = 0,(x),dp;
3185  setring tempRing;
3186  map thetamap = r,x,1;
3187  poly hath = thetamap(hath);
3188  dbprint(p,dbprintWhitespace+" Factorize it using commutative factorization.");
3189  list azeroresult = factorize(hath);
3190  list azeroresult_return_form;
3191  for (i = 1; i<=size(azeroresult[1]);i++)
3192  {//rewrite the result of the commutative factorization
3193    for (j = 1; j <= azeroresult[2][i];j++)
3194    {
3195      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
3196    }
3197  }//rewrite the result of the commutative factorization
3198  setring(r);
3199  map finalmap = tempRing,var(1);
3200  list tempresult = finalmap(azeroresult_return_form);
3201  result = tempresult+result;
3202  return(result);
3203}//proc homogfacFirstShift
3204
3205//==================================================
3206//Computes all possible homogeneous factorizations
3207static proc homogfacFirstShift_all(poly h)
3208{//proc HomogfacFirstShiftAll
3209  int p=printlevel-voice+2; //for dbprint
3210  intvec iv11 = intvec(1,1);
3211  if (deg(h,iv11) <= 0 )
3212  {//h is a constant
3213    return(list(list(h)));
3214  }//h is a constant
3215  def r = basering;
3216  list one_hom_fac; //stands for one homogeneous factorization
3217  int i; int j; int k;
3218  string dbprintWhitespace = "";
3219  for (i = 1; i<=voice;i++)
3220  {dbprintWhitespace = dbprintWhitespace + " ";}
3221  int shiftcounter;
3222  //Compute again a homogeneous factorization
3223  dbprint(p,dbprintWhitespace+" Computing one homog. factorization of the polynomial");
3224  one_hom_fac = homogfacFirstShift(h);
3225  one_hom_fac = delete(one_hom_fac,1);
3226  if (size(one_hom_fac) == 0)
3227  {//there is no homogeneous factorization or the polynomial was not homogeneous
3228    return(list());
3229  }//there is no homogeneous factorization or the polynomial was not homogeneous
3230  dbprint(p,dbprintWhitespace+" Permuting the 0-homogeneous part with the s");
3231  list result = permpp(one_hom_fac);
3232  for (i = 1; i<=size(result);i++)
3233  {
3234    shiftcounter = 0;
3235    for (j = 1; j<=size(result[i]); j++)
3236    {
3237      if (result[i][j]==var(2))
3238      {
3239        shiftcounter++;
3240      }
3241      else
3242      {
3243        result[i][j] = subst(result[i][j], var(1), var(1)-shiftcounter);
3244      }
3245    }
3246    result[i] = insert(result[i],1);
3247  }
3248  dbprint(p,dbprintWhitespace+" Deleting double entries in the resulting list");
3249  result = delete_dublicates_noteval(result);
3250  return(result);
3251}//proc HomogfacFirstShiftAll
3252
3253//==================================================
3254//factorization of the first Shift Algebra
3255proc facFirstShift(poly h)
3256"USAGE: facFirstShift(h); h a polynomial in the first shift algebra
3257RETURN: list
3258PURPOSE: compute all factorizations of a polynomial in the first shift algebra
3259THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle
3260ASSUME: basering is the first shift algebra
3261NOTE: Every entry of the output list is a list with factors for one possible factorization.
3262EXAMPLE: example facFirstShift; shows examples
3263SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl
3264"{//facFirstShift
3265  int p = printlevel - voice + 2;
3266  int i;
3267  string dbprintWhitespace = "";
3268  for (i = 1; i<=voice;i++)
3269  {dbprintWhitespace = dbprintWhitespace + " ";}
3270  dbprint(p,dbprintWhitespace +" Checking the input.");
3271  if(nvars(basering)!=2)
3272  {//Our basering is the Shift algebra, but not the first
3273    ERROR("Basering is not the first shift algebra");
3274    return(list());
3275  }//Our basering is the Shift algebra, but not the first
3276  def r = basering;
3277  setring r;
3278  list LR = ringlist(r);
3279  number @n = leadcoef(LR[5][1,2]);
3280  poly @p = LR[6][1,2];
3281  if  ( @n!=number(1) )
3282  {
3283    ERROR("Basering is not the first shift algebra");
3284    return(list());
3285  }
3286  dbprint(p,dbprintWhitespace +" Done");
3287  list result = list();
3288  int j; int k; int l; //counter
3289  // create a ring with the ordering which makes shift algebra
3290  // graded
3291  // def r = basering; // done before
3292  ring tempRing = LR[1][1],(x,s),(a(0,1),Dp);
3293  def tempRingnc = nc_algebra(1,s);
3294  setring r;
3295  // information on relations
3296  if (@p == -var(1)) // reverted shift algebra
3297  {
3298    dbprint(p,dbprintWhitespace +" Reverted shift algebra. Swaping variables in Ringlist");
3299    setring(tempRingnc);
3300    map transf = r, var(2), var(1);
3301    setring(r);
3302    map transfback = tempRingnc, var(2),var(1);
3303    //    result = transfback(resulttemp);
3304  }
3305  else
3306  {
3307    if ( @p == var(2)) // usual shift algebra
3308    {
3309      setring(tempRingnc);
3310      map transf = r, var(1), var(2);
3311      //    result = facshift(h);
3312      setring(r);
3313      map transfback = tempRingnc, var(1),var(2);
3314    }
3315    else
3316    {
3317      ERROR("Basering is not the first shift algebra");
3318      return(list());
3319    }
3320  }
3321  // main calls
3322  setring(tempRingnc);
3323  dbprint(p,dbprintWhitespace +" Factorize the given polynomial with the subroutine sFacShift");
3324  list resulttemp = sFacShift(transf(h));
3325  dbprint(p,dbprintWhitespace +" Successful");
3326  setring(r);
3327  result = transfback(resulttemp);
3328  return( delete_dublicates_noteval(result) );
3329}//facFirstShift
3330example
3331{
3332  "EXAMPLE:";echo=2;
3333  ring R = 0,(x,s),dp;
3334  def r = nc_algebra(1,s);
3335  setring(r);
3336  poly h = (s^2*x+x)*s;
3337  facFirstShift(h);
3338}
3339
3340static proc sFacShift(poly h)
3341"
3342USAGE: A static procedure to factorize a polynomial in the first Shift algebra, where all the
3343       validity checks were made in advance.
3344INPUT:  A polynomial h in the first Shift Algebra.
3345OUTPUT: A list of different factorizations of h, where the factors are irreducible
3346ASSUMPTIONS:
3347 - The basering is the first Shift algebra and has n as first, and s as second variable, i.e. we
3348   have var(2)*var(1) = var(1)*var(2)+1
3349THEORY: If the given polynomial h is [0,1]-homogeneous, the routines for homogeneous factorizations
3350 are called. Otherwise we map the polynomial into the first Weyl algebra (the first shift
3351 algebra is a subring of the first Weyl algebra), and use facFirstWeyl to factorize it. Later
3352 we map the factors back, if possible.
3353"
3354{//proc sFacShift
3355  int p = printlevel - voice + 2;
3356  int i; int j ;
3357  string dbprintWhitespace = "";
3358  number commonCoefficient = content(h);
3359  for (i = 1; i<=voice;i++)
3360  {dbprintWhitespace = dbprintWhitespace + " ";}
3361  //Checking if given polynomial is homogeneous
3362  if(homogwithorder(h,intvec(0,1)))
3363  {//The given polynomial is [0,1]-homogeneous
3364    dbprint(p,dbprintWhitespace+"The polynomial is [0,1]-homogeneous. Returning the
3365homogeneous factorization");
3366    return(homogfacFirstShift_all(h));
3367  }//The given polynomial is [0,1]-homogeneous
3368
3369  //---------- Start of interesting part ----------
3370
3371  dbprint(p,dbprintWhitespace+"Mapping the polynomial h into the first Weyl algebra.");
3372  poly temph = h/commonCoefficient;
3373  def ourBaseRing = basering;
3374  ring tempWeylAlgebraComm = 0,(x,d),dp;
3375  def tempWeylAlgebra = nc_algebra(1,1);
3376  setring(tempWeylAlgebra);
3377  map shiftMap = ourBaseRing, x*d, d;
3378  poly h = shiftMap(temph);
3379  dbprint(p,dbprintWhitespace+"Successful! The polynomial in the Weyl algebra is "+string(h));
3380  dbprint(p,dbprintWhitespace+"Factorizing the polynomial in the first Weyl algebra");
3381  list factorizationInWeyl = facFirstWeyl(h);
3382  dbprint(p,dbprintWhitespace+"Successful! The factorization is given by:");
3383  dbprint(p,factorizationInWeyl);
3384  list validCombinations;
3385
3386  dbprint(p,dbprintWhitespace+"Now we will map this back to the shift algebra and filter
3387valid results");
3388  //-Now we map the results back to the shift algebra. But first, we need to combine them properly.
3389  for (i = 1; i<=size(factorizationInWeyl); i++)
3390  {//Deleting the first Coefficient factor
3391    factorizationInWeyl[i] = delete(factorizationInWeyl[i],1);
3392    validCombinations = validCombinations + combineNonnegative(factorizationInWeyl[i]);
3393  }//Deleting the first Coefficient factor
3394  if (size(validCombinations) == 0)
3395  {//There are no valid combinations, therefore we can directly say, that h is irreducible
3396    setring(ourBaseRing);
3397    return(list(list(commonCoefficient, h/commonCoefficient)));
3398  }//There are no valid combinations, therefore we can directly say, that h is irreducible
3399  validCombinations = delete_dublicates_noteval(validCombinations);
3400  setring(ourBaseRing);
3401  map backFromWeyl = tempWeylAlgebra, var(1),var(2);
3402  list validCombinations = backFromWeyl(validCombinations);
3403  for (i = 1; i<=size(validCombinations); i++)
3404  {
3405    for (j = 1; j<=size(validCombinations[i]);j++)
3406    {
3407      setring(tempWeylAlgebra);
3408      fromWeylToShiftPoly(validCombinations[i][j],ourBaseRing);
3409      validCombinations[i][j] = result;
3410      kill result;
3411      kill tempResult;
3412      kill zeroPoly;
3413      kill fromWeyl;
3414    }
3415  }
3416  for (i = 1; i<=size(validCombinations); i++)
3417  {//Adding the common factor in the first position of the list
3418    validCombinations[i] = insert(validCombinations[i],commonCoefficient);
3419  }//Adding the common factor in the first position of the list
3420  dbprint(dbprintWhitespace+"Done.");
3421  //mapping
3422  return(validCombinations);
3423}//proc sFacShift
3424
3425static proc combineNonnegative(list l)
3426"
3427USAGE: In sFacShift, when we want to map back the results of the factorization of the polynomial in
3428       the first Weyl algebra to the shift algebra. We need to recombine the factors such that
3429       we can map it back to the shift algebra without any problems.
3430INPUT: A list l containing one factorization of a polynomial in the first Weyl algebra. For example
3431       for the polynomial (1+x)*(1+x+d) we would have the list [1,x+1,x+d+1].
3432OUTPUT:If we can map every factor without a problem back to the shift algebra (i.e. if the smallest
3433 homogeneous summand of every factor is of nonnegative degree), a list containing the same
3434 list as given in the input is returned.
3435 If otherwise some factors cause problems, we consider every possible combination (i.e.
3436 products of the factors) and extract those where all factors have a smallest homogeneous
3437 summand of nonnegative degree.
3438ASSUMPTIONS:
3439 - Weyl algebra is given, and we have var(2)*var(1)=var(1)*var(2) +1
3440"
3441{//combineNonnegative
3442  int p = printlevel - voice + 2;
3443  int i;
3444  string dbprintWhitespace = "";
3445  for (i = 1; i<=voice;i++)
3446  {dbprintWhitespace = dbprintWhitespace + " ";}
3447  //First the easy case: all of the factors fulfill the condition of mapping to shift:
3448  dbprint(p,dbprintWhitespace+"Checking, if the given factors
3449can already be mapped without a problem.");
3450  int isValid = 1;
3451  for (i = 1; i<=size(l);i++)
3452  {//Checking for every entry if the condition is fulfilled.
3453    if (deg(l[i],intvec(1,-1))>0)
3454    {//Found one, where it is not fulfilled
3455      isValid = 0;
3456      break;
3457    }//Found one, where it is not fulfilled
3458  }//Checking for every entry if the condition is fulfilled.
3459  dbprint(p,dbprintWhitespace+"Done.");
3460  if (isValid)
3461  {//We can map every factor to the shift algebra and do not need to combine anything
3462    dbprint(p,dbprintWhitespace+"They can be mapped. Therefore we return them directly.");
3463    return(list(l));
3464  }//We can map every factor to the shift algebra and do not need to combine anything
3465  dbprint(p,dbprintWhitespace+"They cannot be mapped. Looking for valid combinations.");
3466  //Starting with the case, where l only consists of 1 or two elements.
3467  if(size(l)<=2)
3468  {//The case where we won't call the function a second time
3469    if (deg(product(l),intvec(1,-1))>0)
3470    {//No way of a valid combination
3471      return(list());
3472    }//No way of a valid combination
3473    else
3474    {//The product is the only possible and valid combination
3475      return(list(list(product(l))));
3476    }//The product is the only possible and valid combination
3477  }//The case where we won't call the function a second time
3478  //---------- Easy pre-stuff done. now we combine the factors.----------
3479  int pos;
3480  int j; int k;
3481  dbprint(p,dbprintWhitespace+"Making combinations of two.");
3482  list combinationsOfTwo = combinekfinlf(l,2);
3483  dbprint(p,dbprintWhitespace+"Done. Now checking, if there are valid ones in between.");
3484  list result;
3485  list validLHS;
3486  list validRHS;
3487  for (i = 1; i<=size(combinationsOfTwo); i++)
3488  {//go through all combinations and detect the valid ones
3489    if(deg(combinationsOfTwo[i][1],intvec(1,-1))>0 or deg(combinationsOfTwo[i][2],intvec(1,-1))>0)
3490    {//No chance, so no further treatment needed
3491      i++;
3492      continue;
3493    }//No chance, so no further treatment needed
3494    for (pos = 1; pos<=size(l);pos++)
3495    {//find the position where the combination splits
3496      if (product(l[1..pos]) == combinationsOfTwo[i][1])
3497      {//Found the position
3498        break;
3499      }//Found the position
3500    }//find the position where the combination splits
3501    dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " +
3502            string(list(l[1..pos])));
3503    validLHS = combineNonnegative(list(l[1..pos]));
3504    dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " +
3505            string(list(l[pos+1..size(l)])));
3506    validRHS = combineNonnegative(list(l[pos+1..size(l)]));
3507    for (j = 1; j<=size(validLHS); j++)
3508    {//Combining the left hand side valid combnations...
3509      for (k = 1; k<=size(validRHS); k++)
3510      {//... with the right hand side valid combinations
3511        result = insert(result, validLHS[j]+validRHS[k]);
3512      }//... with the right hand side valid combinations
3513    }//Combining the left hand side valid combnations...
3514  }//go through all combinations and detect the valid ones
3515  result = delete_dublicates_noteval(result);
3516  dbprint(p,dbprintWhitespace+"Done.");
3517  return(result);
3518}//combineNonnegative
3519
3520static proc fromWeylToShiftPoly(poly h, sAlgebra)
3521"
3522USAGE: Given a polynomial in the first Weyl algebra, this method returns it -- if possible --
3523       as an element in the first shift algebra, which is given in the method header.
3524INPUT: A polynomial h, and the first shift algebra as a ring
3525OUTPUT: The correct mapping in the shift Algebra
3526ASSUMPTIONS:
3527 - The lowest [-1,1]-homogeneous summand of h is of nonnegative degree
3528 - The shift algebra is given in the way that var(2)*var(1) = (var(1)+1)*var(2)
3529"
3530{//fromWeylToShiftPoly
3531  int p = printlevel - voice + 2;
3532  int i;
3533  string dbprintWhitespace = "";
3534  for (i = 1; i<=voice;i++)
3535  {dbprintWhitespace = dbprintWhitespace + " ";}
3536  if (deg(h,intvec(1,-1))>0)
3537  {//Wrong input polynomial
3538    ERROR("The lowest [-1,1] homogeneous summand of "+string(h)+" is of negative degree.");
3539  }//Wrong input polynomial
3540  def ourHomeBase = basering;
3541  list hDist = homogDistribution(h);
3542  setring(sAlgebra);
3543  poly result = 0;
3544  poly tempResult;
3545  poly zeroPoly;
3546  map fromWeyl = ourHomeBase, var(1), var(2);
3547  setring(ourHomeBase);
3548  poly zeroPoly;
3549  poly tempZeroPoly;
3550  int j; int k;
3551  int derDeg;
3552  for (i = 1; i<=size(hDist);i++)
3553  {
3554    derDeg = hDist[i][1];
3555    setring(sAlgebra);
3556    tempResult = 1;
3557    setring(ourHomeBase);
3558    zeroPoly = lift(d^derDeg, hDist[i][2])[1,1];
3559    for (j = 1; j<=size(zeroPoly); j++)
3560    {
3561      tempZeroPoly = zeroPoly[j];
3562      setring(sAlgebra);
3563      zeroPoly = fromWeyl(tempZeroPoly);
3564      tempResult = tempResult * leadcoef(zeroPoly);
3565      setring(ourHomeBase);
3566      for (k = 1; k<=deg(zeroPoly[j],intvec(0,1));k++)
3567      {
3568        setring(sAlgebra);
3569        tempResult = tempResult*(var(1)-(k-1));
3570        setring(ourHomeBase);
3571      }
3572      setring(sAlgebra);
3573      result = result + tempResult*var(2)^derDeg;
3574      tempResult = 1;
3575      setring(ourHomeBase);
3576    }
3577  }
3578  setring(sAlgebra);
3579  keepring(sAlgebra);
3580}//fromWeylToShiftPoly
3581
3582static proc refineFactList(list L)
3583{
3584  // assume: list L is an output of factorization proc
3585  // doing: remove doubled entries
3586  int s = size(L); int sm;
3587  int i,j,k,cnt;
3588  list M, U, A, B;
3589  A = L;
3590  k = 0;
3591  cnt  = 1;
3592  for (i=1; i<=s; i++)
3593  {
3594    if (size(A[i]) != 0)
3595    {
3596      M = A[i];
3597      //      "probing with"; M; i;
3598      B[cnt] = M; cnt++;
3599      for (j=i+1; j<=s; j++)
3600      {
3601        if ( isEqualList(M,A[j]) )
3602        {
3603          k++;
3604          // U consists of intvecs with equal pairs
3605          U[k] = intvec(i,j);
3606          A[j] = 0;
3607        }
3608      }
3609    }
3610  }
3611  kill A,U,M;
3612  return(B);
3613}
3614example
3615{
3616  "EXAMPLE:";echo=2;
3617  ring R = 0,(x,s),dp;
3618  def r = nc_algebra(1,1);
3619  setring(r);
3620  list l,m;
3621  l = list(1,s2+1,x,s,x+s);
3622  m = l,list(1,s,x,s,x),l;
3623  refineFactList(m);
3624}
3625
3626static proc isEqualList(list L, list M)
3627{
3628  // int boolean: 1=yes, 0 =no : test whether two lists are identical
3629  int s = size(L);
3630  if (size(M)!=s) { return(0); }
3631  int j=1;
3632  while ( (L[j]==M[j]) && (j<s) )
3633  {
3634    j++;
3635  }
3636  if (L[j]==M[j])
3637  {
3638    return(1);
3639  }
3640  return(0);
3641}
3642example
3643{
3644  "EXAMPLE:";echo=2;
3645  ring R = 0,(x,s),dp;
3646  def r = nc_algebra(1,1);
3647  setring(r);
3648  list l,m;
3649  l = list(1,s2+1,x,s,x+s);
3650  m = l;
3651  isEqualList(m,l);
3652}
3653
3654
3655//////////////////////////////////////////////////
3656// Q-WEYL-SECTION
3657//////////////////////////////////////////////////
3658
3659//==================================================
3660//A function to get the i'th triangular number
3661static proc triangNum(int n)
3662{
3663  if (n == 0)
3664  {
3665    return(0);
3666  }
3667  return (n*(n+1) div 2);
3668}
3669
3670//==================================================*
3671//one factorization of a homogeneous polynomial
3672//in the first Q Weyl Algebra
3673proc homogfacFirstQWeyl(poly h)
3674"USAGE: homogfacFirstQWeyl(h); h is a homogeneous polynomial in the
3675 first q-Weyl algebra with respect to the weight vector [-1,1]
3676RETURN: list
3677PURPOSE: Computes a factorization of a homogeneous polynomial h with
3678  respect to the weight vector [-1,1] in the first q-Weyl algebra
3679THEORY: @code{homogfacFirstQWeyl} returns a list with a factorization of the given,
3680 [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with
3681 k positive, the last k entries in the output list are the second
3682 variable. If k is positive, the last k entries will be x. The other
3683 entries will be irreducible polynomials of degree zero or 1 resp. -1.
3684SEE ALSO: homogfacFirstQWeyl_all
3685"{//proc homogfacFirstQWeyl
3686  int p = printlevel-voice+2;//for dbprint
3687  def r = basering;
3688  poly hath;
3689  int i; int j;
3690  string dbprintWhitespace = "";
3691  for (i = 1; i<=voice;i++)
3692  {dbprintWhitespace = dbprintWhitespace + " ";}
3693  intvec ivm11 = intvec(-1,1);
3694  if (!homogwithorder(h,ivm11))
3695  {//The given polynomial is not homogeneous
3696    ERROR("Given polynomial was not [-1,1]-homogeneous");
3697    return(list());
3698  }//The given polynomial is not homogeneous
3699  if (h==0)
3700  {
3701    return(list(0));
3702  }
3703  list result;
3704  int m = deg(h,ivm11);
3705  dbprint(p,dbprintWhitespace+" Splitting the polynomial in A_0 and A_k-Part");
3706  if (m!=0)
3707  {//The degree is not zero
3708    if (m <0)
3709    {//There are more x than y
3710      hath = lift(var(1)^(-m),h)[1,1];
3711      for (i = 1; i<=-m; i++)
3712      {
3713        result = result + list(var(1));
3714      }
3715    }//There are more x than y
3716    else
3717    {//There are more y than x
3718      hath = lift(var(2)^m,h)[1,1];
3719      for (i = 1; i<=m;i++)
3720      {
3721        result = result + list(var(2));
3722      }
3723    }//There are more y than x
3724  }//The degree is not zero
3725  else
3726  {//The degree is zero
3727    hath = h;
3728  }//The degree is zero
3729  dbprint(p," Done");
3730  //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1)
3731  list mons;
3732  dbprint(p," Putting the monomials in the A_0-part in a list.");
3733  for(i = 1; i<=size(hath);i++)
3734  {//Putting the monomials in a list
3735    mons = mons+list(hath[i]);
3736  }//Putting the monomials in a list
3737  dbprint(p," Done");
3738  dbprint(p," Mapping this monomials to K(q)[theta]");
3739  //Now, map to the commutative ring with theta:
3740  list tempRingList = ringlist(r);
3741  tempRingList[2] = insert(tempRingList[2],"theta",2); //New variable theta = x*d
3742  tempRingList = delete(tempRingList,5);
3743  tempRingList = delete(tempRingList,5); //The ring should now be commutative
3744  def tempRing = ring(tempRingList);
3745  setring tempRing;
3746  map thetamap = r,var(1),var(2);
3747  list mons = thetamap(mons);
3748  poly entry;
3749  poly tempSummand;
3750  for (i = 1; i<=size(mons);i++)
3751  {//transforming the monomials as monomials in theta
3752    entry = 1;//leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
3753    for (j = 0; j<leadexp(mons[i])[2];j++)
3754    {
3755      tempSummand = (par(1)^j-1)/(par(1)-1);
3756      entry = entry * theta-tempSummand*entry;
3757    }
3758    //entry;
3759    //leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
3760    mons[i] = entry*leadcoef(mons[i]) * par(1)^(-triangNum(leadexp(mons[i])[2]-1));
3761  }//transforming the monomials as monomials in theta
3762  dbprint(p," Done");
3763  dbprint(p," Factorize the A_0-Part in K[theta]");
3764  list azeroresult = factorize(sum(mons));
3765  dbprint(p," Successful");
3766  list azeroresult_return_form;
3767  for (i = 1; i<=size(azeroresult[1]);i++)
3768  {//rewrite the result of the commutative factorization
3769    for (j = 1; j <= azeroresult[2][i];j++)
3770    {
3771      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
3772    }
3773  }//rewrite the result of the commutative factorization
3774  dbprint(p," Mapping back to A_0.");
3775  setring(r);
3776  map finalmap = tempRing,var(1),var(2),var(1)*var(2);
3777  list tempresult = finalmap(azeroresult_return_form);
3778  dbprint(p,"Successful.");
3779  for (i = 1; i<=size(tempresult);i++)
3780  {//factorizations of theta resp. theta +1
3781    if(tempresult[i]==var(1)*var(2))
3782    {
3783      tempresult = insert(tempresult,var(1),i-1);
3784      i++;
3785      tempresult[i]=var(2);
3786    }
3787    if(tempresult[i]==var(2)*var(1))
3788    {
3789      tempresult = insert(tempresult,var(2),i-1);
3790      i++;
3791      tempresult[i]=var(1);
3792    }
3793  }//factorizations of theta resp. theta +1
3794  result = tempresult+result;
3795  //Correction of the result in the special q-Case:
3796  for (j = 2 ; j<= size(result);j++)
3797  {//Div the whole Term by the leading coefficient and multiply it to the first entry in result[i]
3798    result[1] = result[1] * leadcoef(result[j]);
3799    result[j] = 1/leadcoef(result[j]) * result[j];
3800  }//Div the whole Term by the leading coefficient and multiply it to the first entry in result[i]
3801  return(result);
3802}//proc homogfacFirstQWeyl
3803example
3804{
3805  "EXAMPLE:";echo=2;
3806  ring R = (0,q),(x,d),dp;
3807  def r = nc_algebra (q,1);
3808  setring(r);
3809  poly h = q^25*x^10*d^10+q^16*(q^4+q^3+q^2+q+1)^2*x^9*d^9+
3810    q^9*(q^13+3*q^12+7*q^11+13*q^10+20*q^9+26*q^8+30*q^7+
3811    31*q^6+26*q^5+20*q^4+13*q^3+7*q^2+3*q+1)*x^8*d^8+
3812    q^4*(q^9+2*q^8+4*q^7+6*q^6+7*q^5+8*q^4+6*q^3+
3813     4*q^2+2q+1)*(q^4+q^3+q^2+q+1)*(q^2+q+1)*x^7*d^7+
3814    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+
3815    (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+
3816    6*x^3*d^3+24;
3817  homogfacFirstQWeyl(h);
3818}
3819
3820//==================================================
3821//Computes all possible homogeneous factorizations for an element in the first Q-Weyl Algebra
3822proc homogfacFirstQWeyl_all(poly h)
3823"USAGE: homogfacFirstQWeyl_all(h); h is a homogeneous polynomial in the first q-Weyl algebra
3824 with respect to the weight vector [-1,1]
3825RETURN: list
3826PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect
3827  to the weight vector [-1,1] in the first q-Weyl algebra
3828THEORY: @code{homogfacFirstQWeyl} returns a list with all factorization of the given,
3829 homogeneous polynomial. It uses the output of homogfacFirstQWeyl and permutes
3830 its entries with respect to the commutation rule. Furthermore, if a
3831 factor of degree zero is irreducible in K[  heta], but reducible in
3832 the first q-Weyl algebra, the permutations of this element with the other
3833 entries will also be computed.
3834SEE ALSO: homogfacFirstQWeyl
3835"{//proc HomogfacFirstQWeylAll
3836  int p=printlevel-voice+2;//for dbprint
3837  intvec iv11= intvec(1,1);
3838  if (deg(h,iv11) <= 0 )
3839  {//h is a constant
3840    dbprint(p,"Given polynomial was not homogeneous");
3841    return(list(list(h)));
3842  }//h is a constant
3843  def r = basering;
3844  list one_hom_fac; //stands for one homogeneous factorization
3845  int i; int j; int k;
3846  intvec ivm11 = intvec(-1,1);
3847  dbprint(p," Calculate one homogeneous factorization using homogfacFirstQWeyl");
3848  //Compute again a homogeneous factorization
3849  one_hom_fac = homogfacFirstQWeyl(h);
3850  dbprint(p,"Successful");
3851  if (size(one_hom_fac) == 0)
3852  {//there is no homogeneous factorization or the polynomial was not homogeneous
3853    return(list());
3854  }//there is no homogeneous factorization or the polynomial was not homogeneous
3855  //divide list in A0-Part and a list of x's resp. y's
3856  list list_not_azero = list();
3857  list list_azero;
3858  list k_factor;
3859  int is_list_not_azero_empty = 1;
3860  int is_list_azero_empty = 1;
3861  k_factor = list(one_hom_fac[1]);
3862  if (absValue(deg(h,ivm11))<size(one_hom_fac)-1)
3863  {//There is a nontrivial A_0-part
3864    list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,ivm11)))];
3865    is_list_azero_empty = 0;
3866  }//There is a nontrivial A_0 part
3867  dbprint(p," Combine x,y to xy in the factorization again.");
3868  for (i = 1; i<=size(list_azero)-1;i++)
3869  {//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
3870    if (list_azero[i] == var(1))
3871    {
3872      if (list_azero[i+1]==var(2))
3873      {
3874        list_azero[i] = var(1)*var(2);
3875        list_azero = delete(list_azero,i+1);
3876      }
3877    }
3878    if (list_azero[i] == var(2))
3879    {
3880      if (list_azero[i+1]==var(1))
3881      {
3882        list_azero[i] = var(2)*var(1);
3883        list_azero = delete(list_azero,i+1);
3884      }
3885    }
3886  }//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
3887  dbprint(p," Done");
3888  if(deg(h,ivm11)!=0)
3889  {//list_not_azero is not empty
3890    list_not_azero =
3891      one_hom_fac[(size(one_hom_fac)-absValue(deg(h,ivm11))+1)..size(one_hom_fac)];
3892    is_list_not_azero_empty = 0;
3893  }//list_not_azero is not empty
3894  //Map list_azero in K[theta]
3895  dbprint(p," Map list_azero to K[theta]");
3896  //Now, map to the commutative ring with theta:
3897  list tempRingList = ringlist(r);
3898  tempRingList[2] = insert(tempRingList[2],"theta",2); //New variable theta = x*d
3899  tempRingList = delete(tempRingList,5);
3900  tempRingList = delete(tempRingList,5); //The ring should now be commutative
3901  def tempRing = ring(tempRingList);
3902  setring(tempRing);
3903  poly entry;
3904  map thetamap = r,var(1),var(2);
3905  if(!is_list_not_azero_empty)
3906  {//Mapping in Singular is only possible, if the list before
3907    //contained at least one element of the other ring
3908    list list_not_azero = thetamap(list_not_azero);
3909  }//Mapping in Singular is only possible, if the list before
3910  //contained at least one element of the other ring
3911  if(!is_list_azero_empty)
3912  {//Mapping in Singular is only possible, if the list before
3913    //contained at least one element of the other ring
3914    list list_azero= thetamap(list_azero);
3915  }//Mapping in Singular is only possible, if the list before
3916  //contained at least one element of the other ring
3917  list k_factor = thetamap(k_factor);
3918  list tempmons;
3919  dbprint(p," Done");
3920  for(i = 1; i<=size(list_azero);i++)
3921  {//rewrite the polynomials in A1 as polynomials in K[theta]
3922    tempmons = list();
3923    for (j = 1; j<=size(list_azero[i]);j++)
3924    {
3925      tempmons = tempmons + list(list_azero[i][j]);
3926    }
3927    for (j = 1 ; j<=size(tempmons);j++)
3928    {
3929      //entry = leadcoef(tempmons[j]);
3930      entry = leadcoef(tempmons[j]) * par(1)^(-triangNum(leadexp(tempmons[j])[2]-1));
3931      for (k = 0; k < leadexp(tempmons[j])[2];k++)
3932      {
3933        entry = entry*(theta-(par(1)^k-1)/(par(1)-1));
3934      }
3935      tempmons[j] = entry;
3936    }
3937    list_azero[i] = sum(tempmons);
3938  }//rewrite the polynomials in A1 as polynomials in K[theta]
3939  //Compute all permutations of the A0-part
3940  dbprint(p," Compute all permutations of the A_0-part with the first resp. the snd. variable");
3941  list result;
3942  int shift_sign;
3943  int shift;
3944  poly shiftvar;
3945  if (size(list_not_azero)!=0)
3946  {//Compute all possibilities to permute the x's resp. the y's in the list
3947    if (list_not_azero[1] == var(1))
3948    {//h had a negative weighted degree
3949      shift_sign = 1;
3950      shiftvar = var(1);
3951    }//h had a negative weighted degree
3952    else
3953    {//h had a positive weighted degree
3954      shift_sign = -1;
3955      shiftvar = var(2);
3956    }//h had a positive weighted degree
3957    result = permpp(list_azero + list_not_azero);
3958    for (i = 1; i<= size(result); i++)
3959    {//adjust the a_0-parts
3960      shift = 0;
3961      for (j=1; j<=size(result[i]);j++)
3962      {
3963        if (result[i][j]==shiftvar)
3964        {
3965          shift = shift + shift_sign;
3966        }
3967        else
3968        {
3969          if (shift < 0)
3970          {//We have two distict formulas for x and y. In this case use formula for y
3971            if (shift == -1)
3972            {
3973              result[i][j] = subst(result[i][j],theta,1/par(1)*(theta - 1));
3974            }
3975            else
3976            {
3977              result[i][j] =
3978                subst(result[i][j],
3979                      theta,
3980                      1/par(1)*((theta - 1)/par(1)^(absValue(shift)-1)
3981                                - (par(1)^(shift +2)-par(1))/(1-par(1))));
3982            }
3983          }//We have two distict formulas for x and y. In this case use formula for y
3984          if (shift > 0)
3985          {//We have two distict formulas for x and y. In this case use formula for x
3986            if (shift == 1)
3987            {
3988              result[i][j] = subst(result[i][j],theta,par(1)*theta + 1);
3989            }
3990            else
3991            {
3992              result[i][j] =
3993                subst(result[i][j],
3994                      theta,par(1)^shift*theta+(par(1)^shift-1)/(par(1)-1));
3995            }
3996          }//We have two distict formulas for x and y. In this case use formula for x
3997        }
3998      }
3999    }//adjust the a_0-parts
4000  }//Compute all possibilities to permute the x's resp. the y's in the list
4001  else
4002  {//The result is just all the permutations of the a_0-part
4003    result = permpp(list_azero);
4004  }//The result is just all the permutations of the a_0 part
4005  if (size(result)==0)
4006  {
4007    return(result);
4008  }
4009  dbprint(p," Done");
4010  dbprint(p," Searching for theta resp. theta + 1 in the list and factorize them");
4011  //Now we are going deeper and search for theta resp. theta + 1, substitute
4012  //them by xy resp. yx and go on permuting
4013  int found_theta;
4014  int thetapos;
4015  list leftpart;
4016  list rightpart;
4017  list lparts;
4018  list rparts;
4019  list tempadd;
4020  for (i = 1; i<=size(result) ; i++)
4021  {//checking every entry of result for theta or theta +1
4022    found_theta = 0;
4023    for(j=1;j<=size(result[i]);j++)
4024    {
4025      if (result[i][j]==theta)
4026      {//the jth entry is theta and can be written as x*y
4027        thetapos = j;
4028        result[i]= insert(result[i],var(1),j-1);
4029        j++;
4030        result[i][j] = var(2);
4031        found_theta = 1;
4032        break;
4033      }//the jth entry is theta and can be written as x*y
4034      if(result[i][j] == par(1)*theta +1)
4035      {
4036        thetapos = j;
4037        result[i] = insert(result[i],var(2),j-1);
4038        j++;
4039        result[i][j] = var(1);
4040        found_theta = 1;
4041        break;
4042      }
4043    }
4044    if (found_theta)
4045    {//One entry was theta resp. theta +1
4046      leftpart = result[i];
4047      leftpart = leftpart[1..thetapos];
4048      rightpart = result[i];
4049      rightpart = rightpart[(thetapos+1)..size(rightpart)];
4050      lparts = list(leftpart);
4051      rparts = list(rightpart);
4052      //first deal with the left part
4053      if (leftpart[thetapos] == var(1))
4054      {
4055        shift_sign = 1;
4056        shiftvar = var(1);
4057      }
4058      else
4059      {
4060        shift_sign = -1;
4061        shiftvar = var(2);
4062      }
4063      for (j = size(leftpart); j>1;j--)
4064      {//drip x resp. y
4065        if (leftpart[j-1]==shiftvar)
4066        {//commutative
4067          j--;
4068          continue;
4069        }//commutative
4070        if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
4071        {//stop here
4072          break;
4073        }//stop here
4074        //Here, we can only have a a0- part
4075        if (shift_sign<0)
4076        {
4077          leftpart[j] = subst(leftpart[j-1],theta, 1/par(1)*(theta +shift_sign));
4078        }
4079        if (shift_sign>0)
4080        {
4081          leftpart[j] = subst(leftpart[j-1],theta, par(1)*theta + shift_sign);
4082        }
4083        leftpart[j-1] = shiftvar;
4084        lparts = lparts + list(leftpart);
4085      }//drip x resp. y
4086      //and now deal with the right part
4087      if (rightpart[1] == var(1))
4088      {
4089        shift_sign = 1;
4090        shiftvar = var(1);
4091      }
4092      else
4093      {
4094        shift_sign = -1;
4095        shiftvar = var(2);
4096      }
4097      for (j = 1 ; j < size(rightpart); j++)
4098      {
4099        if (rightpart[j+1] == shiftvar)
4100        {
4101          j++;
4102          continue;
4103        }
4104        if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
4105        {
4106          break;
4107        }
4108        if (shift_sign<0)
4109        {
4110          rightpart[j] = subst(rightpart[j+1], theta, par(1)*theta - shift_sign);
4111        }
4112        if (shift_sign>0)
4113        {
4114          rightpart[j] = subst(rightpart[j+1], theta, 1/par(1)*(theta - shift_sign));
4115        }
4116        rightpart[j+1] = shiftvar;
4117        rparts = rparts + list(rightpart);
4118      }
4119      //And now, we put all possibilities together
4120      tempadd = list();
4121      for (j = 1; j<=size(lparts); j++)
4122      {
4123        for (k = 1; k<=size(rparts);k++)
4124        {
4125          tempadd = tempadd + list(lparts[j]+rparts[k]);
4126        }
4127      }
4128      tempadd = delete(tempadd,1); // The first entry is already in the list
4129      result = result + tempadd;
4130      continue; //We can may be not be done already with the ith entry
4131    }//One entry was theta resp. theta +1
4132  }//checking every entry of result for theta or theta +1
4133  dbprint(p," Done");
4134  //map back to the basering
4135  dbprint(p," Mapping back everything to the basering");
4136  setring(r);
4137  map finalmap = tempRing, var(1), var(2),var(1)*var(2);
4138  list result = finalmap(result);
4139  for (i=1; i<=size(result);i++)
4140  {//adding the K factor
4141    result[i] = k_factor + result[i];
4142  }//adding the k-factor
4143  dbprint(p," Done");
4144  dbprint(p," Delete double entries in the list.");
4145  result = delete_dublicates_noteval(result);
4146  dbprint(p," Done");
4147  return(result);
4148}//proc HomogfacFirstQWeylAll
4149example
4150{
4151  "EXAMPLE:";echo=2;
4152  ring R = (0,q),(x,d),dp;
4153  def r = nc_algebra (q,1);
4154  setring(r);
4155  poly h = q^25*x^10*d^10+q^16*(q^4+q^3+q^2+q+1)^2*x^9*d^9+
4156    q^9*(q^13+3*q^12+7*q^11+13*q^10+20*q^9+26*q^8+30*q^7+
4157    31*q^6+26*q^5+20*q^4+13*q^3+7*q^2+3*q+1)*x^8*d^8+
4158    q^4*(q^9+2*q^8+4*q^7+6*q^6+7*q^5+8*q^4+6*q^3+
4159    4*q^2+2q+1)*(q^4+q^3+q^2+q+1)*(q^2+q+1)*x^7*d^7+
4160    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+
4161    (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+
4162    6*x^3*d^3+24;
4163  homogfacFirstQWeyl_all(h);
4164}
4165
4166//TODO: FirstQWeyl check the parameters...
4167
4168//==================================================
4169// EASY EXAMPLES FOR WEYL ALGEBRA
4170//==================================================
4171/*
4172  Easy and fast example polynomials where one can find factorizations: K<x,d |dx=xd+1>
4173  (x^2+d)*(x^2+d);
4174  (x^2+x)*(x^2+d);
4175  (x^3+x+1)*(x^4+d*x+2);
4176  (x^2*d+d)*(d+x*d);
4177  d^3+x*d^3+2*d^2+2*(x+1)*d^2+d+(x+2)*d; //Example 5 Grigoriev-Schwarz.
4178  (d+1)*(d+1)*(d+x*d); //Landau Example projected to the first dimension.
4179*/
4180
4181//==================================================
4182//Some Bugs(fixed)/hard examples from Martin Lee:
4183//==================================================
4184// ex1, ex2
4185/*
4186ring s = 0,(x,d),Ws(-1,1);
4187def S = nc_algebra(1,1); setring S;
4188poly a = 10x5d4+26x4d5+47x5d2-97x4d3; //Not so hard any more... Done in around 4 minutes
4189def l= facFirstWeyl (a); l;
4190kill l;
4191poly b = -5328x8d5-5328x7d6+720x9d2+720x8d3-16976x7d4-38880x6d5
4192-5184x7d3-5184x6d4-3774x5d5+2080x8d+5760x7d2-6144x6d3-59616x5d4
4193+3108x3d6-4098x6d2-25704x5d3-21186x4d4+8640x6d-17916x4d3+22680x2d5
4194+2040x5d-4848x4d2-9792x3d3+3024x2d4-10704x3d2-3519x2d3+34776xd4
4195+12096xd3+2898d4-5040x2d+8064d3+6048d2; //Still very hard... But it seems to be only because of the
4196//combinatorial explosion
4197def l= facFirstWeyl (b); l;
4198
4199// ex3: there was difference in answers => fixed
4200LIB "ncfactor.lib";
4201ring r = 0,(x,y,z),dp;
4202matrix D[3][3]; D[1,3]=-1;
4203def R = nc_algebra(1,D);
4204setring R;
4205poly g= 7*z4*x+62*z3+26*z;
4206def l1= facSubWeyl (g, x, z);
4207l1;
4208//---- other ring
4209ring s = 0,(x,z),dp;
4210def S = nc_algebra(1,-1); setring S;
4211poly g= 7*z4*x+62*z3+26*z;
4212def l2= facFirstWeyl (g);
4213l2;
4214map F = R,x,0,z;
4215list l1 = F(l1);
4216l1;
4217//---- so the answers look different, check them!
4218testNCfac(l2); // ok
4219testNCfac(l1); // was not ok, but now it's been fixed!!!
4220
4221// selbst D und X so vertauschen dass sie erfuellt ist : ist gemacht
4222
4223*/
4224
4225/*
4226// bug from M Lee
4227LIB "ncfactor.lib";
4228ring s = 0,(z,x),dp;
4229def S = nc_algebra(1,1); setring S;
4230poly f= -60z4x2-54z4-56zx3-59z2x-64;
4231def l= facFirstWeyl (f);
4232l; // before: empty list; after fix: 1 entry, f is irreducible
4233poly g = 75z3x2+92z3+24;
4234def l= facFirstWeyl (g);
4235l; //before: empty list, now: correct
4236*/
4237
4238/* more things from Martin Lee; fixed
4239ring R = 0,(x,s),dp;
4240def r = nc_algebra(1,s);
4241setring(r);
4242poly h = (s2*x+x)*s;
4243h= h* (x+s);
4244def l= facFirstShift(h);
4245l; // contained doubled entries: not anymore, fixed!
4246
4247ring R = 0,(x,s),dp;
4248def r = nc_algebra(1,-1);
4249setring(r);
4250poly h = (s2*x+x)*s;
4251h= h* (x+s);
4252def l= facFirstWeyl(h);
4253l; // contained doubled entries: not anymore, fixed!
4254
4255*/
4256
4257//======================================================================
4258//Examples from TestSuite that are terminating in a reasonable time.
4259//======================================================================
4260
4261//Counter example for old Algorithm, but now working:
4262/*
4263ring R = 0,(x,d),dp;
4264def r = nc_algebra(1,1);
4265setring(r);
4266LIB "ncfactor.lib";
4267poly h = (1+x^2*d)^4;
4268list lsng = facFirstWeyl(h);
4269print(lsng);
4270*/
4271
4272//Example 2.7. from Master thesis
4273/*
4274ring R = 0,(x,d),dp;
4275def r = nc_algebra(1,1);
4276setring(r);
4277LIB "ncfactor.lib";
4278poly h = (xdd + xd+1+ (xd+5)*x)*(((x*d)^2+1)*d + xd+3+ (xd+7)*x);
4279list lsng = facFirstWeyl(h);
4280print(lsng);
4281 */
4282
4283//Example with high combinatorial income
4284/*
4285ring R = 0,(x,d),dp;
4286def r = nc_algebra(1,1);
4287setring(r);
4288LIB "ncfactor.lib";
4289poly 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);
4290list lsng = facFirstWeyl(h);
4291print(lsng);
4292 */
4293
4294//Once a bug, now working
4295/*
4296ring R = 0,(x,d),dp;
4297def r = nc_algebra(1,1);
4298setring(r);
4299LIB "ncfactor.lib";
4300poly h = (x^2*d^2+x)*(x+1);
4301list lsng = facFirstWeyl(h);
4302print(lsng);
4303*/
4304
4305//Another one of that kind
4306/*
4307ring R = 0,(x,d),dp;
4308def r = nc_algebra(1,1);
4309setring(r);
4310LIB "ncfactor.lib";
4311poly h = (x*d*d + (x*d)^5 +x)*((x*d+1)*d-(x*d-1)^5+x);
4312list lsng = facFirstWeyl(h);
4313print(lsng);
4314 */
4315
4316//Example of Victor for Shift Algebra
4317/*
4318ring s = 0,(n,Sn),dp;
4319def S = nc_algebra(1,Sn); setring S;
4320LIB "ncfactor.lib";
4321list lsng = facFirstShift(n^2*Sn^2+3*n*Sn^2-n^2+2*Sn^2-3*n-2);
4322print(lsng);
4323 */
4324
4325//Interesting example, as there are actually also some complex solutions to it:
4326/*
4327ring R = 0,(x,d),dp;
4328def r = nc_algebra(1,1);
4329setring(r);
4330LIB "/Users/albertheinle/Studium/forschung/ncfactor/versionen/ncfactor.lib";
4331poly h =(x^3+x+1)*(x^4+d*x+2);//Example for finitely many, but more than one solution in between.
4332list lsng = facFirstWeyl(h);
4333print(lsng);
4334 */
4335
4336//Another one of that kind:
4337/*
4338ring R = 0,(x,d),dp;
4339def r = nc_algebra(1,1);
4340setring(r);
4341LIB "ncfactor.lib";
4342poly h =(x^2+d)*(x^2+d);//Example for finitely many, but more than one solution in between.
4343list lsng = facFirstWeyl(h);
4344print(lsng);
4345*/
4346
4347//Example by W. Koepf:
4348/*
4349ring R = 0,(x,d),dp;
4350def r = nc_algebra(1,1);
4351setring(r);
4352LIB "ncfactor.lib";
4353poly h = (x^4-1)*x*d^2+(1+7*x^4)*d+8*x^3;
4354list lsng = facFirstWeyl(h);
4355print(lsng);
4356 */
4357
4358//Shift Example from W. Koepf
4359/*
4360ring R = 0,(n,s),dp;
4361def r = nc_algebra(1,s);
4362setring(r);
4363LIB "ncfactor.lib";
4364poly h = n*(n+1)*s^2-2*n*(n+100)*s+(n+99)*(n+100);
4365list lsng = facFirstShift(h);
4366print(lsng);
4367 */
4368
4369//Tsai Example... Once hard, now easy...
4370/*
4371ring R = 0,(x,d),dp;
4372def r = nc_algebra(1,1);
4373setring(r);
4374LIB "ncfactor.lib";
4375poly 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);
4376list lsng =facFirstWeyl(h);
4377print(lsng);
4378 */
4379
4380//======================================================================
4381// Hard examples not yet calculatable in feasible amount of time
4382//======================================================================
4383
4384//Also a counterexample for REDUCE. Very long Groebner basis computation in between.
4385/*
4386ring R = 0,(x,d),dp;
4387def r = nc_algebra(1,1);
4388setring(r);
4389LIB "ncfactor.lib";
4390poly h = (d^4+x^2+dx+x)*(d^2+x^4+xd+d);
4391list lsng = facFirstWeyl(h);
4392print(lsng);
4393*/
4394
4395//Example from the Mainz-Group
4396/*
4397ring R = 0,(x,d),dp;
4398def r = nc_algebra(1,1);
4399setring(r);
4400poly 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+
4401(-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+
4402(-144/35*(x+1)*(1225*x^5+11025*x^4+37485*x^3+61335*x^2+50138*x+16584)-6912/35*(x+2)*
4403(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)*
4404(x+2)^2)*d^3;
4405LIB "ncfactor.lib";
4406printlevel = 5;
4407facFirstWeyl(dop6);
4408$;*/
4409
4410//Another Mainz Example:
4411/*
4412LIB "ncfactor.lib";
4413ring R = 0,(x,d),dp;
4414def r = nc_algebra(1,1);
4415setring(r);
4416poly dopp = 82547*x^4*d^4+60237*x^3*d^3+26772*x^5*d^5+2231*x^6*d^6+x*(1140138*
4417x^2*d^2-55872*x*d-3959658*x^3*d^3-8381805*x^4*d^4-3089576*x^5*d^5-274786*
4418x^6*d^6)+x^2*(-16658622*x*d-83427714*x^2*d^2-19715033*x^3*d^3+78915395*x^4
4419*d^4+35337930*x^5*d^5+3354194*x^6*d^6)+x^3*(-99752472-1164881352*x*d+
44204408536996*x^2*d^2+11774185985*x^3*d^3+5262196786*x^4*d^4+1046030561/2*x^5*
4421d^5-10564451/2*x^6*d^6)+x^4*(-1925782272+21995375398*x*d+123415803356*x^2*
4422d^2+302465300831/2*x^3*d^3+34140803907/2*x^4*d^4-15535653409*x^5*d^5-\
44232277687768*x^6*d^6)+x^5*(71273525520+691398212366*x*d+901772633569*x^2*d^2+
44242281275427069*x^3*d^3+2944352819911/2*x^4*d^4+836872370039/4*x^5*d^5+
44259066399237/4*x^6*d^6)+x^6*(2365174430376+9596715855542*x*d+29459572469704*x^
44262*d^2+92502197003786*x^3*d^3+65712473180525*x^4*d^4+13829360193674*x^5*d^5
4427+3231449477251/4*x^6*d^6)+x^7*(26771079436836+117709870166226*x*d+
4428821686455179082*x^2*d^2+1803972139232179*x^3*d^3+1083654460691481*x^4*d^4+
4429858903621851785/4*x^5*d^5+50096565802957/4*x^6*d^6)+x^8*(179341727601960+
44302144653944040630*x*d+13123246960284302*x^2*d^2+41138357917778169/2*x^3*d^3+
443120605819587976401/2*x^4*d^4+3677396642905423/2*x^5*d^5+402688260229369/4*x^6
4432*d^6)+x^9*(2579190935961288+43587063726809764*x*d+157045086382352387*x^2*d^
44332+172175668477370223*x^3*d^3+138636285385875407/2*x^4*d^4+10707836398626232*
4434x^5*d^5+529435530567584*x^6*d^6)+x^10*(41501953525903392+558336731465626084*
4435x*d+1407267553543222268*x^2*d^2+1153046693323226808*x^3*d^3+
4436372331468563656085*x^4*d^4+48654019090240214*x^5*d^5+2114661191282167*x^6*d
4437^6)+x^11*(364526077273381884+4158060401095928464*x*d+8646807662899324262*x^2*
4438d^2+5914675753405705400*x^3*d^3+1631934058875116005*x^4*d^4+
4439187371894330537204*x^5*d^5+7366806367019734*x^6*d^6)+x^12*(
44401759850321214603648+18265471270535733520*x*d+34201910114871110912*x^2*d^2+
444121265221434709398152*x^3*d^3+5437363546219595036*x^4*d^4+594029113431041060*
4442x^5*d^5+22881659624561644*x^6*d^6)+x^13*(4648382639403200688+
444345699084277107816096*x*d+81049061578449009384*x^2*d^2+48858488665016574368*x
4444^3*d^3+12515362110098721444*x^4*d^4+1412152747420021048*x^5*d^5+
444557196947123984972*x^6*d^6)+x^14*(5459369397960020544+55837825300341621824*x*
4446d+105671876924055409696*x^2*d^2+71551727420848766624*x^3*d^3+
444721094786205096577808*x^4*d^4+2695663190297032192*x^5*d^5+118791751565613264*
4448x^6*d^6)+x^15*(1023333653580043776+47171127937488813824*x*d+
4449157258351906685700352*x^2*d^2+145765192195300531840*x^3*d^3+
445049876215785510342176*x^4*d^4+6647374188802036864*x^5*d^5+287310278455067312*
4451x^6*d^6)+x^16*(11960091747366236160+250326608568269289472*x*d+
4452677587171115580981248*x^2*d^2+538246374825683603456*x^3*d^3+
4453161380433451548754048*x^4*d^4+19149099315354950144*x^5*d^5+
4454746433247985092544*x^6*d^6)+x^17*(42246252365448668160+657220532737851248640*
4455x*d+1531751689216283911680*x^2*d^2+1090829514212206064640*x^3*d^3+
4456299280728709430851840*x^4*d^4+32932767387222323200*x^5*d^5+
44571202281367574179840*x^6*d^6)+x^18*(6239106101942784000+320638742839606579200*
4458x*d+873857213570556364800*x^2*d^2+645649080101933721600*x^3*d^3+
4459177008238160627276800*x^4*d^4+19165088507111475200*x^5*d^5+
4460683600826675660800*x^6*d^6)+x^19*(-60440251454613504000-476055211197689856000
4461*x*d-733497382597635072000*x^2*d^2-386038662982742016000*x^3*d^3-\
446283361486778142976000*x^4*d^4-7524999543181824000*x^5*d^5-232189492987008000*
4463x^6*d^6)+x^20*(1578562930483200000+12628503443865600000*x*d+
446419732036631040000000*x^2*d^2+10523752869888000000*x^3*d^3+
44652302070940288000000*x^4*d^4+210475057397760000*x^5*d^5+6577345543680000*x^6*
4466d^6);
4467printlevel = 3;
4468facFirstWeyl(dopp);
4469*/
4470
4471
4472
4473//Hard Example by Viktor:
4474/*
4475  ring r = 0,(x,d), (dp);
4476def R = nc_algebra(1,1);
4477setring R;
4478LIB "/Users/albertheinle/Studium/forschung/ncfactor/versionen/ncfactor.lib";
4479poly t = x; poly D =d;
4480poly 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
4481-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
4482-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
4483+12*t^3*D+2*t*D+4*t^3-8;
4484list lsng = facFirstWeyl(p);
4485print(lsng);
4486*/
Note: See TracBrowser for help on using the repository browser.