//////////////////////////////////////////////////////////
version="version ncfactor.lib 4.0.0.0 Jun_2013 "; // $Id$
category="Noncommutative";
info="
LIBRARY: ncfactor.lib  Tools for factorization in some noncommutative algebras
AUTHORS: Albert Heinle,     aheinle@uwaterloo.ca
@*       Viktor Levandovskyy,     levandov@math.rwth-aachen.de

OVERVIEW: In this library, new methods for factorization on polynomials
  are implemented for two algebras, both generated by two generators (Weyl and
  shift algebras) over a field K. Recall, that   the first Weyl algebra over K
  is generated by x,d obeying the relation d*x=x*d+1.
@* The first shift algebra over K is generated by x,s obeying the relation s*x=x*s+s.
@* More detailled description of the algorithms can be found at
@url{http://www.math.rwth-aachen.de/\~Albert.Heinle}.

Guide: We are interested in computing a tree of factorizations, that is at the moment
a list of all found factorizations is returned. It may contain factorizations, which
are further reducible.

PROCEDURES:
  facFirstWeyl(h);           factorization in the first Weyl algebra
  testNCfac(l[,h[,1]]);      tests factorizations from a given list for correctness
  facSubWeyl(h,X,D);         factorization in the first Weyl algebra as a subalgebra
  facFirstShift(h);          factorization in the first shift algebra
  homogfacFirstQWeyl(h);     [-1,1]-homogeneous factorization in the first Q-Weyl algebra
  homogfacFirstQWeyl_all(h); [-1,1] homogeneous factorization(complete) in the first Q-Weyl algebra
  tst_ncfactor();            Runs the examples of all contained not static functions. Test thing.
";

LIB "general.lib";
LIB "nctools.lib";
LIB "involut.lib";
LIB "freegb.lib"; // for isVar
LIB "crypto.lib"; //for introot
LIB "matrix.lib"; //for submatrix
LIB "solve.lib"; //right now not needed
LIB "poly.lib"; //for content

proc tst_ncfactor()
"
A little test if the library works correct.
Runs simply all examples of non-static functions.
"
{
  example facFirstWeyl;
  example facFirstShift;
  example facSubWeyl;
  example testNCfac;
  example homogfacFirstQWeyl;
  example homogfacFirstQWeyl_all;
}
example
{
  "EXAMPLE:";echo=2;
   tst_ncfactor();
}

/////////////////////////////////////////////////////
//==================================================*
//deletes double-entries in a list of factorization
//without evaluating the product.
static proc delete_dublicates_noteval(list l)
"
INPUT: A list of lists; Output same as e.g. FacFirstWeyl. Containing different factorizations
       of a polynomial
OUTPUT: If there are dublicates in this list, this procedure deletes them and returns the list
 without double entries
"
{//proc delete_dublicates_noteval
  list result= l;
  int j; int k; int i;
  int deleted = 0;
  int is_equal;
  for (i = 1; i<= size(l); i++)
  {//Iterate over the different factorizations
    for (j = i+1; j<= size(l); j++)
    {//Compare the i'th factorization to the j'th
      if (size(l[i])!= size(l[j]))
      {//different sizes => not equal
        j++;
        continue;
      }//different sizes => not equal
      is_equal = 1;
      for (k = 1; k <= size(l[i]);k++)
      {//Compare every entry
        if (l[i][k]!=l[j][k])
        {
          is_equal = 0;
          break;
        }
      }//Compare every entry
      if (is_equal == 1)
      {//Delete this entry, because there is another equal one int the list
        result = delete(result, i-deleted);
        deleted = deleted+1;
        break;
      }//Delete this entry, because there is another equal one int the list
    }//Compare the i'th factorization to the j'th
  }//Iterate over the different factorizations
  return(result);
}//proc delete_dublicates_noteval

//==================================================
//deletes the double-entries in a list with
//evaluating the products
static proc delete_dublicates_eval(list l)
"
DEPRECATED
"
{//proc delete_dublicates_eval
  list result=l;
  int j; int k; int i;
  int deleted = 0;
  int is_equal;
  for (i = 1; i<= size(result); i++)
  {//Iterating over all elements in result
    for (j = i+1; j<= size(result); j++)
    {//comparing with the other elements
      if (product(result[i]) == product(result[j]))
      {//There are two equal results; throw away that one with the smaller size
        if (size(result[i])>=size(result[j]))
        {//result[i] has more entries
          result = delete(result,j);
          continue;
        }//result[i] has more entries
        else
        {//result[j] has more entries
          result = delete(result,i);
          i--;
          break;
        }//result[j] has more entries
      }//There are two equal results; throw away that one with the smaller size
    }//comparing with the other elements
  }//Iterating over all elements in result
  return(result);
}//proc delete_dublicates_eval


//==================================================*

static proc combinekfinlf(list g, int nof) //nof stands for "number of factors"
"
given a list of factors g and a desired size nof, this
procedure combines the factors, such that we recieve a
list of the length nof.
INPUT: A list of containing polynomials or any type where the *-operator is existent
OUTPUT: All possibilities (without permutation of the given list) to combine the polynomials
 into nof polynomials given by the user.
"
{//Procedure combinekfinlf
  list result;
  int i; int j; int k; //iteration variables
  list fc; //fc stands for "factors combined"
  list temp; //a temporary store for factors
  def nofgl = size(g); //nofgl stands for "number of factors of the given list"
  if (nofgl == 0)
  {//g was the empty list
    return(result);
  }//g was the empty list
  if (nof <= 0)
  {//The user wants to recieve a negative number or no element as a result
    return(result);
  }//The user wants to recieve a negative number or no element as a result
  if (nofgl == nof)
  {//There are no factors to combine
    result = result + list(g);
    return(result);
  }//There are no factors to combine
  if (nof == 1)
  {//User wants to get just one factor
    result = result + list(list(product(g)));
    return(result);
  }//User wants to get just one factor
  for (i = nof; i > 1; i--)
  {//computing the possibilities that have at least one original factor from g
    for (j = i; j>=1; j--)
    {//shifting the window of combinable factors to the left
      //fc below stands for "factors combined"
      fc = combinekfinlf(list(g[(j)..(j+nofgl - i)]),nof - i + 1);
      for (k = 1; k<=size(fc); k++)
      {//iterating over the different solutions of the smaller problem
        if (j>1)
        {//There are g_i before the combination
          if (j+nofgl -i < nofgl)
          {//There are g_i after the combination
            temp = list(g[1..(j-1)]) + fc[k] + list(g[(j+nofgl-i+1)..nofgl]);
          }//There are g_i after the combination
          else
          {//There are no g_i after the combination
            temp = list(g[1..(j-1)]) + fc[k];
          }//There are no g_i after the combination
        }//There are g_i before the combination
        if (j==1)
        {//There are no g_i before the combination
          if (j+ nofgl -i <nofgl)
          {//There are g_i after the combination
            temp = fc[k]+ list(g[(j + nofgl - i +1)..nofgl]);
          }//There are g_i after the combination
        }//There are no g_i before the combination
        result = result + list(temp);
      }//iterating over the different solutions of the smaller problem
    }//shifting the window of combinable factors to the left
  }//computing the possibilities that have at least one original factor from g
  for (i = 2; i<=nofgl div nof;i++)
  {//getting the other possible results
    result = result + combinekfinlf(list(product(list(g[1..i])))+list(g[(i+1)..nofgl]),nof);
  }//getting the other possible results
  result = delete_dublicates_noteval(result);
  return(result);
}//Procedure combinekfinlf


//==================================================*
//merges two sets of factors ignoring common
//factors
static proc merge_icf(list l1, list l2, intvec limits)
"
DEPRECATED
"
{//proc merge_icf
  list g;
  list f;
  int i; int j;
  if (size(l1)==0)
  {
    return(list());
  }
  if (size(l2)==0)
  {
    return(list());
  }
  if (size(l2)<=size(l1))
  {//l1 will be our g, l2 our f
    g = l1;
    f = l2;
  }//l1 will be our g, l2 our f
  else
  {//l1 will be our f, l2 our g
    g = l2;
    f = l1;
  }//l1 will be our f, l2 our g
  def result = combinekfinlf(g,size(f),limits);
  for (i = 1 ; i<= size(result); i++)
  {//Adding the factors of f to every possibility listed in temp
    for (j = 1; j<= size(f); j++)
    {
      result[i][j] = result[i][j]+f[j];
    }
    if(!limitcheck(result[i],limits))
    {
      result = delete(result,i);
      continue;
    }
    for (j = 1; j<=size(f);j++)
    {//Delete entry if there is a zero or an integer as a factor
      if (deg(result[i][j]) <= 0)
      {//found one
        result = delete(result,i);
        i--;
        break;
      }//found one
    }//Delete entry if there is a zero as factor
  }//Adding the factors of f to every possibility listed in temp
  return(result);
}//proc merge_icf

//==================================================*
//merges two sets of factors with respect to the occurrence
//of common factors
static proc merge_cf(list l1, list l2, intvec limits)
"
DEPRECATED
"
{//proc merge_cf
  list g;
  list f;
  int i; int j;
  list pre;
  list post;
  list candidate;
  list temp;
  int temppos;
  if (size(l1)==0)
  {//the first list is empty
    return(list());
  }//the first list is empty
  if(size(l2)==0)
  {//the second list is empty
    return(list());
  }//the second list is empty
  if (size(l2)<=size(l1))
  {//l1 will be our g, l2 our f
    g = l1;
    f = l2;
  }//l1 will be our g, l2 our f
  else
  {//l1 will be our f, l2 our g
    g = l2;
    f = l1;
  }//l1 will be our f, l2 our g
  list M;
  for (i = 2; i<size(f); i++)
  {//finding common factors of f and g...
    for (j=2; j<size(g);j++)
    {//... with g
      if (f[i] == g[j])
      {//we have an equal pair
        M = M + list(list(i,j));
      }//we have an equal pair
    }//... with g
  }//finding common factors of f and g...
  if (g[1]==f[1])
  {//Checking for the first elements to be equal
    M = M + list(list(1,1));
  }//Checking for the first elements to be equal
  if (g[size(g)]==f[size(f)])
  {//Checking for the last elements to be equal
    M = M + list(list(size(f),size(g)));
  }//Checking for the last elements to be equal
  list result;//= list(list());
  while(size(M)>0)
  {//set of equal pairs is not empty
    temp = M[1];
    temppos = 1;
    for (i = 2; i<=size(M); i++)
    {//finding the minimal element of M
      if (M[i][1]<=temp[1])
      {//a possible candidate that is smaller than temp could have been found
        if (M[i][1]==temp[1])
        {//In this case we must look at the second number
          if (M[i][2]< temp[2])
          {//the candidate is smaller
            temp = M[i];
            temppos = i;
          }//the candidate is smaller
        }//In this case we must look at the second number
        else
        {//The candidate is definately smaller
          temp = M[i];
          temppos = i;
        }//The candidate is definately smaller
      }//a possible candidate that is smaller than temp could have been found
    }//finding the minimal element of M
    M = delete(M, temppos);
    if(temp[1]>1)
    {//There are factors to combine before the equal factor
      if (temp[1]<size(f))
      {//The most common case
        //first the combinations ignoring common factors
        pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
        post = merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
        for (i = 1; i <= size(pre); i++)
        {//all possible pre's...
          for (j = 1; j<= size(post); j++)
          {//...combined with all possible post's
            candidate = pre[i]+list(f[temp[1]])+post[j];
            if (limitcheck(candidate,limits))
            {
              result = result + list(candidate);
            }
          }//...combined with all possible post's
        }//all possible pre's...
        //Now the combinations with respect to common factors
        post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
        if (size(post)>0)
        {//There are factors to combine
          for (i = 1; i <= size(pre); i++)
          {//all possible pre's...
            for (j = 1; j<= size(post); j++)
            {//...combined with all possible post's
              candidate= pre[i]+list(f[temp[1]])+post[j];
              if (limitcheck(candidate,limits))
              {
                result = result + list(candidate);
              }
            }//...combined with all possible post's
          }//all possible pre's...
        }//There are factors to combine
      }//The most common case
      else
      {//the last factor is the common one
        pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
        for (i = 1; i<= size(pre); i++)
        {//iterating over the possible pre-factors
          candidate = pre[i]+list(f[temp[1]]);
          if (limitcheck(candidate,limits))
          {
            result = result + list(candidate);
          }
        }//iterating over the possible pre-factors
      }//the last factor is the common one
    }//There are factors to combine before the equal factor
    else
    {//There are no factors to combine before the equal factor
      if (temp[1]<size(f))
      {//Just a check for security
        //first without common factors
        post=merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
        for (i = 1; i<=size(post); i++)
        {
          candidate = list(f[temp[1]])+post[i];
          if (limitcheck(candidate,limits))
          {
            result = result + list(candidate);
          }
        }
        //Now with common factors
        post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
        if(size(post)>0)
        {//we could find other combinations
          for (i = 1; i<=size(post); i++)
          {
            candidate = list(f[temp[1]])+post[i];
            if (limitcheck(candidate,limits))
            {
              result = result + list(candidate);
            }
          }
        }//we could find other combinations
      }//Just a check for security
    }//There are no factors to combine before the equal factor
  }//set of equal pairs is not empty
  for (i = 1; i <= size(result); i++)
  {//delete those combinations, who have an entry with degree less or equal 0
    for (j = 1; j<=size(result[i]);j++)
    {//Delete entry if there is a zero or an integer as a factor
      if (deg(result[i][j]) <= 0)
      {//found one
        result = delete(result,i);
        i--;
        break;
      }//found one
    }//Delete entry if there is a zero as factor
  }//delete those combinations, who have an entry with degree less or equal 0
  return(result);
}//proc merge_cf


//==================================================*
//merges two sets of factors

static proc mergence(list l1, list l2, intvec limits)
"
DEPRECATED
"
{//Procedure mergence
  list g;
  list f;
  int k; int i; int j;
  list F = list();
  list G = list();
  list tempEntry;
  list comb;
  if (size(l2)<=size(l1))
  {//l1 will be our g, l2 our f
    g = l1;
    f = l2;
  }//l1 will be our g, l2 our f
  else
  {//l1 will be our f, l2 our g
    g = l2;
    f = l1;
  }//l1 will be our f, l2 our g
  if (size(f)==1 or size(g)==1)
  {//One of them just has one entry
    if (size(f)== 1) {f = list(1) + f;}
    if (size(g) == 1) {g = list(1) + g;}
  }//One of them just has one entry
  //first, we need to add some latent -1's to the list f and to the list g in order
  //to get really all possibilities of combinations later
  for (i=1;i<=size(f)-1;i++)
  {//first iterator
    for (j=i+1;j<=size(f);j++)
    {//second iterator
      tempEntry = f;
      tempEntry[i] = (-1)*tempEntry[i];
      tempEntry[j] = (-1)*tempEntry[j];
      F = F + list(tempEntry);
    }//secont iterator
  }//first iterator
  F = F + list(f);
  //And now same game with g
  for (i=1;i<=size(g)-1;i++)
  {//first iterator
    for (j=i+1;j<=size(g);j++)
    {//second iterator
      tempEntry = g;
      tempEntry[i] = (-1)*tempEntry[i];
      tempEntry[j] = (-1)*tempEntry[j];
      G = G + list(tempEntry);
    }//secont iterator
  }//first iterator
  G = G + list(g);
  //Done with that

  list result;
  for (i = 1; i<=size(F); i++)
  {//Iterate over all entries in F
    for (j = 1;j<=size(G);j++)
    {//Same with G
      comb = combinekfinlf(F[i],2,limits);
      for (k = 1; k<= size(comb);k++)
      {//for all possibilities of combinations of the factors of f
        result = result + merge_cf(comb[k],G[j],limits);
        result = result + merge_icf(comb[k],G[j],limits);
        result = delete_dublicates_noteval(result);
      }//for all possibilities of combinations of the factors of f
    }//Same with G
  }//Iterate over all entries in F
  return(result);
}//Procedure mergence


//==================================================
//Checks, whether a list of factors doesn't exceed the given limits
static proc limitcheck(list g, intvec limits)
"
DEPRECATED
"
{//proc limitcheck
  int i;
  if (size(limits)!=3)
  {//check the input
    return(0);
  }//check the input
  if(size(g)==0)
  {
    return(0);
  }
  def prod = product(g);
  intvec iv11 = intvec(1,1);
  intvec iv10 = intvec(1,0);
  intvec iv01 = intvec(0,1);
  def limg = intvec(deg(prod,iv11) ,deg(prod,iv10),deg(prod,iv01));
  for (i = 1; i<=size(limg);i++)
  {//the final check
    if(limg[i]>limits[i])
    {
      return(0);
    }
  }//the final check
  return(1);
}//proc limitcheck


//==================================================*
//one factorization of a homogeneous polynomial
//in the first Weyl Algebra
static proc homogfacFirstWeyl(poly h)
"USAGE: homogfacFirstWeyl(h); h is a homogeneous polynomial in the
 first Weyl algebra with respect to the weight vector [-1,1]
RETURN: list
PURPOSE: Computes a factorization of a homogeneous polynomial h with
  respect to the weight vector [-1,1] in the first Weyl algebra
THEORY: @code{homogfacFirstWeyl} returns a list with a factorization of the given,
 [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with
 k positive, the last k entries in the output list are the second
 variable. If k is positive, the last k entries will be x. The other
 entries will be irreducible polynomials of degree zero or 1 resp. -1.
SEE ALSO: homogfacFirstWeyl_all
"{//proc homogfacFirstWeyl
  int p = printlevel-voice+2;//for dbprint
  def r = basering;
  poly hath;
  int i; int j;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  intvec ivm11 = intvec(-1,1);
  if (!homogwithorder(h,ivm11))
  {//The given polynomial is not homogeneous
    ERROR("Given polynomial was not [-1,1]-homogeneous");
    return(list());
  }//The given polynomial is not homogeneous
  if (h==0)
  {
    return(list(0));
  }
  list result;
  int m = deg(h,ivm11);
  dbprint(p,dbprintWhitespace +" Splitting the polynomial in A_0 and A_k-Part");
  if (m!=0)
  {//The degree is not zero
    if (m <0)
    {//There are more x than y
      hath = lift(var(1)^(-m),h)[1,1];
      for (i = 1; i<=-m; i++)
      {
        result = result + list(var(1));
      }
    }//There are more x than y
    else
    {//There are more y than x
      hath = lift(var(2)^m,h)[1,1];
      for (i = 1; i<=m;i++)
      {
        result = result + list(var(2));
      }
    }//There are more y than x
  }//The degree is not zero
  else
  {//The degree is zero
    hath = h;
  }//The degree is zero
  dbprint(p,dbprintWhitespace+" Done");
  //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1)
  list mons;
  dbprint(p,dbprintWhitespace+" Putting the monomials in the A_0-part in a list.");
  for(i = 1; i<=size(hath);i++)
  {//Putting the monomials in a list
    mons = mons+list(hath[i]);
  }//Putting the monomials in a list
  dbprint(p,dbprintWhitespace+" Done");
  dbprint(p,dbprintWhitespace+" Mapping this monomials to K[theta]");
  ring tempRing = 0,(x,y,theta),dp;
  setring tempRing;
  map thetamap = r,x,y;
  list mons = thetamap(mons);
  poly entry;
  for (i = 1; i<=size(mons);i++)
  {//transforming the monomials as monomials in theta
    entry = leadcoef(mons[i]);
    for (j = 0; j<leadexp(mons[i])[2];j++)
    {
      entry = entry * (theta-j);
    }
    mons[i] = entry;
  }//transforming the monomials as monomials in theta
  dbprint(p,dbprintWhitespace+" Done");
  dbprint(p,dbprintWhitespace+" Factorize the A_0-Part in K[theta]");
  list azeroresult = factorize(sum(mons));
  dbprint(p,dbprintWhitespace+" Successful");
  list azeroresult_return_form;
  for (i = 1; i<=size(azeroresult[1]);i++)
  {//rewrite the result of the commutative factorization
    for (j = 1; j <= azeroresult[2][i];j++)
    {
      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
    }
  }//rewrite the result of the commutative factorization
  dbprint(p,dbprintWhitespace+" Mapping back to A_0.");
  setring(r);
  map finalmap = tempRing,var(1),var(2),var(1)*var(2);
  list tempresult = finalmap(azeroresult_return_form);
  dbprint(p,dbprintWhitespace+"Successful.");
  for (i = 1; i<=size(tempresult);i++)
  {//factorizations of theta resp. theta +1
    if(tempresult[i]==var(1)*var(2))
    {
      tempresult = insert(tempresult,var(1),i-1);
      i++;
      tempresult[i]=var(2);
    }
    if(tempresult[i]==var(2)*var(1))
    {
      tempresult = insert(tempresult,var(2),i-1);
      i++;
      tempresult[i]=var(1);
    }
  }//factorizations of theta resp. theta +1
  result = tempresult+result;
  return(result);
}//proc homogfacFirstWeyl
/* example */
/* { */
/*      "EXAMPLE:";echo=2; */
/*      ring R = 0,(x,y),Ws(-1,1); */
/*      def r = nc_algebra(1,1); */
/*      setring(r); */
/*      poly h = (x^2*y^2+1)*(x^4); */
/*      homogfacFirstWeyl(h); */
/* } */

//==================================================
//Computes all possible homogeneous factorizations
static proc homogfacFirstWeyl_all(poly h)
"USAGE: homogfacFirstWeyl_all(h); h is a homogeneous polynomial in the first Weyl algebra
 with respect to the weight vector [-1,1]
RETURN: list
PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect
  to the weight vector [-1,1] in the first Weyl algebra
THEORY: @code{homogfacFirstWeyl} returns a list with all factorization of the given,
 homogeneous polynomial. It uses the output of homogfacFirstWeyl and permutes
 its entries with respect to the commutation rule. Furthermore, if a
 factor of degree zero is irreducible in K[  heta], but reducible in
 the first Weyl algebra, the permutations of this element with the other
 entries will also be computed.
SEE ALSO: homogfacFirstWeyl
"{//proc HomogfacFirstWeylAll
  int p=printlevel-voice+2;//for dbprint
  intvec iv11= intvec(1,1);
  if (deg(h,iv11) <= 0 )
  {//h is a constant
    dbprint(p,"Given polynomial was not homogeneous");
    return(list(list(h)));
  }//h is a constant
  def r = basering;
  list one_hom_fac; //stands for one homogeneous factorization
  int i; int j; int k;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  intvec ivm11 = intvec(-1,1);
  dbprint(p,dbprintWhitespace +" Calculate one homogeneous factorization using homogfacFirstWeyl");
  //Compute again a homogeneous factorization
  one_hom_fac = homogfacFirstWeyl(h);
  dbprint(p,dbprintWhitespace +"Successful");
  if (size(one_hom_fac) == 0)
  {//there is no homogeneous factorization or the polynomial was not homogeneous
    return(list());
  }//there is no homogeneous factorization or the polynomial was not homogeneous
  //divide list in A0-Part and a list of x's resp. y's
  list list_not_azero = list();
  list list_azero;
  list k_factor;
  int is_list_not_azero_empty = 1;
  int is_list_azero_empty = 1;
  k_factor = list(one_hom_fac[1]);
  if (absValue(deg(h,ivm11))<size(one_hom_fac)-1)
  {//There is a nontrivial A_0-part
    list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,ivm11)))];
    is_list_azero_empty = 0;
  }//There is a nontrivial A_0 part
  dbprint(p,dbprintWhitespace +" Combine x,y to xy in the factorization again.");
  for (i = 1; i<=size(list_azero)-1;i++)
  {//in homogfacFirstWeyl, we factorized theta, and this will be made undone
    if (list_azero[i] == var(1))
    {
      if (list_azero[i+1]==var(2))
      {
        list_azero[i] = var(1)*var(2);
        list_azero = delete(list_azero,i+1);
      }
    }
    if (list_azero[i] == var(2))
    {
      if (list_azero[i+1]==var(1))
      {
        list_azero[i] = var(2)*var(1);
        list_azero = delete(list_azero,i+1);
      }
    }
  }//in homogfacFirstWeyl, we factorized theta, and this will be made undone
  dbprint(p,dbprintWhitespace +" Done");
  if(deg(h,ivm11)!=0)
  {//list_not_azero is not empty
    list_not_azero =
      one_hom_fac[(size(one_hom_fac)-absValue(deg(h,ivm11))+1)..size(one_hom_fac)];
    is_list_not_azero_empty = 0;
  }//list_not_azero is not empty
  //Map list_azero in K[theta]
  dbprint(p,dbprintWhitespace +" Map list_azero to K[theta]");
  ring tempRing = 0,(x,y,theta), dp;
  setring(tempRing);
  poly entry;
  map thetamap = r,x,y;
  if(!is_list_not_azero_empty)
  {//Mapping in Singular is only possible, if the list before
    //contained at least one element of the other ring
    list list_not_azero = thetamap(list_not_azero);
  }//Mapping in Singular is only possible, if the list before
  //contained at least one element of the other ring
  if(!is_list_azero_empty)
  {//Mapping in Singular is only possible, if the list before
    //contained at least one element of the other ring
    list list_azero= thetamap(list_azero);
  }//Mapping in Singular is only possible, if the list before
  //contained at least one element of the other ring
  list k_factor = thetamap(k_factor);
  list tempmons;
  dbprint(p,dbprintWhitespace +" Done");
  for(i = 1; i<=size(list_azero);i++)
  {//rewrite the polynomials in A1 as polynomials in K[theta]
    tempmons = list();
    for (j = 1; j<=size(list_azero[i]);j++)
    {
      tempmons = tempmons + list(list_azero[i][j]);
    }
    for (j = 1 ; j<=size(tempmons);j++)
    {
      entry = leadcoef(tempmons[j]);
      for (k = 0; k < leadexp(tempmons[j])[2];k++)
      {
        entry = entry*(theta-k);
      }
      tempmons[j] = entry;
    }
    list_azero[i] = sum(tempmons);
  }//rewrite the polynomials in A1 as polynomials in K[theta]
  //Compute all permutations of the A0-part
  dbprint(p,dbprintWhitespace +" Compute all permutations of the A_0-part with the first resp.
the snd. variable");
  list result;
  int shift_sign;
  int shift;
  poly shiftvar;
  if (size(list_not_azero)!=0)
  {//Compute all possibilities to permute the x's resp. the y's in the list
    if (list_not_azero[1] == x)
    {//h had a negative weighted degree
      shift_sign = 1;
      shiftvar = x;
    }//h had a negative weighted degree
    else
    {//h had a positive weighted degree
      shift_sign = -1;
      shiftvar = y;
    }//h had a positive weighted degree
    result = permpp(list_azero + list_not_azero);
    for (i = 1; i<= size(result); i++)
    {//adjust the a_0-parts
      shift = 0;
      for (j=1; j<=size(result[i]);j++)
      {
        if (result[i][j]==shiftvar)
        {
          shift = shift + shift_sign;
        }
        else
        {
          result[i][j] = subst(result[i][j],theta,theta + shift);
        }
      }
    }//adjust the a_0-parts
  }//Compute all possibilities to permute the x's resp. the y's in the list
  else
  {//The result is just all the permutations of the a_0-part
    result = permpp(list_azero);
  }//The result is just all the permutations of the a_0 part
  if (size(result)==0)
  {
    return(result);
  }
  dbprint(p,dbprintWhitespace +" Done");
  dbprint(p,dbprintWhitespace +" Searching for theta resp. theta+1 in the list and fact. them");
  //Now we are going deeper and search for theta resp. theta + 1, substitute
  //them by xy resp. yx and go on permuting
  int found_theta;
  int thetapos;
  list leftpart;
  list rightpart;
  list lparts;
  list rparts;
  list tempadd;
  for (i = 1; i<=size(result) ; i++)
  {//checking every entry of result for theta or theta +1
    found_theta = 0;
    for(j=1;j<=size(result[i]);j++)
    {
      if (result[i][j]==theta)
      {//the jth entry is theta and can be written as x*y
        thetapos = j;
        result[i]= insert(result[i],x,j-1);
        j++;
        result[i][j] = y;
        found_theta = 1;
        break;
      }//the jth entry is theta and can be written as x*y
      if(result[i][j] == theta +1)
      {
        thetapos = j;
        result[i] = insert(result[i],y,j-1);
        j++;
        result[i][j] = x;
        found_theta = 1;
        break;
      }
    }
    if (found_theta)
    {//One entry was theta resp. theta +1
      leftpart = result[i];
      leftpart = leftpart[1..thetapos];
      rightpart = result[i];
      rightpart = rightpart[(thetapos+1)..size(rightpart)];
      lparts = list(leftpart);
      rparts = list(rightpart);
      //first deal with the left part
      if (leftpart[thetapos] == x)
      {
        shift_sign = 1;
        shiftvar = x;
      }
      else
      {
        shift_sign = -1;
        shiftvar = y;
      }
      for (j = size(leftpart); j>1;j--)
      {//drip x resp. y
        if (leftpart[j-1]==shiftvar)
        {//commutative
          j--;
          continue;
        }//commutative
        if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
        {//stop here
          break;
        }//stop here
        //Here, we can only have a a0- part
        leftpart[j] = subst(leftpart[j-1],theta, theta + shift_sign);
        leftpart[j-1] = shiftvar;
        lparts = lparts + list(leftpart);
      }//drip x resp. y
      //and now deal with the right part
      if (rightpart[1] == x)
      {
        shift_sign = 1;
        shiftvar = x;
      }
      else
      {
        shift_sign = -1;
        shiftvar = y;
      }
      for (j = 1 ; j < size(rightpart); j++)
      {
        if (rightpart[j+1] == shiftvar)
        {
          j++;
          continue;
        }
        if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
        {
          break;
        }
        rightpart[j] = subst(rightpart[j+1], theta, theta - shift_sign);
        rightpart[j+1] = shiftvar;
        rparts = rparts + list(rightpart);
      }
      //And now, we put all possibilities together
      tempadd = list();
      for (j = 1; j<=size(lparts); j++)
      {
        for (k = 1; k<=size(rparts);k++)
        {
          tempadd = tempadd + list(lparts[j]+rparts[k]);
        }
      }
      tempadd = delete(tempadd,1); // The first entry is already in the list
      result = result + tempadd;
      continue; //We can may be not be done already with the ith entry
    }//One entry was theta resp. theta +1
  }//checking every entry of result for theta or theta +1
  dbprint(p,dbprintWhitespace +" Done");
  //map back to the basering
  dbprint(p,dbprintWhitespace +" Mapping back everything to the basering");
  setring(r);
  map finalmap = tempRing, var(1), var(2),var(1)*var(2);
  list result = finalmap(result);
  for (i=1; i<=size(result);i++)
  {//adding the K factor
    result[i] = k_factor + result[i];
  }//adding the k-factor
  dbprint(p,dbprintWhitespace +" Done");
  dbprint(p,dbprintWhitespace +" Delete double entries in the list.");
  result = delete_dublicates_noteval(result);
  dbprint(p,dbprintWhitespace +" Done");
  return(result);
}//proc HomogfacFirstWeylAll
/* example */
/* { */
/*      "EXAMPLE:";echo=2; */
/*      ring R = 0,(x,y),Ws(-1,1); */
/*      def r = nc_algebra(1,1); */
/*      setring(r); */
/*      poly h = (x^2*y^2+1)*(x^4); */
/*      homogfacFirstWeyl_all(h); */
/* } */

//==================================================*
//Computes all permutations of a given list
static proc perm(list l)
"
DEPRECATED
"
{//proc perm
  int i; int j;
  list tempresult;
  list result;
  if (size(l)==0)
  {
    return(list());
  }
  if (size(l)==1)
  {
    return(list(l));
  }
  for (i = 1; i<=size(l); i++ )
  {
    tempresult = perm(delete(l,i));
    for (j = 1; j<=size(tempresult);j++)
    {
      tempresult[j] = list(l[i])+tempresult[j];
    }
    result = result+tempresult;
  }
  return(result);
}//proc perm

//==================================================
//computes all permutations of a given list by
//ignoring equal entries (faster than perm)
static proc permpp(list l)
"
INPUT: A list with entries of a type, where the ==-operator is defined
OUTPUT: A list with all permutations of this given list.
"
{//proc permpp
  int i; int j;
  list tempresult;
  list l_without_double;
  list l_without_double_pos;
  int double_entry;
  list result;
  if (size(l)==0)
  {
    return(list());
  }
  if (size(l)==1)
  {
    return(list(l));
  }
  for (i = 1; i<=size(l);i++)
  {//Filling the list with unique entries
    double_entry = 0;
    for (j = 1; j<=size(l_without_double);j++)
    {
      if (l_without_double[j] == l[i])
      {
        double_entry = 1;
        break;
      }
    }
    if (!double_entry)
    {
      l_without_double = l_without_double + list(l[i]);
      l_without_double_pos = l_without_double_pos + list(i);
    }
  }//Filling the list with unique entries
  for (i = 1; i<=size(l_without_double); i++ )
  {
    tempresult = permpp(delete(l,l_without_double_pos[i]));
    for (j = 1; j<=size(tempresult);j++)
    {
      tempresult[j] = list(l_without_double[i])+tempresult[j];
    }
    result = result+tempresult;
  }
  return(result);
}//proc permpp

//==================================================
//factorization of the first Weyl Algebra

//The following procedure just serves the purpose to
//transform the input into an appropriate input for
//the procedure sfacwa, where the ring must contain the
//variables in a certain order.
proc facFirstWeyl(poly h)
"USAGE: facFirstWeyl(h); h a polynomial in the first Weyl algebra
RETURN: list
PURPOSE: compute all factorizations of a polynomial in the first Weyl algebra
THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle
ASSUME: basering is the first Weyl algebra
NOTE: Every entry of the output list is a list with factors for one possible factorization.
The first factor is always a constant (1, if no nontrivial constant could be excluded).
EXAMPLE: example facFirstWeyl; shows examples
SEE ALSO: facSubWeyl, testNCfac, facFirstShift
"{//proc facFirstWeyl
  //Definition of printlevel variable
  int p = printlevel-voice+2;
  int i;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  dbprint(p,dbprintWhitespace +" Checking if the given algebra is a Weyl algebra");
  //Redefine the ring in my standard form
  if (!isWeyl())
  {//Our basering is not the Weyl algebra
    ERROR("Ring was not the first Weyl algebra");
    return(list());
  }//Our basering is not the Weyl algebra
  dbprint(p,dbprintWhitespace +" Successful");
  dbprint(p,dbprintWhitespace +" Checking, if the given ring is the first Weyl algebra");
  if(nvars(basering)!=2)
  {//Our basering is the Weyl algebra, but not the first
    ERROR("Ring is not the first Weyl algebra");
    return(list());
  }//Our basering is the Weyl algebra, but not the first

  //A last check before we start the real business: Is h already given as a polynomial just
  //in one variable?
  if (deg(h,intvec(1,0))== 0 or deg(h,intvec(0,1)) == 0)
  {//h is in K[x] or in K[d]
    if (deg(h,intvec(1,0))== 0 and deg(h,intvec(0,1)) == 0)
    {//We just have a constant
      return(list(list(h)));
    }//We just have a constant
    dbprint(p,dbprintWhitespace+"Polynomial was given in one variable.
Performing commutative factorization.");
    int theCommVar;
    if (deg(h,intvec(1,0)) == 0)
    {//The second variable is the variable to factorize
      theCommVar = 2;
    }//The second variable is the variable to factorize
    else{theCommVar = 1;}
    def r = basering;
    ring tempRing = 0,(var(theCommVar)),dp;
    if (theCommVar == 1){map mapToCommutative = r,var(1),1;}
    else {map mapToCommutative = r,1,var(1);}
    poly h = mapToCommutative(h);
    list tempResult = factorize(h);
    list result = list(list());
    int j;
    for (i = 1; i<=size(tempResult[1]); i++)
    {
      for (j = 1; j<=tempResult[2][i]; j++)
      {
        result[1] = result[1] + list(tempResult[1][i]);
      }
    }
    //mapping back
    setring(r);
    map mapBackFromCommutative = tempRing,var(theCommVar);
    def result = mapBackFromCommutative(result);
    dbprint(p,dbprintWhitespace+"result:");
    dbprint(p,result);
    dbprint(p,dbprintWhitespace+"Computing all permutations of this factorization");
    poly constantFactor = result[1][1];
    result[1] = delete(result[1],1);//Deleting the constant factor
    result=permpp(result[1]);
    for (i = 1; i<=size(result);i++)
    {//Insert constant factor
      result[i] = insert(result[i],constantFactor);
    }//Insert constant factor
    dbprint(p,dbprintWhitespace+"Done.");
    return(result);
  }//h is in K[x] or in K[d]


  dbprint(p,dbprintWhitespace +" Successful");
  list result = list();
  int j; int k; int l; //counter
  if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why
  {
    dbprint(p,dbprintWhitespace +" positions of the variables have to be switched");
    def r = basering;
    ring tempRing = ringlist(r)[1][1],(x,y),Ws(-1,1); // very strange:
    // setting Wp(-1,1) leads to SegFault; to clarify why!!!
    def NTR = nc_algebra(1,1);
    setring NTR ;
    map transf = r, var(2), var(1);
    dbprint(p,dbprintWhitespace +" Successful");
    list resulttemp = sfacwa(h);
    setring(r);
    map transfback = NTR, var(2),var(1);
    result = transfback(resulttemp);
  }
  else
  {
    dbprint(p, dbprintWhitespace +" factorization of the polynomial with the routine sfacwa");
    result = sfacwa(h);
    dbprint(p,dbprintWhitespace +" Done");
  }
  if (homogwithorder(h,intvec(-1,1)))
  {
    dbprint(p, dbprintWhitespace + " Polynomial was homogeneous, therefore we have
already a complete factorization and do not have to go through the factors recursively.");
    return(result);
  }
  dbprint(p,dbprintWhitespace + "We have the following intermediate list of inhomogeneous
factorizations:");
  dbprint(p,result);
  dbprint(p,dbprintWhitespace +" recursively check factors for irreducibility");
  list recursivetemp;
  int changedSomething;
  for(i = 1; i<=size(result);i++)
  {//recursively factorize factors
    if(size(result[i])>2)
    {//Nontrivial factorization
      for (j=2;j<=size(result[i]);j++)
      {//Factorize every factor
        recursivetemp = facFirstWeyl(result[i][j]);
        //if(size(recursivetemp)>1)
        //{//we have a nontrivial factorization
        changedSomething = 0;
        for(k=1; k<=size(recursivetemp);k++)
        {//insert factorized factors
          if(size(recursivetemp[k])>2)
          {//nontrivial
            changedSomething = 1;
            result = insert(result,result[i],i);
            for(l = size(recursivetemp[k]);l>=2;l--)
            {
              result[i+1] = insert(result[i+1],recursivetemp[k][l],j);
            }
            result[i+1] = delete(result[i+1],j);
          }//nontrivial
        }//insert factorized factors
        if (changedSomething)
        {
          result = delete(result,i);
        }
        //}//we have a nontrivial factorization
      }//Factorize every factor
    }//Nontrivial factorization
  }//recursively factorize factors
  dbprint(p,dbprintWhitespace +" Done");
  if (size(result)==0)
  {//only the trivial factorization could be found
    result = list(list(1,h));
  }//only the trivial factorization could be found
  list resultWithInterchanges;
  dbprint(p,dbprintWhitespace+ "And the result without interchanges with homogeneous factors is:");
  dbprint(p,result);
  for (i = 1; i <= size(result) ; i++)
  {//applying the interchanges to result
    resultWithInterchanges = resultWithInterchanges +
                             checkForHomogInhomogInterchangability(result[i],2,size(result[i]));
  }//applying the interchanges to result
  dbprint(p,dbprintWhitespace + "With interchanges, the result is:");
  dbprint(p,resultWithInterchanges);
  //now, refine the possible redundant list
  return( delete_dublicates_noteval(resultWithInterchanges) );
}//proc facFirstWeyl
example
{
  "EXAMPLE:";echo=2;
  ring R = 0,(x,y),dp;
  def r = nc_algebra(1,1);
  setring(r);
  poly h = (x^2*y^2+x)*(x+1);
  facFirstWeyl(h);
}

//////////////////////////////////////////////////
/////BRANDNEW!!!!////////////////////
//////////////////////////////////////////////////

static proc checkForHomogInhomogInterchangability(list factors,int posLeft,int posRight)
"
INPUT:  A list consisting of factors of a certain polynomial in the first Weyl
        algebra, factors, and a position from the left and the right, where the last swap was done.
OUTPUT: A list containing lists consisting of factors of a certain polynomial in the first Weyl
        algebra.
The purpose of this function is to check whether we can interchange certain inhomogeneous factors
with homogeneous ones. If it is possible, this function returns a list of lists
of possible interchanges.

The idea came because of an example, where we need an extra swap in the end, otherwise we would
not capture all factorizations. The example was
h = x4d7+11x3d6+x2d7+x2d6+x3d4+29x2d5+xd6+8xd5+d6+5x2d3+14xd4+13d4+5xd2+d3+d;

ASSUMPTIONS:

- All factors are irreducible
"
{//checkForHomogInhomogInterchangability
  int p = printlevel-voice+2;
  string dbprintWhitespace = "";
  int i; int j; int k;
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  if (size(factors) <= 2 || posLeft >= posRight - 1)
  {//easiest case: There is nothing to swap
    return (list(factors));
  }//easiest case: There is nothing to swap
  list result = list(factors);
  list tempResultEntries;
  list tempSwaps;
  list tempSwapsTempEntry;
  list attemptToSwap;
  intvec ivm11 = intvec(-1,1);
  dbprint(p, dbprintWhitespace+"We try to swap elements in the following list:");
  dbprint(p, factors);
  for (i = posLeft; i < posRight; i++)
  {//checking within the window posLeft <--> posRight, if there are interchanges possible
    if (homogwithorder(factors[i],ivm11) && !homogwithorder(factors[i+1],ivm11))
    {//position i is homogeneous, position i+1 is not ==> trying to swap
      attemptToSwap = extractHomogeneousDivisorsRight(factors[i]*factors[i+1]);
      if (size(attemptToSwap[1])>1)
      {//Bingo, we were able to swap this one element
        dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1));
        dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:");
        dbprint(p,attemptToSwap);
        tempSwapsTempEntry = list();
        for (j = size(factors); j >=1; j--)
        {//creating a new entry for the resulting list, replacing the swap in factors
          if (j==i+1)
          {
            for (k = size(attemptToSwap[1]); k >=1 ; k--)
            {
              tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][k]);
            }
            j--; //Because we changed entry i+1 and i
          }
          else
          {
            tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]);
          }
        }//creating a new entry for the resulting list, replacing the swap in factors
        tempSwaps = insert(tempSwaps,list(list(i+1,posRight),tempSwapsTempEntry));
      }//Bingo, we were able to swap this one element
    }//position i is homogeneous, position i+1 is not ==> trying to swap
    else
    {
      if(!homogwithorder(factors[i],ivm11) && homogwithorder(factors[i+1],ivm11))
      {//position i+1 is homogeneous, position i is not ==> trying to swap
        attemptToSwap = extractHomogeneousDivisorsLeft(factors[i]*factors[i+1]);
        if (size(attemptToSwap[1])>1)
        {//Bingo, we were able to swap this one element
          dbprint(p,dbprintWhitespace+"We can swap entry "+string(i)+" and "+ string(i+1));
          dbprint(p,dbprintWhitespace+"The elements look like the following after the swap:");
          dbprint(p,attemptToSwap);
          tempSwapsTempEntry = list();
          for (j = size(factors); j >=1; j--)
          {//creating a new entry for the resulting list, replacing the swap in factors
            if (j==i+1)
            {
              for (k = size(attemptToSwap[1]); k >=1 ; k--)
              {
                tempSwapsTempEntry = insert(tempSwapsTempEntry, attemptToSwap[1][k]);
              }
              j--; //Because we changed entry i+1 and i
            }
            else
            {
              tempSwapsTempEntry = insert(tempSwapsTempEntry,factors[j]);
            }
          }//creating a new entry for the resulting list, replacing the swap in factors
          tempSwaps = insert(tempSwaps,list(list(posLeft,i),tempSwapsTempEntry));
        }//Bingo, we were able to swap this one element
      }//position i+1 is homogeneous, position i is not ==> trying to swap
    }
  }//checking within the window posLeft <--> posRight, if there are interchanges possible
  //Now we will recursively call the function for all swapped entries.
  dbprint(p,dbprintWhitespace+ "Our list of different factorizations is now:");
  dbprint(p,tempSwaps);
  for (i = 1; i<=size(tempSwaps);i++)
  {//recursive call to all formerly attempted swaps.
    tempResultEntries=checkForHomogInhomogInterchangability(tempSwaps[i][2],
                                                            tempSwaps[i][1][1],tempSwaps[i][1][2]);
    result = result + tempResultEntries;
  }//recursive call to all formerly attempted swaps.
  result = delete_dublicates_noteval(result);
  return(result);
}//checkForHomogInhomogInterchangability

static proc sfacwa(poly h)
"INPUT: A polynomial h in the first Weyl algebra
OUTPUT: A list of factorizations, where the factors might still be reducible.
ASSUMPTIONS:
- Our basering is the first Weyl algebra; the x is the first variable,
  the differential operator the second.
"
{//proc sfacwa
  int i; int j; int k;
  int p = printlevel-voice+2;
  string dbprintWhitespace = "";
  number commonCoefficient = content(h);
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  dbprint(p,dbprintWhitespace + " Extracting homogeneous left and right factors");
  if(homogwithorder(h,intvec(-1,1)))
  {//we are already dealing with a -1,1 homogeneous poly
    dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous. Start homog.
fac. and ret. its result");
    return(homogfacFirstWeyl_all(h));
  }//we are already dealing with a -1,1 homogeneous poly
  list resulttemp = extractHomogeneousDivisors(h/commonCoefficient);
  //resulttemp = resulttemp + list(list(h/commonCoefficient));
  list inhomogeneousFactorsToFactorize;
  int isAlreadyInInhomogList;
  dbprint(p,dbprintWhitespace +" Done");
  dbprint(p,dbprintWhitespace +" Making Set of inhomogeneous polynomials we have to factorize.");
  for (i = 1; i<=size(resulttemp); i++)
  {//Going through all different kinds of factorizations where we extracted homogeneous factors
    for (j = 1;j<=size(resulttemp[i]);j++)
    {//searching for the inhomogeneous factor
      if (!homogwithorder(resulttemp[i][j],intvec(-1,1)))
      {//We have found our candidate
        isAlreadyInInhomogList = 0;
        for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++)
        {//Checking if our candidate is already in our tofactorize-list
          if (inhomogeneousFactorsToFactorize[k]==resulttemp[i][j])
          {//The candidate was already in the list
            isAlreadyInInhomogList = 1;
            break;
          }//The candidate was already in the list
        }//Checking if our candidate is already in our tofactorize-list
        if (!isAlreadyInInhomogList)
        {
          inhomogeneousFactorsToFactorize=inhomogeneousFactorsToFactorize + list(resulttemp[i][j]);
        }
      }//We have found our candidate
    }//searching for the inhomogeneous factor
  }//Going through all different kinds of factorizations where we extracted homogeneous factors
  dbprint(p,dbprintWhitespace +" Done");
  dbprint(p,dbprintWhitespace + "The set is:");
  dbprint(p,inhomogeneousFactorsToFactorize);
  dbprint(p,dbprintWhitespace+ "Factorizing the different occuring inhomogeneous factors");
  for (i = 1; i<= size(inhomogeneousFactorsToFactorize); i++)
  {//Factorizing all kinds of inhomogeneous factors
    inhomogeneousFactorsToFactorize[i] = sfacwa2(inhomogeneousFactorsToFactorize[i]);
    for (j = 1; j<=size(inhomogeneousFactorsToFactorize[i]);j++)
    {//Deleting the leading coefficient since we don't need him
      if (deg(inhomogeneousFactorsToFactorize[i][j][1],intvec(1,1))==0)
      {
        inhomogeneousFactorsToFactorize[i][j] = delete(inhomogeneousFactorsToFactorize[i][j],1);
      }
    }//Deleting the leading coefficient since we don't need him
  }//Factorizing all kinds of inhomogeneous factors
  dbprint(p,dbprintWhitespace +" Done");
  dbprint(p,dbprintWhitespace +" Putting the factorizations in the lists");
  list result;
  int posInhomogPoly;
  int posInhomogFac;
  for (i = 1; i<=size(resulttemp); i++)
  {//going through all by now calculated factorizations
    for (j = 1;j<=size(resulttemp[i]); j++)
    {//Finding the inhomogeneous factor
      if (!homogwithorder(resulttemp[i][j],intvec(-1,1)))
      {//Found it
        posInhomogPoly = j;
        break;
      }//Found it
    }//Finding the inhomogeneous factor
    for (k = 1; k<=size(inhomogeneousFactorsToFactorize);k++)
    {//Finding the matching inhomogeneous factorization we already determined
      if(product(inhomogeneousFactorsToFactorize[k][1]) == resulttemp[i][j])
      {//found it
        posInhomogFac = k;
        break;
      }//Found it
    }//Finding the matching inhomogeneous factorization we already determined
    for (j = 1; j <= size(inhomogeneousFactorsToFactorize[posInhomogFac]); j++)
    {
      result = insert(result, resulttemp[i]);
      result[1] = delete(result[1],posInhomogPoly);
      for (k =size(inhomogeneousFactorsToFactorize[posInhomogFac][j]);k>=1; k--)
      {//Inserting factorizations
        result[1] = insert(result[1],inhomogeneousFactorsToFactorize[posInhomogFac][j][k],
                           posInhomogPoly-1);
      }//Inserting factorizations
      dbprint(p,dbprintWhitespace + "Added a factorization to result, namely:");
      dbprint(p, result[1]);
    }
  }//going through all by now calculated factorizations
  dbprint(p,dbprintWhitespace +" Done");
  result = delete_dublicates_noteval(result);
  for (i = 1; i<=size(result);i++)
  {//Putting the content everywhere
    result[i] = insert(result[i],commonCoefficient);
  }//Putting the content everywhere
  return(result);
}//proc sfacwa

static proc sfacwa2(poly h)
"
Subprocedure of sfacwa
Assumptions:
- h is not in K[x] or in K[d], or even in K. These cases are caught by the input
- The coefficients are integer values and the gcd of the coefficients is 1
"
{//proc sfacwa2
  int p=printlevel-voice+2; // for dbprint
  int i;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  intvec ivm11 = intvec(-1,1);
  intvec iv11 = intvec(1,1);
  intvec iv10 = intvec(1,0);
  intvec iv01 = intvec(0,1);
  intvec iv1m1 = intvec(1,-1);
  poly p_max; poly p_min; poly q_max; poly q_min;
  map invo = basering,-var(1),var(2);
  list calculatedRightFactors;
  if(homogwithorder(h,ivm11))
  {//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here
    dbprint(p,dbprintWhitespace+" Given polynomial is -1,1 homogeneous.
Start homog. fac. and ret. its result");
    return(homogfacFirstWeyl_all(h));
  }//Unnecessary how we are using it, but if one wants to use it on its own, we are stating it here
  list result = list();
  int j; int k; int l;
  dbprint(p,dbprintWhitespace+" Computing the degree-limits of the factorization");
  //end finding the limits
  dbprint(p,dbprintWhitespace+" Computing the maximal and the minimal
homogeneous part of the given polynomial");
  list M = computeCombinationsMinMaxHomog(h);
  dbprint(p,dbprintWhitespace+" Done.");
  dbprint(p,dbprintWhitespace+" Filtering invalid combinations in M.");
  for (i = 1 ; i<= size(M); i++)
  {//filter valid combinations
    if (product(M[i]) == h)
    {//We have one factorization
      result = result + divides(M[i][1],h,invo,1);
      dbprint(p,dbprintWhitespace+"Result list updated:");
      dbprint(p,dbprintWhitespace+string(result));
      M = delete(M,i);
      continue;
    }//We have one factorization
  }//filter valid combinations
  dbprint(p,dbprintWhitespace+"Done.");
  dbprint(p,dbprintWhitespace+"The size of M is "+string(size(M)));
  for (i = 1; i<=size(M); i++)
  {//Iterate over all first combinations (p_max + p_min)(q_max + q_min)
    dbprint(p,dbprintWhitespace+" Combination No. "+string(i)+" in M:" );
    p_max = jet(M[i][1],deg(M[i][1],ivm11),ivm11)-jet(M[i][1],deg(M[i][1],ivm11)-1,ivm11);
    p_min = jet(M[i][1],deg(M[i][1],iv1m1),iv1m1)-jet(M[i][1],deg(M[i][1],iv1m1)-1,iv1m1);
    q_max = jet(M[i][2],deg(M[i][2],ivm11),ivm11)-jet(M[i][2],deg(M[i][2],ivm11)-1,ivm11);
    q_min = jet(M[i][2],deg(M[i][2],iv1m1),iv1m1)-jet(M[i][2],deg(M[i][2],iv1m1)-1,iv1m1);
    dbprint(p,dbprintWhitespace+" pmax = "+string(p_max));
    dbprint(p,dbprintWhitespace+" pmin = "+string(p_min));
    dbprint(p,dbprintWhitespace+" qmax = "+string(q_max));
    dbprint(p,dbprintWhitespace+" qmin = "+string(q_min));
    //Check, whether p_max + p_min or q_max and q_min are already left or right divisors.
    if (divides(p_min + p_max,h,invo))
    {
      dbprint(p,dbprintWhitespace+" Got one result.");
      result = result + divides(p_min + p_max,h,invo,1);
    }
    else
    {
      if (divides(q_min + q_max,h,invo))
      {
        dbprint(p,dbprintWhitespace+" Got one result.");
        result = result + divides(q_min + q_max, h , invo, 1);
      }
    }
    //Now the check, if deg(p_max) = deg(p_min)+1 (and the same with q_max and q_min)

    if (deg(p_max, ivm11) == deg(p_min, ivm11) +1 or  deg(q_max, ivm11) == deg(q_min, ivm11) +1 )
    {//Therefore, p_max + p_min must be a left factor or we can dismiss the combination
      dbprint(p,dbprintWhitespace+" There are no homogeneous parts we can put between
pmax and pmin resp. qmax and qmin.");
      //TODO: Prove, that then also a valid right factor is not possible
      M = delete(M,i);
      continue;
    }//Therefore, p_max + p_min must be a left factor or we can dismiss the combination

    //Done with the Check

    //If we come here, there are still homogeneous parts to be added to p_max + p_min
    //AND to q_max and q_min in
    //order to obtain a real factor
    //We use the procedure determineRestOfHomogParts to find our q.
    dbprint(p,dbprintWhitespace+" Solving for the other homogeneous parts in q");
    calculatedRightFactors = determineRestOfHomogParts(p_max,p_min,q_max,q_min,h);
    dbprint(p,dbprintWhitespace+" Done with it. Found "+string(size(calculatedRightFactors))
            +" solutions.");
    for (j = 1; j<=size(calculatedRightFactors);j++)
    {//Check out whether we really have right factors of h in calculatedRightFactors
      if (divides(calculatedRightFactors[j],h,invo))
      {
        result = result + divides(calculatedRightFactors[j],h,invo,1);
      }
      else
      {
        dbprint(p,"Solution for max and min homog found, but not a divisor of h");
        //TODO: Proof, why this can happen.
      }
    }//Check out whether we really have right factors of h in calculatedRightFactors
  }//Iterate over all first combinations (p_max + p_min)(q_max + q_min)


  result = delete_dublicates_noteval(result);
  //print(M);
  if (size(result) == 0)
  {//no factorization found
    result = list(list(h));
  }//no factorization found
  return(result);
}//proc sfacwa2

static proc determineRestOfHomogParts(poly pmax, poly pmin, poly qmax, poly qmin, poly h)
"INPUT: Polynomials p_max, p_min, q_max, q_min and h. The maximum homogeneous part h_max of h is
 given by p_max*pmin, the minimum homogeneous part h_min of h is given by p_min*q_min.
OUTPUT: A list of right factors q of h that have q_max and q_min as their maximum respectively
 minimum homogeneous part. Empty list, if those elements are not existent
ASSUMPTIONS:
 - deg(p_max,intvec(-1,1))>deg(p_min,intvec(-1,1)) +1
 - deg(q_max,intvec(-1,1))>deg(q_min,intvec(-1,1)) +1
 - p_max*q_max = h_max
 - p_min*q_min = h_min
 - The basering is the first Weyl algebra
"
{//proc determineRestOfHomogParts
  int p=printlevel-voice+2; // for dbprint
  string dbprintWhitespace = "";
  int i;
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  int kappa = Min(intvec(deg(h,intvec(1,0)), deg(h,intvec(0,1))));
  def R = basering;
  int n1 = deg(pmax,intvec(-1,1));
  int nk = -deg(pmin,intvec(1,-1));
  int m1 = deg(qmax,intvec(-1,1));
  int ml = -deg(qmin,intvec(1,-1));
  int j; int k;
  ideal mons;

  dbprint(p,dbprintWhitespace+" Extracting zero homog. parts of pmax, qmax, pmin, qmin and h.");
  //Extracting the zero homogeneous part of the given polynomials
  ideal pandqZero = pmax,pmin,qmax,qmin;
  if (n1 > 0){pandqZero[1] = lift(var(2)^n1,pmax)[1,1];}
  else{if (n1 < 0){pandqZero[1] = lift(var(1)^(-n1),pmax)[1,1];}
    else{pandqZero[1] = pmax;}}
  if (nk > 0){pandqZero[2] = lift(var(2)^nk,pmin)[1,1];}
  else{if (nk < 0){pandqZero[2] = lift(var(1)^(-nk),pmin)[1,1];}
    else{pandqZero[2] = pmin;}}
  if (m1 > 0){pandqZero[3] = lift(var(2)^m1,qmax)[1,1];}
  else{if (m1 < 0){pandqZero[3] = lift(var(1)^(-m1),qmax)[1,1];}
    else{pandqZero[3] = qmax;}}
  if (ml > 0){pandqZero[4] = lift(var(2)^ml,qmin)[1,1];}
  else{if (ml < 0){pandqZero[4] = lift(var(1)^(-ml),qmin)[1,1];}
    else{pandqZero[4] = qmin;}}
  list hZeroinR = homogDistribution(h);
  for (i = 1; i<=size(hZeroinR);i++)
  {//Extracting the zero homogeneous parts of the homogeneous summands of h
    if (hZeroinR[i][1] > 0){hZeroinR[i][2] = lift(var(2)^hZeroinR[i][1],hZeroinR[i][2])[1,1];}
    if (hZeroinR[i][1] < 0){hZeroinR[i][2] = lift(var(1)^(-hZeroinR[i][1]),hZeroinR[i][2])[1,1];}
  }//Extracting the zero homogeneous parts of the homogeneous summands of h
  dbprint(p,dbprintWhitespace+" Done!");
  //Moving everything into the ring K[theta]
  dbprint(p,dbprintWhitespace+" Moving everything into the ring K[theta]");
  ring KTheta = 0,(x,d,theta),dp;
  map thetamap = R, x, d;
  poly entry;
  ideal mons;
  ideal pandqZero;
  list hZeroinKTheta;
  setring(R);

  //Starting with p and q
  for (k=1; k<=4; k++)
  {//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys
    mons = ideal();
    for(i = 1; i<=size(pandqZero[k]);i++)
    {//Putting the monomials in a list
      mons[size(mons)+1] = pandqZero[k][i];
    }//Putting the monomials in a list
    setring(KTheta);
    mons = thetamap(mons);
    for (i = 1; i<=size(mons);i++)
    {//transforming the monomials as monomials in theta
      entry = leadcoef(mons[i]);
      for (j = 0; j<leadexp(mons[i])[2];j++)
      {
        entry = entry * (theta-j);
      }
      mons[i] = entry;
    }//transforming the monomials as monomials in theta
    pandqZero[size(pandqZero)+1] = sum(mons);
    setring(R);
  }//Transforming pmax(0),qmax(0),pmin(0),qmin(0) in theta-polys

  //Now hZero
  for (k = size(hZeroinR); k>= 1;k--)
  {//Transforming the different homogeneous parts of h into polys in K[theta]
    mons = ideal();
    for(i = 1; i<=size(hZeroinR[k][2]);i++)
    {//Putting the monomials in a list
      mons[size(mons)+1] = hZeroinR[k][2][i];
    }//Putting the monomials in a list
    setring(KTheta);
    mons = thetamap(mons);
    for (i = 1; i<=size(mons);i++)
    {//transforming the monomials as monomials in theta
      entry = leadcoef(mons[i]);
      for (j = 0; j<leadexp(mons[i])[2];j++)
      {
        entry = entry * (theta-j);
      }
      mons[i] = entry;
    }//transforming the monomials as monomials in theta
    hZeroinKTheta = hZeroinKTheta + list(sum(mons));
    setring(R);
  }//Transforming the different homogeneous parts of h into polys in K[theta]
  dbprint(p,dbprintWhitespace+" Done!");
  //Making the solutionRing
  ring solutionRing = 0,(theta,q(0..(kappa+1)*(m1-ml-1)-1)),lp;
  dbprint(p,dbprintWhitespace+" Our solution ring is given by "+ string(solutionRing));
  //mapping the different ps and qs and HZeros
  dbprint(p,dbprintWhitespace+" Setting up our solution system.");
  list ps;
  ideal pandqZero = imap(KTheta,pandqZero);
  ps[1] = list(n1,pandqZero[1]);
  ps[n1-nk+1] = list(nk,pandqZero[2]);
  for (i = 2; i<=n1-nk; i++)
  {
    ps[i] = list(n1-i+1,0);
  }
  list qs;
  qs[1] = list(m1,pandqZero[3]);
  qs[m1-ml+1] = list(ml,pandqZero[4]);
  for (i = 2; i<=m1-ml; i++)
  {
    qs[i] = list(m1-i+1,0);
    for (j = 0; j<=kappa; j++)
    {
      qs[i][2] = qs[i][2] + q((i-2)*(kappa+1)+j)*theta^j;
    }
  }
  list hZero = imap(KTheta,hZeroinKTheta);
  for (i = 1; i<=size(hZero); i++)
  {
    hZero[i] = list(n1+m1-i+1,hZero[i]);
  }

  //writing and solving the system
  list lhs;
  lhs[1] = list(ps[1][2],1);
  list rhs;
  rhs[n1-nk+1] = list(ps[size(ps)][2],1);
  for (i = 2; i<= n1-nk+1;i++)
  {
    lhs[i] = list(hZero[i][2],subst(qs[1][2],theta,theta+n1-i+1)*gammaForTheta(ps[i][1],qs[1][1]));
    rhs[n1-nk-i+2] = list(hZero[size(hZero)-i+1][2],subst(qs[m1-ml+1][2],theta,
                          theta+nk+i-1)*gammaForTheta(ps[n1-nk-i+2][1],qs[m1-ml+1][1]));
    for (j = 1; j<i; j++)
    {
      for (k = 1; k<=size(qs);k++)
      {
        if(ps[j][1]+qs[k][1] == hZero[i][1])
        {
          lhs[i][1] = lhs[i][1]*lhs[j][2];
          lhs[i][2] = lhs[i][2]*lhs[j][2];
          lhs[i][1] = lhs[i][1]-lhs[j][1]*subst(qs[k][2],theta, theta + n1- j +1)
            *gammaForTheta(ps[j][1],qs[k][1]);
        }
        if(ps[n1-nk+2-j][1] + qs[m1-ml+2-k][1] ==hZero[size(hZero)-i+1][1])
        {
          rhs[n1-nk-i+2][1] = rhs[n1-nk-i+2][1]*rhs[n1-nk+2-j][2];
          rhs[n1-nk-i+2][2] = rhs[n1-nk-i+2][2]*rhs[n1-nk+2-j][2];
          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,
                             theta + nk -j+1)*gammaForTheta(ps[n1-nk+2-j][1],qs[m1-ml+2-k][1]);
        }
      }
    }
  }
  list eqs;
  poly tempgcd;
  poly templhscoeff;
  poly temprhscoeff;
  for (i = 2; i<=n1-nk;i++)
  {
    if (gcd(rhs[i][2],lhs[i][2]) == 1)
    {
      eqs = eqs + list(lhs[i][1]*rhs[i][2] -rhs[i][1]*lhs[i][2]);
    }
    else
    {
      tempgcd = gcd(rhs[i][2],lhs[i][2]);
      templhscoeff = quotient(rhs[i][2],tempgcd)[1];
      temprhscoeff = quotient(lhs[i][2],tempgcd)[1];
      eqs = eqs + list(lhs[i][1]*templhscoeff -rhs[i][1]*temprhscoeff);
    }
  }
  matrix tempCoefMatrix;
  ideal solutionSystemforqs;
  for (i= 1; i<=size(eqs); i++)
  {
    tempCoefMatrix = coef(eqs[i],theta);
    solutionSystemforqs = solutionSystemforqs + ideal(submat(tempCoefMatrix,intvec(2),
                                                             intvec(1..ncols(tempCoefMatrix))));
  }
  dbprint(p,dbprintWhitespace+" Solution system for the coefficients of q is given by:");
  dbprint(p,solutionSystemforqs);
  option(redSB);
  dbprint(p,dbprintWhitespace+" Calculating reduced Groebner Basis of that system.");
  solutionSystemforqs = slimgb(solutionSystemforqs);
  dbprint(p,dbprintWhitespace+" Done!, the solution for the system is:");
  dbprint(p,dbprintWhitespace+string(solutionSystemforqs));
  if(vdim(slimgb(solutionSystemforqs+theta))==0)
  {//No solution in this case. Return the empty list
    dbprint(p,dbprintWhitespace+"The Groebner Basis of the solution system was <1>.");
    setring(R);
    return(list());
  }//No solution in this case. Return the empty list
  if(vdim(slimgb(solutionSystemforqs+theta))==-1)
  {//My conjecture is that this would never happen
    //ERROR("This is an counterexample to your conjecture. We have infinitely many solutions");
    //TODO: See, what we would do here
    dbprint(p,dbprintWhitespace+"There are infinitely many solution to this system.
We will return the empty list.");
    setring(R);
    return(list());
  }//My conjecture is that this would never happen
  else
  {//We have finitely many solutions
    if(vdim(slimgb(solutionSystemforqs+theta))==1)
    {//exactly one solution
      for (i = 2; i<= size(qs)-1;i++)
      {
        qs[i][2] = NF(qs[i][2],solutionSystemforqs);
      }
      setring(R);
      map backFromSolutionRing = solutionRing,var(1)*var(2);
      list qs = backFromSolutionRing(qs);
      list result = list(0);
      for (i = 1; i<=size(qs); i++)
      {
        if (qs[i][1]>0){qs[i][2] = qs[i][2]*var(2)^qs[i][1];}
        if (qs[i][1]<0){qs[i][2] = qs[i][2]*var(1)^(-qs[i][1]);}
        result[1] = result[1] + qs[i][2];
      }
      dbprint(p,dbprintWhitespace+"Found one unique solution. Returning the result.");
      return(result);
    }//exactly one solution
    else
    {//We have more than one solution, but finitely many
      dbprint(p,dbprintWhitespace+"Finitely many, but more than one solution.
Right now the conjecture is that this cannot happen.");
      dbprint(p,dbprintWhitespace+"Yay to counterexample :-) We can have more than one solution.");
      //TODO: Some theoretical work why this can never happan
      return(list());
    }//We have more than one solution, but finitely many
  }//We have finitely many solutions
}//proc determineRestOfHomogParts

static proc gammaForTheta(int j1,int j2)
"
INPUT: Two integers j1 and j2
OUTPUT: A polynomial in the first variable of the given ring. It calculates the following function:
   / 1,                                                if j1,j2>0 or j1,j2 <= 0
   | prod_{kappa = 0}^{|j1|-1}(var(1)-kappa),          if j1<0, j2>0, |j1|<=|j2|
gamma_{j1,j2}:= <  prod_{kappa = 0}^{j2-1}(var(1)-kappa-|j1|+|j2|),  if j1<0, j2>0, |j1|>|j2|
   | prod_{kappa = 1}^{j1}(var(1)+kappa),              if j1>0, j2<0, |j1|<=|j2|
   \ prod_{kappa = 1}^{|j2|}(\var(1)+kappa+|j1|-|j2|), if j1>0, j2<0, |j1|>|j2|
ASSUMPTION:
- Ring has at least one variable
"
{//gammaForTheta
  if (j1<=0 && j2 <=0) {return(1);}
  if (j1>=0 && j2 >=0){return(1);}
  poly result;
  int i;
  if (j1<0 && j2>0)
  {//case 2 or 3 from description above
    if (absValue(j1)<=absValue(j2))
    {//Case 2 holds here
      result = 1;
      for (i = 0;i<absValue(j1);i++)
      {
        result = result*(var(1)-i);
      }
      return(result);
    }//Case 2 holds here
    else
    {//Case 3 holds here
      result = 1;
      for (i = 0; i<j2; i++)
      {
        result = result*(var(1)-i-absValue(j1)+absValue(j2));
      }
      return(result);
    }//Case 3 holds here
  }//case 2 or 3 from description above
  else
  {//Case 4 or 5 from description above hold
    if (absValue(j1)<=absValue(j2))
    {//Case 4 holds
      result = 1;
      for (i = 1; i<=j1; i++)
      {
        result = result*(var(1)+i);
      }
      return(result);
    }//Case 4 holds
    else
    {//Case 5 holds
      result = 1;
      for (i = 1; i<=absValue(j2); i++)
      {
        result = result*(var(1)+i+absValue(j1)-absValue(j2));
      }
      return(result);
    }//Case 5 holds
  }//Case 4 or 5 from description above hold
}//gammaForTheta

static proc extractHomogeneousDivisors(poly h)
"INPUT: A polynomial h in the first Weyl algebra
OUTPUT: If h is homogeneous, then all factorizations of h are returned.
 If h is inhomogeneous, then a list l is returned whose entries
 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h and there
 exists an i in {1,...,n}, such that k_i is inhomogeneous and there
 is no homogeneous polynomial that divides this k_i neither from the
 left nor from the right. All the other entries in k are homogeneous
 polynomials.
"
{//extractHomogeneousDivisors
  int p=printlevel-voice+2; // for dbprint
  string dbprintWhitespace = "";
  int i; int j; int k; int l;
  list result;
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  if (homogwithorder(h,intvec(-1,1)))
  {//given polynomial was homogeneous already
    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
    result = homogfacFirstWeyl_all(h);
    for (i = 1; i<=size(result);i++)
    {//removing the first entry (coefficient) from the list result
      result[i] = delete(result[i],1);
    }//removing the first entry (coefficient) from the list result
    return(result);
  }//given polynomial was homogeneous already
  dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous left divisors extracted");
  list leftDivisionPossibilities = extractHomogeneousDivisorsLeft(h);
  dbprint(p,dbprintWhitespace+"Done. The result is:");
  dbprint(p,leftDivisionPossibilities);
  dbprint(p,dbprintWhitespace+"Calculating list with all homogeneous Right divisors extracted");
  list rightDivisionPossibilities = extractHomogeneousDivisorsRight(h);
  dbprint(p,dbprintWhitespace+"Done. The result is:");
  dbprint(p,rightDivisionPossibilities);
  list tempList;
  dbprint(p,dbprintWhitespace+"Calculating remaining right and left homogeneous divisors");
  for (i = 1; i<=size(leftDivisionPossibilities); i++)
  {//iterating through the list with extracted left divisors
    tempList = extractHomogeneousDivisorsRight
      (leftDivisionPossibilities[i][size(leftDivisionPossibilities[i])]);
    leftDivisionPossibilities[i] = delete(leftDivisionPossibilities[i],
                                          size(leftDivisionPossibilities[i]));
    for (j=1;j<=size(tempList);j++)
    {//Updating the list for Result
      tempList[j] = leftDivisionPossibilities[i] + tempList[j];
    }//Updating the list for Result
    result = result + tempList;
  }//iterating through the list with extracted left divisors
  for (i = 1; i<=size(rightDivisionPossibilities); i++)
  {//iterating through the list with extracted left divisors
    tempList = extractHomogeneousDivisorsLeft(rightDivisionPossibilities[i][1]);
    rightDivisionPossibilities[i] = delete(rightDivisionPossibilities[i],1);
    for (j=1;j<=size(tempList);j++)
    {//Updating the list for Result
      tempList[j] =  tempList[j]+rightDivisionPossibilities[i];
    }//Updating the list for Result
    result = result + tempList;
  }//iterating through the list with extracted left divisors
  dbprint(p,dbprintWhitespace+"Done");
  int posInhomog;
  poly hath = 1;
  list tempResult;
  dbprint(p,dbprintWhitespace+"Checking if we can swap left resp. right
divisors and updating result.");
  for (i = 1; i<= size(result); i++)
  {//Checking if we can swap left resp. right divisors
    for (j = 1; j<=size(result[i]);j++)
    {//finding the position of the inhomogeneous element in the list
      if(!homogwithorder(result[i][j],intvec(-1,1)))
      {
        posInhomog = j;
        break;
      }
    }//finding the position of the inhomogeneous element in the list
    hath = result[i][posInhomog];
    for(j=posInhomog-1;j>=1;j--)
    {
      hath = result[i][j]*hath;
      tempList = extractHomogeneousDivisorsRight(hath);
      if(size(tempList[1])==1)
      {//We could not swap this element to the right
        break;
      }//We could not swap this element to the right
      dbprint(p,dbprintWhitespace+"A swapping (left) of an element was possible");
      for(k = 1; k<=size(tempList);k++)
      {
        tempResult = insert(tempResult,result[i]);
        for (l = j;l<=posInhomog;l++)
        {
          tempResult[1] = delete(tempResult[1],j);
        }
        for (l = size(tempList[k]);l>=1;l--)
        {
          tempResult[1] = insert(tempResult[1],tempList[k][l],j-1);
        }
      }
    }
    hath = result[i][posInhomog];
    for(j=posInhomog+1;j<=size(result[i]);j++)
    {
      hath = hath*result[i][j];
      tempList = extractHomogeneousDivisorsLeft(hath);
      if(size(tempList[1])==1)
      {//We could not swap this element to the right
        break;
      }//We could not swap this element to the right
      dbprint(p,dbprintWhitespace+"A swapping (right) of an element was possible");
      for(k = 1; k<=size(tempList);k++)
      {
        tempResult = insert(tempResult,result[i]);
        for (l=posInhomog; l<=j;l++)
        {
          tempResult[1] = delete(tempResult[1],posInhomog);
        }
        for (l = size(tempList[k]);l>=1;l--)
        {
          tempResult[1] = insert(tempResult[1],tempList[k][l],posInhomog-1);
        }
      }
    }
  }//Checking if we can swap left resp. right divisors
  result = result + tempResult;
  result = delete_dublicates_noteval(result);
  return(result);
}//extractHomogeneousDivisors

static proc extractHomogeneousDivisorsLeft(poly h)
"INPUT: A polynomial h in the first Weyl algebra
OUTPUT: If h is homogeneous, then all factorizations of h are returned.
 If h is inhomogeneous, then a list l is returned whose entries
 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h.
 The entry k_n is inhomogeneous and has no other homogeneous
 left divisors any more.
 All the other entries in k are homogeneous
 polynomials.
"
{//extractHomogeneousDivisorsLeft
  int p=printlevel-voice+2; // for dbprint
  string dbprintWhitespace = "";
  int i;int j; int k;
  list result;
  poly hath;
  list recResult;
  map invo = basering,-var(1),var(2);
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  if (homogwithorder(h,intvec(-1,1)))
  {//given polynomial was homogeneous already
    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
    result = homogfacFirstWeyl_all(h);
    for (i = 1; i<=size(result);i++)
    {//removing the first entry (coefficient) from the list result
      result[i] = delete(result[i],1);
    }//removing the first entry (coefficient) from the list result
    return(result);
  }//given polynomial was homogeneous already
  list hlist = homogDistribution(h);
  dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands.");
  for (i = 1; i<= size(hlist); i++)
  {
    hlist[i] = homogfacFirstWeyl_all(hlist[i][2]);
    if (size(hlist[i][1])==1)
    {//One homogeneous part just has a trivial factorization
      if(hlist[i][1][1] == 0)
      {
        hlist = delete(hlist,i);
        continue;
      }
      else
      {
        return(list(list(h)));
      }
    }//One homogeneous part just has a trivial factorization
  }
  dbprint(p,dbprintWhitespace+ " Done.");
  dbprint(p,dbprintWhitespace+ " Trying to find Left divisors");
  list alreadyConsideredCandidates;
  poly candidate;
  int isCandidate;
  for (i = 1; i<=size(hlist[1]);i++)
  {//Finding candidates for homogeneous left divisors of h
    candidate = hlist[1][i][2];
    isCandidate = 0;
    for (j=1;j<=size(alreadyConsideredCandidates);j++)
    {
      if(alreadyConsideredCandidates[j] == candidate)
      {
        isCandidate =1;
        break;
      }
    }
    if(isCandidate)
    {
      i++;
      continue;
    }
    else
    {
      alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate);
    }
    dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous left divisor");
    for (j = 2; j<=size(hlist);j++)
    {//Iterating through the other homogeneous parts
      isCandidate = 0;
      for(k=1; k<=size(hlist[j]);k++)
      {
        if(hlist[j][k][2]==candidate)
        {
          isCandidate = 1;
          break;
        }
      }
      if(!isCandidate)
      {
        break;
      }
    }//Iterating through the other homogeneous parts
    if(isCandidate)
    {//candidate was really a left divisor
      dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous left divisor");
      hath = involution(lift(involution(candidate,invo),involution(h,invo))[1,1],invo);
      recResult = extractHomogeneousDivisorsLeft(hath);
      for (j = 1; j<=size(recResult); j++)
      {
        recResult[j] = insert(recResult[j],candidate);
      }
      result = result + recResult;
    }//Candidate was really a left divisor
  }//Finding candidates for homogeneous left divisors of h
  if (size(result)==0)
  {
    return(list(list(h)));
  }
  return(result);
}//extractHomogeneousDivisorsLeft

static proc extractHomogeneousDivisorsRight(poly h)
"INPUT: A polynomial h in the first Weyl algebra
OUTPUT: If h is homogeneous, then all factorizations of h are returned.
 If h is inhomogeneous, then a list l is returned whose entries
 are again lists k = [k_1,...,k_n], where k_1*...*k_n = h.
 The entry k_1 is inhomogeneous and has no other homogeneous
 right divisors any more.
 All the other entries in k are homogeneous
 polynomials.
"
{//extractHomogeneousDivisorsRight
  int p=printlevel-voice+2; // for dbprint
  string dbprintWhitespace = "";
  int i;int j; int k;
  list result;
  poly hath;
  list recResult;
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  if (homogwithorder(h,intvec(-1,1)))
  {//given polynomial was homogeneous already
    dbprint(p,dbprintWhitespace+"Polynomial was homogeneous. Just returning all factorizations.");
    result = homogfacFirstWeyl_all(h);
    for (i = 1; i<=size(result);i++)
    {//removing the first entry (coefficient) from the list result
      result[i] = delete(result[i],1);
    }//removing the first entry (coefficient) from the list result
    return(result);
  }//given polynomial was homogeneous already
  list hlist = homogDistribution(h);
  dbprint(p,dbprintWhitespace+ " Computing factorizations of all homogeneous summands.");
  for (i = 1; i<= size(hlist); i++)
  {
    hlist[i] = homogfacFirstWeyl_all(hlist[i][2]);
    if (size(hlist[i][1])==1)
    {//One homogeneous part just has a trivial factorization
      if(hlist[i][1][1] == 0)
      {
        hlist = delete(hlist,i);
        continue;
      }
      else
      {
        return(list(list(h)));
      }
    }//One homogeneous part just has a trivial factorization
  }
  dbprint(p,dbprintWhitespace+ " Done.");
  dbprint(p,dbprintWhitespace+ " Trying to find right divisors");
  list alreadyConsideredCandidates;
  poly candidate;
  int isCandidate;
  for (i = 1; i<=size(hlist[1]);i++)
  {//Finding candidates for homogeneous left divisors of h
    candidate = hlist[1][i][size(hlist[1][i])];
    isCandidate = 0;
    for (j=1;j<=size(alreadyConsideredCandidates);j++)
    {
      if(alreadyConsideredCandidates[j] == candidate)
      {
        isCandidate =1;
        break;
      }
    }
    if(isCandidate)
    {
      i++;
      continue;
    }
    else
    {
      alreadyConsideredCandidates = alreadyConsideredCandidates + list(candidate);
    }
    dbprint(p,dbprintWhitespace+"Checking if "+string(candidate)+" is a homogeneous r-divisor");
    for (j = 2; j<=size(hlist);j++)
    {//Iterating through the other homogeneous parts
      isCandidate = 0;
      for(k=1; k<=size(hlist[j]);k++)
      {
        if(hlist[j][k][size(hlist[j][k])]==candidate)
        {
          isCandidate = 1;
          break;
        }
      }
      if(!isCandidate)
      {
        break;
      }
    }//Iterating through the other homogeneous parts
    if(isCandidate)
    {//candidate was really a left divisor
      dbprint(p,dbprintWhitespace+string(candidate)+" is a homogeneous right divisor");
      hath = lift(candidate,h)[1,1];
      recResult = extractHomogeneousDivisorsRight(hath);
      for (j = 1; j<=size(recResult); j++)
      {
        recResult[j] = insert(recResult[j],candidate,size(recResult[j]));
      }
      result = result + recResult;
    }//Candidate was really a left divisor
  }//Finding candidates for homogeneous left divisors of h
  if (size(result)==0)
  {
    result = list(list(h));
  }
  return(result);
}//extractHomogeneousDivisorsRight

static proc fromZeroHomogToThetaPoly(poly h)
"
//DEPRECATED
INPUT: A polynomial h in the first Weyl algebra, homogeneous of degree 0
 OUTPUT: The ring Ktheta with a polynomial result representing h as polynomial in
  theta
 ASSUMPTIONS:
- h is homogeneous of degree 0 with respect to the [-1,1] weight vector
- The basering is the first Weyl algebra
"
{//proc fromZeroHomogToThetaPoly
  int i; int j;
  list mons;
  if(!homogwithorder(h,intvec(-1,1)))
  {//Input not a homogeneous polynomial causes an error
    ERROR("The input was not a homogeneous polynomial");
  }//Input not a homogeneous polynomial causes an error
  if(deg(h,intvec(-1,1))!=0)
  {//Input does not have degree 0
    ERROR("The input did not have degree 0");
  }//Input does not have degree 0
  for(i = 1; i<=size(h);i++)
  {//Putting the monomials in a list
    mons = mons+list(h[i]);
  }//Putting the monomials in a list
  ring KTheta = 0,(x,y,theta),dp;
  setring KTheta;
  map thetamap = r,x,y;
  list mons = thetamap(mons);
  poly entry;
  for (i = 1; i<=size(mons);i++)
  {//transforming the monomials as monomials in theta
    entry = leadcoef(mons[i]);
    for (j = 0; j<leadexp(mons[i])[2];j++)
    {
      entry = entry * (theta-j);
    }
    mons[i] = entry;
  }//transforming the monomials as monomials in theta
  poly result = sum(mons);
  keepring KTheta;
}//proc fromZeroHomogToThetaPoly

static proc fromHomogDistributionListToPoly(list l)
"returns the corresponding polynomial to a given output of homogDistribution"
{//proc fromHomogDistributionListToPoly
  poly result = 0;
  int i;
  for (i = 1; i<=size(l);i++)
  {
    result = result + l[i][2];
  }
  return(result);
}//proc fromHomogDistributionListToPoly

static proc divides(poly p1, poly p2,map invo, list #)
  "Tests, whether p1 divides p2 either from left or from right. The involution invo is needed
for checking both sides. The optional argument is needed in order to also return the other factor.
RETURN: If no optional argument is given, it will just return 1 or 0.
 Otherwise a list with at least one element
 Case 1: p1 does not divide p2 from any side. Then the output will be the empty list.
 Case 2: p2 does divide p2 from one side at least.
  Then it returns a list with tuples p,q, such that p or q equals p1 and
  pq = p2.
ASSUMPTIONS: - The map invo is an involution on the basering."
{//proc divides
  list result = list();
  poly tempfactor;
  if (involution(reduce(involution(p2,invo),std(involution(ideal(p1),invo))),invo)==0)
  {//p1 is a left divisor
    if(size(#)==0){return(1);}
    tempfactor = involution(lift(involution(p1,invo),involution(p2,invo))[1,1],invo);
    if (leadcoef(p1)<0 && leadcoef(tempfactor)<0)
    {//both have a negative leading coefficient
      result = result +list(list(1,-p1, -tempfactor));
    }//both have a negative leading coefficient
    else
    {
      if (leadcoef(p1)*leadcoef(tempfactor)<0)
      {//One of them has a negative leading coefficient
        if (leadcoef(p1)<0)
        {
          result = result + list(list(-1,-p1,tempfactor));
        }
        else
        {
          result = result + list(list(-1,p1,-tempfactor));
        }
      }//One of them has a negative leading coefficient
      else
      {//no negative coefficient at all
        result = result + list(list(1,p1,tempfactor));
      }//no negative coefficient at all
    }
  }//p1 is a left divisor

  if (reduce(p2,std(ideal(p1))) == 0)
  {//p1 is  a right divisor
    if(size(#)==0){return(1);}
    tempfactor = lift(p1, p2)[1,1];
    if (leadcoef(p1)<0 && leadcoef(tempfactor)<0)
    {//both have a negative leading coefficient
      result = result +list(list(1, -tempfactor,-p1));
    }//both have a negative leading coefficient
    else
    {
      if (leadcoef(p1)*leadcoef(tempfactor)<0)
      {//One of them has a negative leading coefficient
        if (leadcoef(p1)<0)
        {
          result = result + list(list(-1,tempfactor,-p1));
        }
        else
        {
          result = result + list(list(-1,-tempfactor,p1));
        }
      }//One of them has a negative leading coefficient
      else
      {//no negative coefficient at all
        result = result + list(list(1,tempfactor,p1));
      }//no negative coefficient at all
    }
  }//p1 is already a right divisor
  if (size(#)==0){return(0);}
  return(result);
}//proc divides

static proc computeCombinationsMinMaxHomog(poly h)
"Input: A polynomial h in the first Weyl Algebra
Output: Combinations of the form (p_max + p_min)(q_max + q_min), such that p_max, p_min,
 q_max and q_min are homogeneous and p_max*q_max equals the maximal homogeneous
 part in h, and p_max * q_max equals the minimal homogeneous part in h.
 h is not homogeneous.
"{//proc computeCombinationsMinMaxHomog
  intvec ivm11 = intvec(-1,1);
  intvec iv11 = intvec(1,1);
  intvec iv10 = intvec(1,0);
  intvec iv01 = intvec(0,1);
  intvec iv1m1 = intvec(1,-1);
  poly maxh = jet(h,deg(h,ivm11),ivm11)-jet(h,deg(h,ivm11)-1,ivm11);
  poly minh = jet(h,deg(h,iv1m1),iv1m1)-jet(h,deg(h,iv1m1)-1,iv1m1);
  list f1 = homogfacFirstWeyl_all(maxh);
  list f2 = homogfacFirstWeyl_all(minh);
  list result = list();
  int i;
  int j;
  //TODOs: Nur die kombinieren, die auch wirklich eine brauchbare Kombination sind.
  //Also, wir duerfen nicht aus den Geraden herausfliegen
  list pqmax = list();
  list pqmin = list();
  list tempList = list();

  //First, we are going to deal with our most hated guys: The Coefficients.
  //
  list coeffTuplesMax = getAllCoeffTuplesComb(factorizeInt(number(f1[1][1])));
  //We can assume without loss of generality, that p_max has a
  //nonnegative leading coefficient
  for (i = 1; i<=size(coeffTuplesMax);i++)
  {//Deleting all tuples with negative entries for p_max
    if (coeffTuplesMax[i][1]<0)
    {
      coeffTuplesMax = delete(coeffTuplesMax,i);
      continue;
    }
  }//Deleting all tuples with negative entries for p_max
  list coeffTuplesMin = getAllCoeffTuplesComb(factorizeInt(number(f2[1][1])));

  //Now, we will be actally dealing with the Combinations.
  //Let's start with the pqmax
  if (size(f1[1]) == 1)
  {//the maximal homogeneous factor is a constant
    pqmax = coeffTuplesMax;
  }//the maximal homogeneous factor is a constant
  else
  {//the maximal homogeneous factor is not a constant
    for (i = 1; i<=size(f1); i++)
    {//We can forget about the first coefficient now. Therefore we will delete him from the list.
      f1[i] = delete(f1[i],1);
      if(size(f1[i])==1)
      {//trivial thing
        for (j = 1; j<=size(coeffTuplesMax); j++)
        {
          pqmax = pqmax + list(list(coeffTuplesMax[j][1],coeffTuplesMax[j][2]*f1[i][1]));
          pqmax = pqmax + list(list(coeffTuplesMax[j][1]*f1[i][1],coeffTuplesMax[j][2]));
        }
        f1 = delete(f1,i);
        continue;
      }//trivial thing
    }//We can forget about the first coefficient now. Therefore we will delete him from the list.
    if (size(f1)>0)
    {
      tempList = getAllCombOfHomogFact(f1);
      for (i = 1; i<=size(tempList); i++)
      {//Every combination combined with the coefficient possibilities
        for (j = 1; j<=size(coeffTuplesMax); j++)
        {//iterating through the possible coefficient choices
          pqmax = pqmax + list(list(coeffTuplesMax[j][1]*tempList[i][1],
                                    coeffTuplesMax[j][2]*tempList[i][2]));
        }//iterating through the possible coefficient choices
      }//Every combination combined with the coefficient possibilities
      for (i = 1; i<=size(coeffTuplesMax); i++)
      {
        pqmax = pqmax + list(list(coeffTuplesMax[i][1],maxh/coeffTuplesMax[i][1]));
        pqmax = pqmax + list(list(maxh/coeffTuplesMax[i][2],coeffTuplesMax[i][2]));
      }
    }
  }//the maximal homogeneous factor is not a constant
  //Now we go to pqmin
  if (size(f2[1]) == 1)
  {//the minimal homogeneous factor is a constant
    pqmin = coeffTuplesMin;
  }//the minimal homogeneous factor is a constant
  else
  {//the minimal homogeneous factor is not a constant
    for (i = 1; i<=size(f2); i++)
    {//We can forget about the first coefficient now. Therefore we will delete him from the list.
      f2[i] = delete(f2[i],1);
      if(size(f2[i])==1)
      {//trivial thing
        for (j = 1; j<=size(coeffTuplesMin); j++)
        {
          pqmin = pqmin + list(list(coeffTuplesMin[j][1],coeffTuplesMin[j][2]*f2[i][1]));
          pqmin = pqmin + list(list(coeffTuplesMin[j][1]*f2[i][1],coeffTuplesMin[j][2]));
        }
        f2 = delete(f2,i);
        continue;
      }
    }//We can forget about the first coefficient now. Therefore we will delete him from the list.
    if(size(f2)>0)
    {
      tempList = getAllCombOfHomogFact(f2);
      for (i = 1; i<=size(tempList); i++)
      {//Every combination combined with the coefficient possibilities
        for (j = 1; j<=size(coeffTuplesMin); j++)
        {//iterating through the possible coefficient choices
          pqmin = pqmin + list(list(coeffTuplesMin[j][1]*tempList[i][1],
                                    coeffTuplesMin[j][2]*tempList[i][2]));
        }//iterating through the possible coefficient choices
      }//Every combination combined with the coefficient possibilities
      for (i = 1; i<=size(coeffTuplesMin); i++)
      {
        pqmin = pqmin + list(list(coeffTuplesMin[i][1],minh/coeffTuplesMin[i][1]));
        pqmin = pqmin + list(list(minh/coeffTuplesMin[i][2],coeffTuplesMin[i][2]));
      }
    }
  }//the minimal homogeneous factor is not a constant

  //and now we combine them together to obtain all possibilities.
  for (i = 1; i<=size(pqmax); i++)
  {//iterate over the maximal homogeneous combination possibilities
    for (j = 1; j<=size(pqmin); j++)
    {//iterate over the minimal homogeneous combiniation possibilities
      if (deg(pqmax[i][1], ivm11)>=deg(pqmin[j][1],ivm11) and deg(pqmax[i][2],
                                                                  ivm11)>=deg(pqmin[j][2],ivm11))
      {
        if (pqmax[i][1]+pqmin[j][1]!=0 and pqmax[i][2]+pqmin[j][2]!=0)
        {

          if (deg(h,ivm11)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),ivm11))
          {
            j++;
            continue;
          }
          if (deg(h,iv1m1)<=deg(h-(pqmax[i][1]+pqmin[j][1])*(pqmax[i][2]+pqmin[j][2]),iv1m1))
          {
            j++;
            continue;
          }
          result = result +list(list(pqmax[i][1]+pqmin[j][1],pqmax[i][2]+pqmin[j][2]));
        }
      }
    }//iterate over the minimal homogeneous combiniation possibilities
  }//iterate over the maximal homogeneous combination possibilities
  //Now deleting double entries
  result = delete_dublicates_noteval(result);
  return(result);
}//proc computeCombinationsMinMaxHomog

static proc getAllCombOfHomogFact(list l)
"Gets called in computeCombinationsMinMaxHomog. It gets a list of different homogeneous
factorizations of
one homogeneous polynomial and returns the possibilities to combine them into two factors.
Assumptions:
- The list does not contain the first coefficient.
- The list contains at least one list with two elements."
{//proc getAllCombOfHomogFact
  list result;
  list leftAndRightHandSides;
  int i; int j;
  list tempset;
  if (size(l)==1 and size(l[1])==2)
  {
    result = result + list(list(l[1][1],l[1][2]));
    return(result);
  }
  leftAndRightHandSides = getPossibilitiesForRightSides(l);
  for (i = 1; i<=size(leftAndRightHandSides); i++)
  {
    result =result+list(list(leftAndRightHandSides[i][1],product(leftAndRightHandSides[i][2][1])));
    //tidy up the right hand sides, because, if it is just one irreducible factor, we are done
    for (j = 1; j<=size(leftAndRightHandSides[i][2]);j++)
    {//Tidy up right hand sides
      if (size(leftAndRightHandSides[i][2][j])<2)
      {//Element can be dismissed
        leftAndRightHandSides[i][2] = delete(leftAndRightHandSides[i][2],j);
        continue;
      }//Element can be dismissed
    }//Tidy up right hand sides
    if (size(leftAndRightHandSides[i][2])>0)
    {
      tempset = getAllCombOfHomogFact(leftAndRightHandSides[i][2]);
      for (j = 1; j<=size(tempset);j++)
      {//multiplying the first factor with the left hand side
        result = result + list(list(leftAndRightHandSides[i][1]*tempset[j][1],tempset[j][2]));
      }//multiplying the first factor with the left hand side
    }
  }
  return(result);
}//proc getAllCombOfHomogFact

static proc getPossibilitiesForRightSides(list l)
"Given a list of different factorizations l, this function returns a list of the form
(a,{(a_2,...,a_n)| (a,a_2,...,a_n) in A})"
{//getPossibilitiesForRightSide
  list templ = l;
  list result;
  poly firstElement;
  list rightSides;
  list tempRightSide;
  int i; int j;
  while (size(templ)>0)
  {
    firstElement = templ[1][1];
    rightSides = list();
    for (i = 1; i<= size(templ); i++)
    {
      if (templ[i][1] == firstElement)
      {//save the right sides
        tempRightSide = list();
        for (j = 2; j<=size(templ[i]);j++)
        {
          tempRightSide = tempRightSide + list(templ[i][j]);
        }
        if (size(tempRightSide)!=0)
        {
          rightSides = rightSides + list(tempRightSide);
        }
        templ = delete(templ,i);
        continue;
      }//save the right sides
    }
    result = result + list(list(firstElement,rightSides));
  }
  return(result);
}//getPossibilitiesForRightSide

static proc getAllCoeffTuplesComb(list l)"
Given the output of factorizeInt ((a_1,...,a_n),(i_1,...,i_n)) , it returns all possible tuples
of the set {(a,b) | There exists an real N!=emptyset subset of {1,...,n}, such that
a = prod_{i \in N}a_i, b=prod_{i \not\in N} a_i}
Assumption: The list is sorted from smallest integer to highest.
- it is not the factorization of 0.
"
{//proc getAllCoeffTuplesComb
  list result;
  if (l[1][1] == 0)
  {
    ERROR("getAllCoeffTuplesComb: Zero Coefficients as leading and Tail Coeffs?
That is not possible. Something went wrong.");
  }
  if (size(l[1]) == 1)
  {//Trivial Factorization, just 1
    if (l[1][1] == 1)
    {
      return(list(list(1,1),list(-1,-1)));
    }
    else
    {
      return(list(list(-1,1),list(1,-1)));
    }
  }//Trivial Factorization, just 1
  if (size(l[1]) == 2 and l[2][2]==1)
  {//Just a prime number
    if (l[1][1] == 1)
    {
      result = list(list(l[1][2],1),list(1,l[1][2]));
      result = result + list(list(-l[1][2],-1),list(-1,-l[1][2]));
      return(result);
    }
    else
    {
      result = list(list(l[1][2],-1),list(1,-l[1][2]));
      result = result + list(list(-l[1][2],1),list(-1,l[1][2]));
      return(result);
    }
  }//Just a prime number
  //Now comes the interesting case: a product of primes
  list tempPrimeFactors;
  list tempPowersOfThem;
  int i;
  for (i = 2; i<=size(l[1]);i++)
  {//Removing the starting 1 or -1 to get the N's
    tempPrimeFactors[i-1] = l[1][i];
    tempPowersOfThem[i-1] = l[2][i];
  }//Removing the starting 1 or -1 to get the N's
  list Ns = getAllSubsetsN(list(tempPrimeFactors,tempPowersOfThem));
  list tempTuples;
  number productOfl = multiplyFactIntOutput(l);
  if (productOfl<0){productOfl = -productOfl;}
  tempTuples = tempTuples + list(list(1,productOfl),list(productOfl,1));
  for (i = 1; i<=size(Ns); i++)
  {
    if (productOfl/Ns[i]>Ns[i])
    {//TODO: BEWEISEN, dass das die einzigen Combos sind
      tempTuples = tempTuples + list(list(Ns[i],productOfl/Ns[i]),list(productOfl/Ns[i],Ns[i]));
    }//TODO: BEWEISEN, dass das die einzigen Combos sind
    if (productOfl/Ns[i]==Ns[i])
    {
      tempTuples = tempTuples + list(list(Ns[i],Ns[i]));
    }
  }
  //And now, it just remains to get the -1s and 1-s correctly to the tuples
  list tempEntry;
  if (l[1][1] == 1)
  {
    for (i = 1; i<=size(tempTuples);i++)
    {//Adding everything to result
      tempEntry = tempTuples[i];
      result = result + list(tempEntry);
      result = result + list(list(-tempEntry[1], -tempEntry[2]));
    }//Adding everyThing to Result
  }
  else
  {
    for (i = 1; i<=size(tempTuples);i++)
    {//Adding everything to result
      tempEntry = tempTuples[i];
      result = result + list(list(tempEntry[1],-tempEntry[2]));
      result = result + list(list(-tempEntry[1], tempEntry[2]));
    }//Adding everyThing to Result
  }
  return(result);
}//proc getAllCoeffTuplesComb

static proc contains(list l, int elem)
"Assumption: l is sorted"
{//Binary Search in list
  if (size(l)<=1)
  {
    if(size(l) == 0){return(0);}
    if (l[1]!=elem){return(0);}
    else{return(1);}
  }
  int imax = size(l);
  int imin = 1;
  int imid;
  while(imax >= imin)
  {
    imid = (imin + imax)/2;
    if (l[imid] == elem){return(1);}
    if (l[imid] <elem) {imin = imid +1;}
    else{imax = imid -1;}
  }
  return(0)
    }//Binary Search in list

static proc getAllSubsetsN(list l)
"
Assumptions:
- The list is containing two lists. They can be assumed to be outputs of the function
factorizeInt. They have at least one entry. If it is exactly one entry, the second intvec should
contain a value at least 2.
  "
{
  list primeFactors=l[1];
  list powersOfThem = l[2];
  int i;int j;
  //Casting the entries to be numbers
  for (i=1; i<=size(primeFactors); i++)
  {
    primeFactors[i] = number(primeFactors[i]);
    powersOfThem[i] = number(powersOfThem[i]);
  }

  //Done
  list result;
  list tempPrimeFactors;
  list tempPowersOfThem;
  list tempset;
  if (sum(powersOfThem) <=2)
  {//Easy Case
    return(list(primeFactors[1]));
  }//Easy Case
  if (size(primeFactors)==1)
  {//Also Easy Case
    for (j = 1; j<powersOfThem[1]; j++)
    {
      result = result + list(primeFactors[1]^j);
    }
    return(result);
  }//Also Easy Case
  for (i = 1; i<= size(primeFactors); i++)
  {//Going through every entry
    result = result + list(primeFactors[i]);
    if (i == size(primeFactors))
    {
      for (j = 1;j<powersOfThem[i];j++)
      {
        result = result + list (primeFactors[i]^j);
      }
      break;
    }
    if (powersOfThem[i]==1)
    {
      for (j = i+1;j<=size(primeFactors);j++)
      {
        tempPrimeFactors[j-i] = primeFactors[j];
        tempPowersOfThem[j-i] = powersOfThem[j];
      }
    }
    else
    {
      for (j = i; j<=size(primeFactors);j++)
      {
        tempPrimeFactors[j-i+1] = primeFactors[j];
        tempPowersOfThem[j-i+1] = powersOfThem[j];
        tempPowersOfThem[1] = tempPowersOfThem[1]-1;
      }
    }
    tempset = getAllSubsetsN(list(tempPrimeFactors,tempPowersOfThem));
    for (j = 1; j<=size(tempset); j++)
    {
      result = result +list((tempset[j])*(primeFactors[i]));
    }
  }//Going through every entry
  result = sort(result)[1];
  result = delete_dublicates_noteval(result);
  return(result);
}

static proc multiplyFactIntOutput(list l)
"Given the output of factorizeInt, this method computes the product of it."
{//proc multiplyFactIntOutput
  int i;
  number result = 1;
  for (i = 1; i<=size(l[1]); i++)
  {
    result = result*(l[1][i])^(l[2][i]);
  }
  return(result);
}//proc multiplyFactIntOutput

static proc fromListToIntvec(list l)
"Converter from List to intvec"
{
  intvec result; int i;
  for (i = 1; i<=size(l); i++)
  {
    result[i] = l[i];
  }
  return(result);
}

static proc fromIntvecToList(intvec l)"
Converter from intvec to list"
{//proc fromIntvecToList
  list result = list();
  int i;
  for (i = size(l); i>=1; i--)
  {
    result = insert(result, l[i]);
  }
  return(result);
}//proc fromIntvecToList


static proc factorizeInt(number n)
"Given an integer n, factorizeInt computes its factorization. The output is a list
containing two intvecs. The first contains the prime factors, the second its powers.
ASSUMPTIONS:
- n is given as integer number
"{
  if (n==0)
  {return(list(list(0),list(1)));}
  int i;
  list temp = primefactors(n);
  if (n<0)
  {list result = list(list(-1),list(1));}
  else
  {list result = list(list(1),list(1));}
  result[1] = result[1] + temp[1];
  result[2] = result[2] + temp[2];
  return(result);
}


static proc homogDistribution(poly h)
"Input: A polynomial in the first Weyl Algebra.
  Output: A two-dimensional list of the following form. Every sublist contains exactly two entries.
   One for the Z-degree of the corresponding homogeneous part (integer), and the homogeneous
   polynomial itself, and those sublists are oredered by ascending degree.
   For example a call of homogDistribution(x+d+1) would have the output
   [1]:
     [1]:
       -1
     [2]:
       x
   [2]:
     [1]:
       0
     [2]:
       1
   [3]:
     [1]:
       1
     [2]:
       d
"{//homogDistribution
  if (h == 0)
  {//trivial case where input is 0
    return(list(list(0,0)));
  }//trivial case where input is 0
  if (!isWeyl())
  {//Our basering is not the Weyl algebra
    ERROR("Ring was not the first Weyl algebra");
    return(list());
  }//Our basering is not the Weyl algebra
  if(nvars(basering)!=2)
  {//Our basering is the Weyl algebra, but not the first
    ERROR("Ring is not the first Weyl algebra");
    return(list());
  }//Our basering is the Weyl algebra, but not the first
  intvec ivm11 = intvec(-1,1);
  intvec iv1m1 = intvec(1,-1);
  poly tempH = h;
  poly minh;
  list result = list();
  int nextExpectedDegree = -deg(tempH,iv1m1);
  while (tempH != 0)
  {
    minh = jet(tempH,deg(tempH,iv1m1),iv1m1)-jet(tempH,deg(tempH,iv1m1)-1,iv1m1);
    while (deg(minh,ivm11)>nextExpectedDegree)
    {//filling empty homogeneous spaces with 0
      result = result + list(list(nextExpectedDegree,0));
      nextExpectedDegree = nextExpectedDegree +1;
    }//filling empty homogeneous spaces with 0
    result = result + list(list(deg(minh,ivm11),minh));
    tempH = tempH - minh;
    nextExpectedDegree = nextExpectedDegree +1;
  }
  return(result);
}//homogDistribution

static proc countHomogParts(poly h)
"Counts the homogeneous parts of a given polynomial h"
{
  int i;
  list outPutHD = homogDistribution(h);
  int result = 0;
  for (i = 1; i <=size(outPutHD); i++)
  {
    if (outPutHD[i][2] != 0){result++;}
  }
  return(result);
}

//////////////////////////////////////////////////
/////BRANDNEW!!!!////////////////////
//////////////////////////////////////////////////

//==================================================
/*Singular has no way implemented to test polynomials
  for homogenity with respect to a weight vector.
  The following procedure does exactly this*/
static proc homogwithorder(poly h, intvec weights)
{//proc homogwithorder
  if(size(weights) != nvars(basering))
  {//The user does not know how many variables the current ring has
    return(0);
  }//The user does not know how many variables the current ring has
  int i;
  int dofp = deg(h,weights); //degree of polynomial
  for (i = 1; i<=size(h);i++)
  {
    if (deg(h[i],weights)!=dofp)
    {
      return(0);
    }
  }
  return(1);
}//proc homogwithorder

//==================================================
//Testfac: Given a list with different factorizations of
// one polynomial, the following procedure checks
// whether they all refer to the same polynomial.
// If they do, the output will be a list, that contains
// the product of each factorization. If not, the empty
// list will be returned.
// If the optional argument # is given (i.e. the polynomial
// which is factorized by the elements of the given list),
// then we look, if the entries are factorizations of p
// and if not, a list with the products subtracted by p
// will be returned
proc testNCfac(list l, list #)
"USAGE: testNCfac(l[,p,b]); l is a list, p is an optional poly, b is 1 or 0
RETURN: Case 1: No optional argument. In this case the output is 1, if the
  entries in the given list represent the same polynomial or 0
  otherwise.
 Case 2: One optional argument p is given. In this case it returns 1,
  if all the entries in l are factorizations of p, otherwise 0.
 Case 3: Second optional b is given. In this case a list is returned
  containing the difference between the product of each entry in
  l and p.
ASSUME: basering is the first Weyl algebra, the entries of l are polynomials
PURPOSE: Checks whether a list of factorizations contains factorizations of
  the same element in the first Weyl algebra
THEORY: @code{testNCfac} multiplies out each factorization and checks whether
 each factorization was a factorization of the same element.
@* - if there is only a list given, the output will be 0, if it
     does not contain factorizations of the same element. Otherwise the output
     will be 1.
@* - if there is a polynomial in the second argument, then the procedure checks
     whether the given list contains factorizations of this polynomial. If it
     does, then the output depends on the third argument. If it is not given,
     the procedure will check whether the factorizations in the list
     l are associated to this polynomial and return either 1 or 0, respectively.
     If the third argument is given, the output will be a list with
     the length of the given one and in each entry is the product of one
     entry in l subtracted by the polynomial.
EXAMPLE: example testNCfac; shows examples
SEE ALSO: facFirstWeyl, facSubWeyl, facFirstShift
"{//proc testfac
  int p = printlevel - voice + 2;
  int i;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  dbprint(p,dbprintWhitespace + " Checking the input");
  if (size(l)==0)
  {//The empty list is given
    dbprint(p,dbprintWhitespace + " Given list was empty");
    return(list());
  }//The empty list is given
  if (size(#)>2)
  {//We want max. two optional arguments
    dbprint(p,dbprintWhitespace + " More than two optional arguments");
    return(list());
  }//We want max. two optional arguments
  dbprint(p,dbprintWhitespace + " Done");
  list result;
  int j;
  if (size(#)==0)
  {//No optional argument is given
    dbprint(p,dbprintWhitespace + " No optional arguments");
    int valid = 1;
    for (i = size(l);i>=1;i--)
    {//iterate over the elements of the given list
      if (size(result)>0)
      {
        if (product(l[i])!=result[size(l)-i])
        {
          valid = 0;
          break;
        }
      }
      result = insert(result, product(l[i]));
    }//iterate over the elements of the given list
    return(valid);
  }//No optional argument is given
  else
  {
    dbprint(p,dbprintWhitespace + " Optional arguments are given.");
    int valid = 1;
    for (i = size(l);i>=1;i--)
    {//iterate over the elements of the given list
      if (product(l[i])!=#[1])
      {
        valid = 0;
      }
      result = insert(result, product(l[i])-#[1]);
    }//iterate over the elements of the given list
    if(size(#)==2)
    {
      dbprint(p,dbprintWhitespace + " A third argument is given. Output is a list now.");
      return(result);
    }
    return(valid);
  }
}//proc testfac
example
{
  "EXAMPLE:";echo=2;
  ring r = 0,(x,y),dp;
  def R = nc_algebra(1,1);
  setring R;
  poly h = (x^2*y^2+1)*(x^2);
  def t1 = facFirstWeyl(h);
  //fist a correct list
  testNCfac(t1);
  //now a correct list with the factorized polynomial
  testNCfac(t1,h);
  //now we put in an incorrect list without a polynomial
  t1[3][3] = y;
  testNCfac(t1);
  // take h as additional input
  testNCfac(t1,h);
  // take h as additional input and output list of differences
  testNCfac(t1,h,1);
}
//==================================================
//Procedure facSubWeyl:
//This procedure serves the purpose to compute a
//factorization of a given polynomial in a ring, whose subring
//is the first Weyl algebra. The polynomial must only contain
//the two arguments, which are also given by the user.

proc facSubWeyl(poly h, def X, def D)
"USAGE:  facSubWeyl(h,x,D); h, X, D polynomials
RETURN: list
ASSUME: X and D are variables of a basering, which satisfy DX = XD +1.
@* That is,  they generate the copy of the first Weyl algebra in a basering.
@* Moreover, h is a polynomial in X and D only.
PURPOSE: compute factorizations of the polynomial, which depends on X and D.
EXAMPLE: example facSubWeyl; shows examples
SEE ALSO: facFirstWeyl, testNCfac, facFirstShift
"{
  int p = printlevel - voice + 2;
  dbprint(p," Start initial Checks of the input.");
  // basering can be anything having a Weyl algebra as subalgebra
  def @r = basering;
  //We begin to check the input for assumptions
  // which are: X,D are vars of the basering,
  if ( (isVar(X)!=1) || (isVar(D)!=1) || (size(X)>1) || (size(D)>1) ||
       (leadcoef(X) != number(1)) || (leadcoef(D) != number(1)) )
  {
    ERROR("expected pure variables as generators of a subalgebra");
  }
  // Weyl algebra:
  poly w = D*X-X*D-1; // [D,X]=1
  poly u = D*X-X*D+1; // [X,D]=1
  if (u*w!=0)
  {
    // that is no combination gives Weyl
    ERROR("2nd and 3rd argument do not generate a Weyl algebra");
  }
  // one of two is correct
  int isReverted = 0; // Reverted Weyl if dx=xd-1 holds
  if (u==0)
  {
    isReverted = 1;
  }
  // else: do nothing
  // DONE with assumptions, Input successfully checked
  dbprint(p," Successful");
  intvec lexpofX = leadexp(X);
  intvec lexpofD = leadexp(D);
  int varnumX=1;
  int varnumD=1;
  while(lexpofX[varnumX] != 1)
  {
    varnumX++;
  }
  while(lexpofD[varnumD] != 1)
  {
    varnumD++;
  }
  /* VL : to add printlevel stuff */
  dbprint(p," Change positions of the two variables in the list, if needed");
  if (isReverted)
  {
    ring firstweyl = 0,(var(varnumD),var(varnumX)),dp;
    def Firstweyl = nc_algebra(1,1);
    setring Firstweyl;
    ideal M = 0:nvars(@r);
    M[varnumX]=var(2);
    M[varnumD]=var(1);
    map Q = @r,M;
    poly h= Q(h);
  }
  else
  { // that is unReverted
    ring firstweyl = 0,(var(varnumX),var(varnumD)),dp;
    def Firstweyl = nc_algebra(1,1);
    setring Firstweyl;
    poly h= imap(@r,h);
  }
  dbprint(p," Done!");
  list result = facFirstWeyl(h);
  setring @r;
  list result;
  if (isReverted)
  {
    // map swap back
    ideal M; M[1] = var(varnumD); M[2] = var(varnumX);
    map S = Firstweyl, M;
    result = S(result);
  }
  else
  {
    // that is unReverted
    result = imap(Firstweyl,result);
  }
  return(result);
}//proc facSubWeyl
example
{
  "EXAMPLE:";echo=2;
  ring r = 0,(x,y,z),dp;
  matrix D[3][3]; D[1,3]=-1;
  def R = nc_algebra(1,D); // x,z generate Weyl subalgebra
  setring R;
  poly h = (x^2*z^2+x)*x;
  list fact1 = facSubWeyl(h,x,z);
  // compare with facFirstWeyl:
  ring s = 0,(z,x),dp;
  def S = nc_algebra(1,1); setring S;
  poly h = (x^2*z^2+x)*x;
  list fact2 = facFirstWeyl(h);
  map F = R,x,0,z;
  list fact1 = F(fact1); // it is identical to list fact2
  testNCfac(fact1); // check the correctness again
}
//==================================================

//==================================================
//************From here: Shift-Algebra**************
//==================================================
//==================================================*
//one factorization of a homogeneous polynomial
//in the first Shift Algebra
static proc homogfacFirstShift(poly h)
{//proc homogfacFirstShift
  int p=printlevel-voice+2; //for dbprint
  def r = basering;
  poly hath;
  intvec iv01 = intvec(0,1);
  int i; int j;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  if (!homogwithorder(h,iv01))
  {//The given polynomial is not homogeneous
    ERROR("The given polynomial is not homogeneous.");
    return(list());
  }//The given polynomial is not homogeneous
  if (h==0)
  {
    return(list(0));
  }
  list result;
  int m = deg(h,iv01);
  dbprint(p,dbprintWhitespace+" exclude the homogeneous part of deg. 0");
  if (m>0)
  {//The degree is not zero
    hath = lift(var(2)^m,h)[1,1];
    for (i = 1; i<=m;i++)
    {
      result = result + list(var(2));
    }
  }//The degree is not zero
  else
  {//The degree is zero
    hath = h;
  }//The degree is zero
  ring tempRing = 0,(x),dp;
  setring tempRing;
  map thetamap = r,x,1;
  poly hath = thetamap(hath);
  dbprint(p,dbprintWhitespace+" Factorize it using commutative factorization.");
  list azeroresult = factorize(hath);
  list azeroresult_return_form;
  for (i = 1; i<=size(azeroresult[1]);i++)
  {//rewrite the result of the commutative factorization
    for (j = 1; j <= azeroresult[2][i];j++)
    {
      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
    }
  }//rewrite the result of the commutative factorization
  setring(r);
  map finalmap = tempRing,var(1);
  list tempresult = finalmap(azeroresult_return_form);
  result = tempresult+result;
  return(result);
}//proc homogfacFirstShift

//==================================================
//Computes all possible homogeneous factorizations
static proc homogfacFirstShift_all(poly h)
{//proc HomogfacFirstShiftAll
  int p=printlevel-voice+2; //for dbprint
  intvec iv11 = intvec(1,1);
  if (deg(h,iv11) <= 0 )
  {//h is a constant
    return(list(list(h)));
  }//h is a constant
  def r = basering;
  list one_hom_fac; //stands for one homogeneous factorization
  int i; int j; int k;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  int shiftcounter;
  //Compute again a homogeneous factorization
  dbprint(p,dbprintWhitespace+" Computing one homog. factorization of the polynomial");
  one_hom_fac = homogfacFirstShift(h);
  one_hom_fac = delete(one_hom_fac,1);
  if (size(one_hom_fac) == 0)
  {//there is no homogeneous factorization or the polynomial was not homogeneous
    return(list());
  }//there is no homogeneous factorization or the polynomial was not homogeneous
  dbprint(p,dbprintWhitespace+" Permuting the 0-homogeneous part with the s");
  list result = permpp(one_hom_fac);
  for (i = 1; i<=size(result);i++)
  {
    shiftcounter = 0;
    for (j = 1; j<=size(result[i]); j++)
    {
      if (result[i][j]==var(2))
      {
        shiftcounter++;
      }
      else
      {
        result[i][j] = subst(result[i][j], var(1), var(1)-shiftcounter);
      }
    }
    result[i] = insert(result[i],1);
  }
  dbprint(p,dbprintWhitespace+" Deleting double entries in the resulting list");
  result = delete_dublicates_noteval(result);
  return(result);
}//proc HomogfacFirstShiftAll

//==================================================
//factorization of the first Shift Algebra
proc facFirstShift(poly h)
"USAGE: facFirstShift(h); h a polynomial in the first shift algebra
RETURN: list
PURPOSE: compute all factorizations of a polynomial in the first shift algebra
THEORY: Implements the new algorithm by A. Heinle and V. Levandovskyy, see the thesis of A. Heinle
ASSUME: basering is the first shift algebra
NOTE: Every entry of the output list is a list with factors for one possible factorization.
EXAMPLE: example facFirstShift; shows examples
SEE ALSO: testNCfac, facFirstWeyl, facSubWeyl
"{//facFirstShift
  int p = printlevel - voice + 2;
  int i;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  dbprint(p,dbprintWhitespace +" Checking the input.");
  if(nvars(basering)!=2)
  {//Our basering is the Shift algebra, but not the first
    ERROR("Basering is not the first shift algebra");
    return(list());
  }//Our basering is the Shift algebra, but not the first
  def r = basering;
  setring r;
  list LR = ringlist(r);
  number @n = leadcoef(LR[5][1,2]);
  poly @p = LR[6][1,2];
  if  ( @n!=number(1) )
  {
    ERROR("Basering is not the first shift algebra");
    return(list());
  }
  dbprint(p,dbprintWhitespace +" Done");
  list result = list();
  int j; int k; int l; //counter
  // create a ring with the ordering which makes shift algebra
  // graded
  // def r = basering; // done before
  ring tempRing = LR[1][1],(x,s),(a(0,1),Dp);
  def tempRingnc = nc_algebra(1,s);
  setring r;
  // information on relations
  if (@p == -var(1)) // reverted shift algebra
  {
    dbprint(p,dbprintWhitespace +" Reverted shift algebra. Swaping variables in Ringlist");
    setring(tempRingnc);
    map transf = r, var(2), var(1);
    setring(r);
    map transfback = tempRingnc, var(2),var(1);
    //    result = transfback(resulttemp);
  }
  else
  {
    if ( @p == var(2)) // usual shift algebra
    {
      setring(tempRingnc);
      map transf = r, var(1), var(2);
      //    result = facshift(h);
      setring(r);
      map transfback = tempRingnc, var(1),var(2);
    }
    else
    {
      ERROR("Basering is not the first shift algebra");
      return(list());
    }
  }
  // main calls
  setring(tempRingnc);
  dbprint(p,dbprintWhitespace +" Factorize the given polynomial with the subroutine sFacShift");
  list resulttemp = sFacShift(transf(h));
  dbprint(p,dbprintWhitespace +" Successful");
  setring(r);
  result = transfback(resulttemp);
  return( delete_dublicates_noteval(result) );
}//facFirstShift
example
{
  "EXAMPLE:";echo=2;
  ring R = 0,(x,s),dp;
  def r = nc_algebra(1,s);
  setring(r);
  poly h = (s^2*x+x)*s;
  facFirstShift(h);
}

static proc sFacShift(poly h)
"
USAGE: A static procedure to factorize a polynomial in the first Shift algebra, where all the
       validity checks were made in advance.
INPUT:  A polynomial h in the first Shift Algebra.
OUTPUT: A list of different factorizations of h, where the factors are irreducible
ASSUMPTIONS:
 - The basering is the first Shift algebra and has n as first, and s as second variable, i.e. we
   have var(2)*var(1) = var(1)*var(2)+1
THEORY: If the given polynomial h is [0,1]-homogeneous, the routines for homogeneous factorizations
 are called. Otherwise we map the polynomial into the first Weyl algebra (the first shift
 algebra is a subring of the first Weyl algebra), and use facFirstWeyl to factorize it. Later
 we map the factors back, if possible.
"
{//proc sFacShift
  int p = printlevel - voice + 2;
  int i; int j ;
  string dbprintWhitespace = "";
  number commonCoefficient = content(h);
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  //Checking if given polynomial is homogeneous
  if(homogwithorder(h,intvec(0,1)))
  {//The given polynomial is [0,1]-homogeneous
    dbprint(p,dbprintWhitespace+"The polynomial is [0,1]-homogeneous. Returning the
homogeneous factorization");
    return(homogfacFirstShift_all(h));
  }//The given polynomial is [0,1]-homogeneous

  //---------- Start of interesting part ----------

  dbprint(p,dbprintWhitespace+"Mapping the polynomial h into the first Weyl algebra.");
  poly temph = h/commonCoefficient;
  def ourBaseRing = basering;
  ring tempWeylAlgebraComm = 0,(x,d),dp;
  def tempWeylAlgebra = nc_algebra(1,1);
  setring(tempWeylAlgebra);
  map shiftMap = ourBaseRing, x*d, d;
  poly h = shiftMap(temph);
  dbprint(p,dbprintWhitespace+"Successful! The polynomial in the Weyl algebra is "+string(h));
  dbprint(p,dbprintWhitespace+"Factorizing the polynomial in the first Weyl algebra");
  list factorizationInWeyl = facFirstWeyl(h);
  dbprint(p,dbprintWhitespace+"Successful! The factorization is given by:");
  dbprint(p,factorizationInWeyl);
  list validCombinations;

  dbprint(p,dbprintWhitespace+"Now we will map this back to the shift algebra and filter
valid results");
  //-Now we map the results back to the shift algebra. But first, we need to combine them properly.
  for (i = 1; i<=size(factorizationInWeyl); i++)
  {//Deleting the first Coefficient factor
    factorizationInWeyl[i] = delete(factorizationInWeyl[i],1);
    validCombinations = validCombinations + combineNonnegative(factorizationInWeyl[i]);
  }//Deleting the first Coefficient factor
  if (size(validCombinations) == 0)
  {//There are no valid combinations, therefore we can directly say, that h is irreducible
    setring(ourBaseRing);
    return(list(list(commonCoefficient, h/commonCoefficient)));
  }//There are no valid combinations, therefore we can directly say, that h is irreducible
  validCombinations = delete_dublicates_noteval(validCombinations);
  setring(ourBaseRing);
  map backFromWeyl = tempWeylAlgebra, var(1),var(2);
  list validCombinations = backFromWeyl(validCombinations);
  for (i = 1; i<=size(validCombinations); i++)
  {
    for (j = 1; j<=size(validCombinations[i]);j++)
    {
      setring(tempWeylAlgebra);
      fromWeylToShiftPoly(validCombinations[i][j],ourBaseRing);
      validCombinations[i][j] = result;
      kill result;
      kill tempResult;
      kill zeroPoly;
      kill fromWeyl;
    }
  }
  for (i = 1; i<=size(validCombinations); i++)
  {//Adding the common factor in the first position of the list
    validCombinations[i] = insert(validCombinations[i],commonCoefficient);
  }//Adding the common factor in the first position of the list
  dbprint(dbprintWhitespace+"Done.");
  //mapping
  return(validCombinations);
}//proc sFacShift

static proc combineNonnegative(list l)
"
USAGE: In sFacShift, when we want to map back the results of the factorization of the polynomial in
       the first Weyl algebra to the shift algebra. We need to recombine the factors such that
       we can map it back to the shift algebra without any problems.
INPUT: A list l containing one factorization of a polynomial in the first Weyl algebra. For example
       for the polynomial (1+x)*(1+x+d) we would have the list [1,x+1,x+d+1].
OUTPUT:If we can map every factor without a problem back to the shift algebra (i.e. if the smallest
 homogeneous summand of every factor is of nonnegative degree), a list containing the same
 list as given in the input is returned.
 If otherwise some factors cause problems, we consider every possible combination (i.e.
 products of the factors) and extract those where all factors have a smallest homogeneous
 summand of nonnegative degree.
ASSUMPTIONS:
 - Weyl algebra is given, and we have var(2)*var(1)=var(1)*var(2) +1
"
{//combineNonnegative
  int p = printlevel - voice + 2;
  int i;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  //First the easy case: all of the factors fulfill the condition of mapping to shift:
  dbprint(p,dbprintWhitespace+"Checking, if the given factors
can already be mapped without a problem.");
  int isValid = 1;
  for (i = 1; i<=size(l);i++)
  {//Checking for every entry if the condition is fulfilled.
    if (deg(l[i],intvec(1,-1))>0)
    {//Found one, where it is not fulfilled
      isValid = 0;
      break;
    }//Found one, where it is not fulfilled
  }//Checking for every entry if the condition is fulfilled.
  dbprint(p,dbprintWhitespace+"Done.");
  if (isValid)
  {//We can map every factor to the shift algebra and do not need to combine anything
    dbprint(p,dbprintWhitespace+"They can be mapped. Therefore we return them directly.");
    return(list(l));
  }//We can map every factor to the shift algebra and do not need to combine anything
  dbprint(p,dbprintWhitespace+"They cannot be mapped. Looking for valid combinations.");
  //Starting with the case, where l only consists of 1 or two elements.
  if(size(l)<=2)
  {//The case where we won't call the function a second time
    if (deg(product(l),intvec(1,-1))>0)
    {//No way of a valid combination
      return(list());
    }//No way of a valid combination
    else
    {//The product is the only possible and valid combination
      return(list(list(product(l))));
    }//The product is the only possible and valid combination
  }//The case where we won't call the function a second time
  //---------- Easy pre-stuff done. now we combine the factors.----------
  int pos;
  int j; int k;
  dbprint(p,dbprintWhitespace+"Making combinations of two.");
  list combinationsOfTwo = combinekfinlf(l,2);
  dbprint(p,dbprintWhitespace+"Done. Now checking, if there are valid ones in between.");
  list result;
  list validLHS;
  list validRHS;
  for (i = 1; i<=size(combinationsOfTwo); i++)
  {//go through all combinations and detect the valid ones
    if(deg(combinationsOfTwo[i][1],intvec(1,-1))>0 or deg(combinationsOfTwo[i][2],intvec(1,-1))>0)
    {//No chance, so no further treatment needed
      i++;
      continue;
    }//No chance, so no further treatment needed
    for (pos = 1; pos<=size(l);pos++)
    {//find the position where the combination splits
      if (product(l[1..pos]) == combinationsOfTwo[i][1])
      {//Found the position
        break;
      }//Found the position
    }//find the position where the combination splits
    dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " +
            string(list(l[1..pos])));
    validLHS = combineNonnegative(list(l[1..pos]));
    dbprint(p,dbprintWhitespace+"Calling combineNonnegative recursively with argument " +
            string(list(l[pos+1..size(l)])));
    validRHS = combineNonnegative(list(l[pos+1..size(l)]));
    for (j = 1; j<=size(validLHS); j++)
    {//Combining the left hand side valid combnations...
      for (k = 1; k<=size(validRHS); k++)
      {//... with the right hand side valid combinations
        result = insert(result, validLHS[j]+validRHS[k]);
      }//... with the right hand side valid combinations
    }//Combining the left hand side valid combnations...
  }//go through all combinations and detect the valid ones
  result = delete_dublicates_noteval(result);
  dbprint(p,dbprintWhitespace+"Done.");
  return(result);
}//combineNonnegative

static proc fromWeylToShiftPoly(poly h,def sAlgebra)
"
USAGE: Given a polynomial in the first Weyl algebra, this method returns it -- if possible --
       as an element in the first shift algebra, which is given in the method header.
INPUT: A polynomial h, and the first shift algebra as a ring
OUTPUT: The correct mapping in the shift Algebra
ASSUMPTIONS:
 - The lowest [-1,1]-homogeneous summand of h is of nonnegative degree
 - The shift algebra is given in the way that var(2)*var(1) = (var(1)+1)*var(2)
"
{//fromWeylToShiftPoly
  int p = printlevel - voice + 2;
  int i;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  if (deg(h,intvec(1,-1))>0)
  {//Wrong input polynomial
    ERROR("The lowest [-1,1] homogeneous summand of "+string(h)+" is of negative degree.");
  }//Wrong input polynomial
  def ourHomeBase = basering;
  list hDist = homogDistribution(h);
  setring(sAlgebra);
  poly result = 0;
  poly tempResult;
  poly zeroPoly;
  map fromWeyl = ourHomeBase, var(1), var(2);
  setring(ourHomeBase);
  poly zeroPoly;
  poly tempZeroPoly;
  int j; int k;
  int derDeg;
  for (i = 1; i<=size(hDist);i++)
  {
    derDeg = hDist[i][1];
    setring(sAlgebra);
    tempResult = 1;
    setring(ourHomeBase);
    zeroPoly = lift(d^derDeg, hDist[i][2])[1,1];
    for (j = 1; j<=size(zeroPoly); j++)
    {
      tempZeroPoly = zeroPoly[j];
      setring(sAlgebra);
      zeroPoly = fromWeyl(tempZeroPoly);
      tempResult = tempResult * leadcoef(zeroPoly);
      setring(ourHomeBase);
      for (k = 1; k<=deg(zeroPoly[j],intvec(0,1));k++)
      {
        setring(sAlgebra);
        tempResult = tempResult*(var(1)-(k-1));
        setring(ourHomeBase);
      }
      setring(sAlgebra);
      result = result + tempResult*var(2)^derDeg;
      tempResult = 1;
      setring(ourHomeBase);
    }
  }
  setring(sAlgebra);
  keepring(sAlgebra);
}//fromWeylToShiftPoly

static proc refineFactList(list L)
{
  // assume: list L is an output of factorization proc
  // doing: remove doubled entries
  int s = size(L); int sm;
  int i,j,k,cnt;
  list M, U, A, B;
  A = L;
  k = 0;
  cnt  = 1;
  for (i=1; i<=s; i++)
  {
    if (size(A[i]) != 0)
    {
      M = A[i];
      //      "probing with"; M; i;
      B[cnt] = M; cnt++;
      for (j=i+1; j<=s; j++)
      {
        if ( isEqualList(M,A[j]) )
        {
          k++;
          // U consists of intvecs with equal pairs
          U[k] = intvec(i,j);
          A[j] = 0;
        }
      }
    }
  }
  kill A,U,M;
  return(B);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = 0,(x,s),dp;
  def r = nc_algebra(1,1);
  setring(r);
  list l,m;
  l = list(1,s2+1,x,s,x+s);
  m = l,list(1,s,x,s,x),l;
  refineFactList(m);
}

static proc isEqualList(list L, list M)
{
  // int boolean: 1=yes, 0 =no : test whether two lists are identical
  int s = size(L);
  if (size(M)!=s) { return(0); }
  int j=1;
  while ( (L[j]==M[j]) && (j<s) )
  {
    j++;
  }
  if (L[j]==M[j])
  {
    return(1);
  }
  return(0);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = 0,(x,s),dp;
  def r = nc_algebra(1,1);
  setring(r);
  list l,m;
  l = list(1,s2+1,x,s,x+s);
  m = l;
  isEqualList(m,l);
}


//////////////////////////////////////////////////
// Q-WEYL-SECTION
//////////////////////////////////////////////////

//==================================================
//A function to get the i'th triangular number
static proc triangNum(int n)
{
  if (n == 0)
  {
    return(0);
  }
  return (n*(n+1) div 2);
}

//==================================================*
//one factorization of a homogeneous polynomial
//in the first Q Weyl Algebra
proc homogfacFirstQWeyl(poly h)
"USAGE: homogfacFirstQWeyl(h); h is a homogeneous polynomial in the
 first q-Weyl algebra with respect to the weight vector [-1,1]
RETURN: list
PURPOSE: Computes a factorization of a homogeneous polynomial h with
  respect to the weight vector [-1,1] in the first q-Weyl algebra
THEORY: @code{homogfacFirstQWeyl} returns a list with a factorization of the given,
 [-1,1]-homogeneous polynomial. If the degree of the polynomial is k with
 k positive, the last k entries in the output list are the second
 variable. If k is positive, the last k entries will be x. The other
 entries will be irreducible polynomials of degree zero or 1 resp. -1.
SEE ALSO: homogfacFirstQWeyl_all
"{//proc homogfacFirstQWeyl
  int p = printlevel-voice+2;//for dbprint
  def r = basering;
  poly hath;
  int i; int j;
  string dbprintWhitespace = "";
  for (i = 1; i<=voice;i++)
  {dbprintWhitespace = dbprintWhitespace + " ";}
  intvec ivm11 = intvec(-1,1);
  if (!homogwithorder(h,ivm11))
  {//The given polynomial is not homogeneous
    ERROR("Given polynomial was not [-1,1]-homogeneous");
    return(list());
  }//The given polynomial is not homogeneous
  if (h==0)
  {
    return(list(0));
  }
  list result;
  int m = deg(h,ivm11);
  dbprint(p,dbprintWhitespace+" Splitting the polynomial in A_0 and A_k-Part");
  if (m!=0)
  {//The degree is not zero
    if (m <0)
    {//There are more x than y
      hath = lift(var(1)^(-m),h)[1,1];
      for (i = 1; i<=-m; i++)
      {
        result = result + list(var(1));
      }
    }//There are more x than y
    else
    {//There are more y than x
      hath = lift(var(2)^m,h)[1,1];
      for (i = 1; i<=m;i++)
      {
        result = result + list(var(2));
      }
    }//There are more y than x
  }//The degree is not zero
  else
  {//The degree is zero
    hath = h;
  }//The degree is zero
  dbprint(p," Done");
  //beginning to transform x^i*y^i in theta(theta-1)...(theta-i+1)
  list mons;
  dbprint(p," Putting the monomials in the A_0-part in a list.");
  for(i = 1; i<=size(hath);i++)
  {//Putting the monomials in a list
    mons = mons+list(hath[i]);
  }//Putting the monomials in a list
  dbprint(p," Done");
  dbprint(p," Mapping this monomials to K(q)[theta]");
  //Now, map to the commutative ring with theta:
  list tempRingList = ringlist(r);
  tempRingList[2] = insert(tempRingList[2],"theta",2); //New variable theta = x*d
  tempRingList = delete(tempRingList,5);
  tempRingList = delete(tempRingList,5); //The ring should now be commutative
  def tempRing = ring(tempRingList);
  setring tempRing;
  map thetamap = r,var(1),var(2);
  list mons = thetamap(mons);
  poly entry;
  poly tempSummand;
  for (i = 1; i<=size(mons);i++)
  {//transforming the monomials as monomials in theta
    entry = 1;//leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
    for (j = 0; j<leadexp(mons[i])[2];j++)
    {
      tempSummand = (par(1)^j-1)/(par(1)-1);
      entry = entry * theta-tempSummand*entry;
    }
    //entry;
    //leadcoef(mons[i]) * q^(-triangNum(leadexp(mons[i])[2]-1));
    mons[i] = entry*leadcoef(mons[i]) * par(1)^(-triangNum(leadexp(mons[i])[2]-1));
  }//transforming the monomials as monomials in theta
  dbprint(p," Done");
  dbprint(p," Factorize the A_0-Part in K[theta]");
  list azeroresult = factorize(sum(mons));
  dbprint(p," Successful");
  list azeroresult_return_form;
  for (i = 1; i<=size(azeroresult[1]);i++)
  {//rewrite the result of the commutative factorization
    for (j = 1; j <= azeroresult[2][i];j++)
    {
      azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
    }
  }//rewrite the result of the commutative factorization
  dbprint(p," Mapping back to A_0.");
  setring(r);
  map finalmap = tempRing,var(1),var(2),var(1)*var(2);
  list tempresult = finalmap(azeroresult_return_form);
  dbprint(p,"Successful.");
  for (i = 1; i<=size(tempresult);i++)
  {//factorizations of theta resp. theta +1
    if(tempresult[i]==var(1)*var(2))
    {
      tempresult = insert(tempresult,var(1),i-1);
      i++;
      tempresult[i]=var(2);
    }
    if(tempresult[i]==var(2)*var(1))
    {
      tempresult = insert(tempresult,var(2),i-1);
      i++;
      tempresult[i]=var(1);
    }
  }//factorizations of theta resp. theta +1
  result = tempresult+result;
  //Correction of the result in the special q-Case:
  for (j = 2 ; j<= size(result);j++)
  {//Div the whole Term by the leading coefficient and multiply it to the first entry in result[i]
    result[1] = result[1] * leadcoef(result[j]);
    result[j] = 1/leadcoef(result[j]) * result[j];
  }//Div the whole Term by the leading coefficient and multiply it to the first entry in result[i]
  return(result);
}//proc homogfacFirstQWeyl
example
{
  "EXAMPLE:";echo=2;
  ring R = (0,q),(x,d),dp;
  def r = nc_algebra (q,1);
  setring(r);
  poly h = q^25*x^10*d^10+q^16*(q^4+q^3+q^2+q+1)^2*x^9*d^9+
    q^9*(q^13+3*q^12+7*q^11+13*q^10+20*q^9+26*q^8+30*q^7+
    31*q^6+26*q^5+20*q^4+13*q^3+7*q^2+3*q+1)*x^8*d^8+
    q^4*(q^9+2*q^8+4*q^7+6*q^6+7*q^5+8*q^4+6*q^3+
     4*q^2+2q+1)*(q^4+q^3+q^2+q+1)*(q^2+q+1)*x^7*d^7+
    q*(q^2+q+1)*(q^5+2*q^4+2*q^3+3*q^2+2*q+1)*(q^4+q^3+q^2+q+1)*(q^2+1)*(q+1)*x^6*d^6+
    (q^10+5*q^9+12*q^8+21*q^7+29*q^6+33*q^5+31*q^4+24*q^3+15*q^2+7*q+12)*x^5*d^5+
    6*x^3*d^3+24;
  homogfacFirstQWeyl(h);
}

//==================================================
//Computes all possible homogeneous factorizations for an element in the first Q-Weyl Algebra
proc homogfacFirstQWeyl_all(poly h)
"USAGE: homogfacFirstQWeyl_all(h); h is a homogeneous polynomial in the first q-Weyl algebra
 with respect to the weight vector [-1,1]
RETURN: list
PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect
  to the weight vector [-1,1] in the first q-Weyl algebra
THEORY: @code{homogfacFirstQWeyl} returns a list with all factorization of the given,
 homogeneous polynomial. It uses the output of homogfacFirstQWeyl and permutes
 its entries with respect to the commutation rule. Furthermore, if a
 factor of degree zero is irreducible in K[  heta], but reducible in
 the first q-Weyl algebra, the permutations of this element with the other
 entries will also be computed.
SEE ALSO: homogfacFirstQWeyl
"{//proc HomogfacFirstQWeylAll
  int p=printlevel-voice+2;//for dbprint
  intvec iv11= intvec(1,1);
  if (deg(h,iv11) <= 0 )
  {//h is a constant
    dbprint(p,"Given polynomial was not homogeneous");
    return(list(list(h)));
  }//h is a constant
  def r = basering;
  list one_hom_fac; //stands for one homogeneous factorization
  int i; int j; int k;
  intvec ivm11 = intvec(-1,1);
  dbprint(p," Calculate one homogeneous factorization using homogfacFirstQWeyl");
  //Compute again a homogeneous factorization
  one_hom_fac = homogfacFirstQWeyl(h);
  dbprint(p,"Successful");
  if (size(one_hom_fac) == 0)
  {//there is no homogeneous factorization or the polynomial was not homogeneous
    return(list());
  }//there is no homogeneous factorization or the polynomial was not homogeneous
  //divide list in A0-Part and a list of x's resp. y's
  list list_not_azero = list();
  list list_azero;
  list k_factor;
  int is_list_not_azero_empty = 1;
  int is_list_azero_empty = 1;
  k_factor = list(one_hom_fac[1]);
  if (absValue(deg(h,ivm11))<size(one_hom_fac)-1)
  {//There is a nontrivial A_0-part
    list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,ivm11)))];
    is_list_azero_empty = 0;
  }//There is a nontrivial A_0 part
  dbprint(p," Combine x,y to xy in the factorization again.");
  for (i = 1; i<=size(list_azero)-1;i++)
  {//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
    if (list_azero[i] == var(1))
    {
      if (list_azero[i+1]==var(2))
      {
        list_azero[i] = var(1)*var(2);
        list_azero = delete(list_azero,i+1);
      }
    }
    if (list_azero[i] == var(2))
    {
      if (list_azero[i+1]==var(1))
      {
        list_azero[i] = var(2)*var(1);
        list_azero = delete(list_azero,i+1);
      }
    }
  }//in homogfacFirstQWeyl, we factorized theta, and this will be made undone
  dbprint(p," Done");
  if(deg(h,ivm11)!=0)
  {//list_not_azero is not empty
    list_not_azero =
      one_hom_fac[(size(one_hom_fac)-absValue(deg(h,ivm11))+1)..size(one_hom_fac)];
    is_list_not_azero_empty = 0;
  }//list_not_azero is not empty
  //Map list_azero in K[theta]
  dbprint(p," Map list_azero to K[theta]");
  //Now, map to the commutative ring with theta:
  list tempRingList = ringlist(r);
  tempRingList[2] = insert(tempRingList[2],"theta",2); //New variable theta = x*d
  tempRingList = delete(tempRingList,5);
  tempRingList = delete(tempRingList,5); //The ring should now be commutative
  def tempRing = ring(tempRingList);
  setring(tempRing);
  poly entry;
  map thetamap = r,var(1),var(2);
  if(!is_list_not_azero_empty)
  {//Mapping in Singular is only possible, if the list before
    //contained at least one element of the other ring
    list list_not_azero = thetamap(list_not_azero);
  }//Mapping in Singular is only possible, if the list before
  //contained at least one element of the other ring
  if(!is_list_azero_empty)
  {//Mapping in Singular is only possible, if the list before
    //contained at least one element of the other ring
    list list_azero= thetamap(list_azero);
  }//Mapping in Singular is only possible, if the list before
  //contained at least one element of the other ring
  list k_factor = thetamap(k_factor);
  list tempmons;
  dbprint(p," Done");
  for(i = 1; i<=size(list_azero);i++)
  {//rewrite the polynomials in A1 as polynomials in K[theta]
    tempmons = list();
    for (j = 1; j<=size(list_azero[i]);j++)
    {
      tempmons = tempmons + list(list_azero[i][j]);
    }
    for (j = 1 ; j<=size(tempmons);j++)
    {
      //entry = leadcoef(tempmons[j]);
      entry = leadcoef(tempmons[j]) * par(1)^(-triangNum(leadexp(tempmons[j])[2]-1));
      for (k = 0; k < leadexp(tempmons[j])[2];k++)
      {
        entry = entry*(theta-(par(1)^k-1)/(par(1)-1));
      }
      tempmons[j] = entry;
    }
    list_azero[i] = sum(tempmons);
  }//rewrite the polynomials in A1 as polynomials in K[theta]
  //Compute all permutations of the A0-part
  dbprint(p," Compute all permutations of the A_0-part with the first resp. the snd. variable");
  list result;
  int shift_sign;
  int shift;
  poly shiftvar;
  if (size(list_not_azero)!=0)
  {//Compute all possibilities to permute the x's resp. the y's in the list
    if (list_not_azero[1] == var(1))
    {//h had a negative weighted degree
      shift_sign = 1;
      shiftvar = var(1);
    }//h had a negative weighted degree
    else
    {//h had a positive weighted degree
      shift_sign = -1;
      shiftvar = var(2);
    }//h had a positive weighted degree
    result = permpp(list_azero + list_not_azero);
    for (i = 1; i<= size(result); i++)
    {//adjust the a_0-parts
      shift = 0;
      for (j=1; j<=size(result[i]);j++)
      {
        if (result[i][j]==shiftvar)
        {
          shift = shift + shift_sign;
        }
        else
        {
          if (shift < 0)
          {//We have two distict formulas for x and y. In this case use formula for y
            if (shift == -1)
            {
              result[i][j] = subst(result[i][j],theta,1/par(1)*(theta - 1));
            }
            else
            {
              result[i][j] =
                subst(result[i][j],
                      theta,
                      1/par(1)*((theta - 1)/par(1)^(absValue(shift)-1)
                                - (par(1)^(shift +2)-par(1))/(1-par(1))));
            }
          }//We have two distict formulas for x and y. In this case use formula for y
          if (shift > 0)
          {//We have two distict formulas for x and y. In this case use formula for x
            if (shift == 1)
            {
              result[i][j] = subst(result[i][j],theta,par(1)*theta + 1);
            }
            else
            {
              result[i][j] =
                subst(result[i][j],
                      theta,par(1)^shift*theta+(par(1)^shift-1)/(par(1)-1));
            }
          }//We have two distict formulas for x and y. In this case use formula for x
        }
      }
    }//adjust the a_0-parts
  }//Compute all possibilities to permute the x's resp. the y's in the list
  else
  {//The result is just all the permutations of the a_0-part
    result = permpp(list_azero);
  }//The result is just all the permutations of the a_0 part
  if (size(result)==0)
  {
    return(result);
  }
  dbprint(p," Done");
  dbprint(p," Searching for theta resp. theta + 1 in the list and factorize them");
  //Now we are going deeper and search for theta resp. theta + 1, substitute
  //them by xy resp. yx and go on permuting
  int found_theta;
  int thetapos;
  list leftpart;
  list rightpart;
  list lparts;
  list rparts;
  list tempadd;
  for (i = 1; i<=size(result) ; i++)
  {//checking every entry of result for theta or theta +1
    found_theta = 0;
    for(j=1;j<=size(result[i]);j++)
    {
      if (result[i][j]==theta)
      {//the jth entry is theta and can be written as x*y
        thetapos = j;
        result[i]= insert(result[i],var(1),j-1);
        j++;
        result[i][j] = var(2);
        found_theta = 1;
        break;
      }//the jth entry is theta and can be written as x*y
      if(result[i][j] == par(1)*theta +1)
      {
        thetapos = j;
        result[i] = insert(result[i],var(2),j-1);
        j++;
        result[i][j] = var(1);
        found_theta = 1;
        break;
      }
    }
    if (found_theta)
    {//One entry was theta resp. theta +1
      leftpart = result[i];
      leftpart = leftpart[1..thetapos];
      rightpart = result[i];
      rightpart = rightpart[(thetapos+1)..size(rightpart)];
      lparts = list(leftpart);
      rparts = list(rightpart);
      //first deal with the left part
      if (leftpart[thetapos] == var(1))
      {
        shift_sign = 1;
        shiftvar = var(1);
      }
      else
      {
        shift_sign = -1;
        shiftvar = var(2);
      }
      for (j = size(leftpart); j>1;j--)
      {//drip x resp. y
        if (leftpart[j-1]==shiftvar)
        {//commutative
          j--;
          continue;
        }//commutative
        if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
        {//stop here
          break;
        }//stop here
        //Here, we can only have a a0- part
        if (shift_sign<0)
        {
          leftpart[j] = subst(leftpart[j-1],theta, 1/par(1)*(theta +shift_sign));
        }
        if (shift_sign>0)
        {
          leftpart[j] = subst(leftpart[j-1],theta, par(1)*theta + shift_sign);
        }
        leftpart[j-1] = shiftvar;
        lparts = lparts + list(leftpart);
      }//drip x resp. y
      //and now deal with the right part
      if (rightpart[1] == var(1))
      {
        shift_sign = 1;
        shiftvar = var(1);
      }
      else
      {
        shift_sign = -1;
        shiftvar = var(2);
      }
      for (j = 1 ; j < size(rightpart); j++)
      {
        if (rightpart[j+1] == shiftvar)
        {
          j++;
          continue;
        }
        if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
        {
          break;
        }
        if (shift_sign<0)
        {
          rightpart[j] = subst(rightpart[j+1], theta, par(1)*theta - shift_sign);
        }
        if (shift_sign>0)
        {
          rightpart[j] = subst(rightpart[j+1], theta, 1/par(1)*(theta - shift_sign));
        }
        rightpart[j+1] = shiftvar;
        rparts = rparts + list(rightpart);
      }
      //And now, we put all possibilities together
      tempadd = list();
      for (j = 1; j<=size(lparts); j++)
      {
        for (k = 1; k<=size(rparts);k++)
        {
          tempadd = tempadd + list(lparts[j]+rparts[k]);
        }
      }
      tempadd = delete(tempadd,1); // The first entry is already in the list
      result = result + tempadd;
      continue; //We can may be not be done already with the ith entry
    }//One entry was theta resp. theta +1
  }//checking every entry of result for theta or theta +1
  dbprint(p," Done");
  //map back to the basering
  dbprint(p," Mapping back everything to the basering");
  setring(r);
  map finalmap = tempRing, var(1), var(2),var(1)*var(2);
  list result = finalmap(result);
  for (i=1; i<=size(result);i++)
  {//adding the K factor
    result[i] = k_factor + result[i];
  }//adding the k-factor
  dbprint(p," Done");
  dbprint(p," Delete double entries in the list.");
  result = delete_dublicates_noteval(result);
  dbprint(p," Done");
  return(result);
}//proc HomogfacFirstQWeylAll
example
{
  "EXAMPLE:";echo=2;
  ring R = (0,q),(x,d),dp;
  def r = nc_algebra (q,1);
  setring(r);
  poly h = q^25*x^10*d^10+q^16*(q^4+q^3+q^2+q+1)^2*x^9*d^9+
    q^9*(q^13+3*q^12+7*q^11+13*q^10+20*q^9+26*q^8+30*q^7+
    31*q^6+26*q^5+20*q^4+13*q^3+7*q^2+3*q+1)*x^8*d^8+
    q^4*(q^9+2*q^8+4*q^7+6*q^6+7*q^5+8*q^4+6*q^3+
    4*q^2+2q+1)*(q^4+q^3+q^2+q+1)*(q^2+q+1)*x^7*d^7+
    q*(q^2+q+1)*(q^5+2*q^4+2*q^3+3*q^2+2*q+1)*(q^4+q^3+q^2+q+1)*(q^2+1)*(q+1)*x^6*d^6+
    (q^10+5*q^9+12*q^8+21*q^7+29*q^6+33*q^5+31*q^4+24*q^3+15*q^2+7*q+12)*x^5*d^5+
    6*x^3*d^3+24;
  homogfacFirstQWeyl_all(h);
}

//TODO: FirstQWeyl check the parameters...

//==================================================
// EASY EXAMPLES FOR WEYL ALGEBRA
//==================================================
/*
  Easy and fast example polynomials where one can find factorizations: K<x,d |dx=xd+1>
  (x^2+d)*(x^2+d);
  (x^2+x)*(x^2+d);
  (x^3+x+1)*(x^4+d*x+2);
  (x^2*d+d)*(d+x*d);
  d^3+x*d^3+2*d^2+2*(x+1)*d^2+d+(x+2)*d; //Example 5 Grigoriev-Schwarz.
  (d+1)*(d+1)*(d+x*d); //Landau Example projected to the first dimension.
*/

//==================================================
//Some Bugs(fixed)/hard examples from Martin Lee:
//==================================================
// ex1, ex2
/*
ring s = 0,(x,d),Ws(-1,1);
def S = nc_algebra(1,1); setring S;
poly a = 10x5d4+26x4d5+47x5d2-97x4d3; //Not so hard any more... Done in around 4 minutes
def l= facFirstWeyl (a); l;
kill l;
poly b = -5328x8d5-5328x7d6+720x9d2+720x8d3-16976x7d4-38880x6d5
-5184x7d3-5184x6d4-3774x5d5+2080x8d+5760x7d2-6144x6d3-59616x5d4
+3108x3d6-4098x6d2-25704x5d3-21186x4d4+8640x6d-17916x4d3+22680x2d5
+2040x5d-4848x4d2-9792x3d3+3024x2d4-10704x3d2-3519x2d3+34776xd4
+12096xd3+2898d4-5040x2d+8064d3+6048d2; //Still very hard... But it seems to be only because of the
//combinatorial explosion
def l= facFirstWeyl (b); l;

// ex3: there was difference in answers => fixed
LIB "ncfactor.lib";
ring r = 0,(x,y,z),dp;
matrix D[3][3]; D[1,3]=-1;
def R = nc_algebra(1,D);
setring R;
poly g= 7*z4*x+62*z3+26*z;
def l1= facSubWeyl (g, x, z);
l1;
//---- other ring
ring s = 0,(x,z),dp;
def S = nc_algebra(1,-1); setring S;
poly g= 7*z4*x+62*z3+26*z;
def l2= facFirstWeyl (g);
l2;
map F = R,x,0,z;
list l1 = F(l1);
l1;
//---- so the answers look different, check them!
testNCfac(l2); // ok
testNCfac(l1); // was not ok, but now it's been fixed!!!

// selbst D und X so vertauschen dass sie erfuellt ist : ist gemacht

*/

/*
// bug from M Lee
LIB "ncfactor.lib";
ring s = 0,(z,x),dp;
def S = nc_algebra(1,1); setring S;
poly f= -60z4x2-54z4-56zx3-59z2x-64;
def l= facFirstWeyl (f);
l; // before: empty list; after fix: 1 entry, f is irreducible
poly g = 75z3x2+92z3+24;
def l= facFirstWeyl (g);
l; //before: empty list, now: correct
*/

/* more things from Martin Lee; fixed
ring R = 0,(x,s),dp;
def r = nc_algebra(1,s);
setring(r);
poly h = (s2*x+x)*s;
h= h* (x+s);
def l= facFirstShift(h);
l; // contained doubled entries: not anymore, fixed!

ring R = 0,(x,s),dp;
def r = nc_algebra(1,-1);
setring(r);
poly h = (s2*x+x)*s;
h= h* (x+s);
def l= facFirstWeyl(h);
l; // contained doubled entries: not anymore, fixed!

*/

//======================================================================
//Examples from TestSuite that are terminating in a reasonable time.
//======================================================================

//Counter example for old Algorithm, but now working:
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (1+x^2*d)^4;
list lsng = facFirstWeyl(h);
print(lsng);
*/

//Example 2.7. from Master thesis
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (xdd + xd+1+ (xd+5)*x)*(((x*d)^2+1)*d + xd+3+ (xd+7)*x);
list lsng = facFirstWeyl(h);
print(lsng);
 */

//Example with high combinatorial income
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (xdddd + (xd+1)*d*d+ (xd+5)*x*d*d)*(((x*d)^2+1)*d*x*x + (xd+3)*x*x+ (xd+7)*x*x*x);
list lsng = facFirstWeyl(h);
print(lsng);
 */

//Once a bug, now working
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (x^2*d^2+x)*(x+1);
list lsng = facFirstWeyl(h);
print(lsng);
*/

//Another one of that kind
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (x*d*d + (x*d)^5 +x)*((x*d+1)*d-(x*d-1)^5+x);
list lsng = facFirstWeyl(h);
print(lsng);
 */

//Example of Victor for Shift Algebra
/*
ring s = 0,(n,Sn),dp;
def S = nc_algebra(1,Sn); setring S;
LIB "ncfactor.lib";
list lsng = facFirstShift(n^2*Sn^2+3*n*Sn^2-n^2+2*Sn^2-3*n-2);
print(lsng);
 */

//Interesting example, as there are actually also some complex solutions to it:
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "/Users/albertheinle/Studium/forschung/ncfactor/versionen/ncfactor.lib";
poly h =(x^3+x+1)*(x^4+d*x+2);//Example for finitely many, but more than one solution in between.
list lsng = facFirstWeyl(h);
print(lsng);
 */

//Another one of that kind:
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h =(x^2+d)*(x^2+d);//Example for finitely many, but more than one solution in between.
list lsng = facFirstWeyl(h);
print(lsng);
*/

//Example by W. Koepf:
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (x^4-1)*x*d^2+(1+7*x^4)*d+8*x^3;
list lsng = facFirstWeyl(h);
print(lsng);
 */

//Shift Example from W. Koepf
/*
ring R = 0,(n,s),dp;
def r = nc_algebra(1,s);
setring(r);
LIB "ncfactor.lib";
poly h = n*(n+1)*s^2-2*n*(n+100)*s+(n+99)*(n+100);
list lsng = facFirstShift(h);
print(lsng);
 */

//Tsai Example... Once hard, now easy...
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (x^6+2*x^4-3*x^2)*d^2-(4*x^5-4*x^4-12*x^2-12*x)*d + (6*x^4-12*x^3-6*x^2-24*x-12);
list lsng =facFirstWeyl(h);
print(lsng);
 */

//======================================================================
// Hard examples not yet calculatable in feasible amount of time
//======================================================================

//Also a counterexample for REDUCE. Very long Groebner basis computation in between.
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
LIB "ncfactor.lib";
poly h = (d^4+x^2+dx+x)*(d^2+x^4+xd+d);
list lsng = facFirstWeyl(h);
print(lsng);
*/

//Example from the Mainz-Group
/*
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
poly dop6 = 1/35*x^4*(27-70*x+35*x^2)+ 1/35*x*(32+152*x+100*x^2-59*x^3+210*x^4+105*x^5)*d+
(-10368/35-67056/35*x-35512/7*x^2-50328/7*x^3-40240/7*x^4-2400*x^5-400*x^6)*d^2+
(-144/35*(x+1)*(1225*x^5+11025*x^4+37485*x^3+61335*x^2+50138*x+16584)-6912/35*(x+2)*
(x+1)*(105*x^4+1155*x^3+4456*x^2+7150*x+4212) -27648/35*(x+3)*(x+1)*(35*x^2+350*x+867)*
(x+2)^2)*d^3;
LIB "ncfactor.lib";
printlevel = 5;
facFirstWeyl(dop6);
$;*/

//Another Mainz Example:
/*
LIB "ncfactor.lib";
ring R = 0,(x,d),dp;
def r = nc_algebra(1,1);
setring(r);
poly dopp = 82547*x^4*d^4+60237*x^3*d^3+26772*x^5*d^5+2231*x^6*d^6+x*(1140138*
x^2*d^2-55872*x*d-3959658*x^3*d^3-8381805*x^4*d^4-3089576*x^5*d^5-274786*
x^6*d^6)+x^2*(-16658622*x*d-83427714*x^2*d^2-19715033*x^3*d^3+78915395*x^4
*d^4+35337930*x^5*d^5+3354194*x^6*d^6)+x^3*(-99752472-1164881352*x*d+
4408536996*x^2*d^2+11774185985*x^3*d^3+5262196786*x^4*d^4+1046030561/2*x^5*
d^5-10564451/2*x^6*d^6)+x^4*(-1925782272+21995375398*x*d+123415803356*x^2*
d^2+302465300831/2*x^3*d^3+34140803907/2*x^4*d^4-15535653409*x^5*d^5-\
2277687768*x^6*d^6)+x^5*(71273525520+691398212366*x*d+901772633569*x^2*d^2+
2281275427069*x^3*d^3+2944352819911/2*x^4*d^4+836872370039/4*x^5*d^5+
9066399237/4*x^6*d^6)+x^6*(2365174430376+9596715855542*x*d+29459572469704*x^
2*d^2+92502197003786*x^3*d^3+65712473180525*x^4*d^4+13829360193674*x^5*d^5
+3231449477251/4*x^6*d^6)+x^7*(26771079436836+117709870166226*x*d+
821686455179082*x^2*d^2+1803972139232179*x^3*d^3+1083654460691481*x^4*d^4+
858903621851785/4*x^5*d^5+50096565802957/4*x^6*d^6)+x^8*(179341727601960+
2144653944040630*x*d+13123246960284302*x^2*d^2+41138357917778169/2*x^3*d^3+
20605819587976401/2*x^4*d^4+3677396642905423/2*x^5*d^5+402688260229369/4*x^6
*d^6)+x^9*(2579190935961288+43587063726809764*x*d+157045086382352387*x^2*d^
2+172175668477370223*x^3*d^3+138636285385875407/2*x^4*d^4+10707836398626232*
x^5*d^5+529435530567584*x^6*d^6)+x^10*(41501953525903392+558336731465626084*
x*d+1407267553543222268*x^2*d^2+1153046693323226808*x^3*d^3+
372331468563656085*x^4*d^4+48654019090240214*x^5*d^5+2114661191282167*x^6*d
^6)+x^11*(364526077273381884+4158060401095928464*x*d+8646807662899324262*x^2*
d^2+5914675753405705400*x^3*d^3+1631934058875116005*x^4*d^4+
187371894330537204*x^5*d^5+7366806367019734*x^6*d^6)+x^12*(
1759850321214603648+18265471270535733520*x*d+34201910114871110912*x^2*d^2+
21265221434709398152*x^3*d^3+5437363546219595036*x^4*d^4+594029113431041060*
x^5*d^5+22881659624561644*x^6*d^6)+x^13*(4648382639403200688+
45699084277107816096*x*d+81049061578449009384*x^2*d^2+48858488665016574368*x
^3*d^3+12515362110098721444*x^4*d^4+1412152747420021048*x^5*d^5+
57196947123984972*x^6*d^6)+x^14*(5459369397960020544+55837825300341621824*x*
d+105671876924055409696*x^2*d^2+71551727420848766624*x^3*d^3+
21094786205096577808*x^4*d^4+2695663190297032192*x^5*d^5+118791751565613264*
x^6*d^6)+x^15*(1023333653580043776+47171127937488813824*x*d+
157258351906685700352*x^2*d^2+145765192195300531840*x^3*d^3+
49876215785510342176*x^4*d^4+6647374188802036864*x^5*d^5+287310278455067312*
x^6*d^6)+x^16*(11960091747366236160+250326608568269289472*x*d+
677587171115580981248*x^2*d^2+538246374825683603456*x^3*d^3+
161380433451548754048*x^4*d^4+19149099315354950144*x^5*d^5+
746433247985092544*x^6*d^6)+x^17*(42246252365448668160+657220532737851248640*
x*d+1531751689216283911680*x^2*d^2+1090829514212206064640*x^3*d^3+
299280728709430851840*x^4*d^4+32932767387222323200*x^5*d^5+
1202281367574179840*x^6*d^6)+x^18*(6239106101942784000+320638742839606579200*
x*d+873857213570556364800*x^2*d^2+645649080101933721600*x^3*d^3+
177008238160627276800*x^4*d^4+19165088507111475200*x^5*d^5+
683600826675660800*x^6*d^6)+x^19*(-60440251454613504000-476055211197689856000
*x*d-733497382597635072000*x^2*d^2-386038662982742016000*x^3*d^3-\
83361486778142976000*x^4*d^4-7524999543181824000*x^5*d^5-232189492987008000*
x^6*d^6)+x^20*(1578562930483200000+12628503443865600000*x*d+
19732036631040000000*x^2*d^2+10523752869888000000*x^3*d^3+
2302070940288000000*x^4*d^4+210475057397760000*x^5*d^5+6577345543680000*x^6*
d^6);
printlevel = 3;
facFirstWeyl(dopp);
*/



//Hard Example by Viktor:
/*
  ring r = 0,(x,d), (dp);
def R = nc_algebra(1,1);
setring R;
LIB "/Users/albertheinle/Studium/forschung/ncfactor/versionen/ncfactor.lib";
poly t = x; poly D =d;
poly p = 2*t^2*D^8-6*t*D^8+2*t^2*D^7+8*t*D^7+12*D^7-2*t^4*D^6+6*t^3*D^6+12*t*D^6-20*D^6
-2*t^4*D^5-8*t^3*D^5-4*t^2*D^5+12*t*D^5-28*D^5-12*t^3*D^4-4*t^2*D^4-4*t*D^4-24*D^4+4*t^4*D^3
-12*t^3*D^3+2*t^2*D^3-18*t*D^3+16*D^3+6*t^4*D^2-2*t^3*D^2+2*t^2*D^2-2*t*D^2+44*D^2+2*t^4*D
+12*t^3*D+2*t*D+4*t^3-8;
list lsng = facFirstWeyl(p);
print(lsng);
*/