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

///////////////////////////////////////////////////////////
version = "$Id$";
category="Noncommutative";
info="
LIBRARY: ncfactor.lib  Tools for factorization in some noncommutative algebras
AUTHORS: Albert Heinle,     albert.heinle@rwth-aachen.de
@*       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]);   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
";

LIB "general.lib";
LIB "nctools.lib";
LIB "involut.lib";
LIB "freegb.lib"; // for isVar

proc tst_ncfactor()
{
  example facFirstWeyl;
  example facFirstShift;
  example facSubWeyl;
  example testNCfac;
}

/////////////////////////////////////////////////////
//==================================================*
//deletes double-entries in a list of factorization
//without evaluating the product.
static proc delete_dublicates_noteval(list l)
{//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)
{//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


//==================================================*
//given a list of factors g and a desired size nof, the following
//procedure combines the factors, such that we recieve a
//list of the length nof.
static proc combinekfinlf(list g, int nof, intvec limits) //nof stands for "number of factors"
{//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
    if (limitcheck(g,limits))
    {
      result = result + list(g);
    }
    return(result);
  }//There are no factors to combine
  if (nof == 1)
  {//User wants to get just one factor
    if (limitcheck(list(product(g)),limits))
    {
      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,limits);
      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
        if (limitcheck(temp,limits))
        {
          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/nof;i++)
  {//getting the other possible results
    result = result + combinekfinlf(list(product(list(g[1..i])))+list(g[(i+1)..nofgl]),nof,limits);
  }//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)
{//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);
      i--;
    }
  }//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)
{//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
  return(result);
}//proc merge_cf


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

static proc mergence(list l1, list l2, intvec limits)
{//Procedure mergence
  list g;
  list f;
  int l; int k;
  list F;
  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 result;
  for (l = size(f); l>=1; l--)
  {//all possibilities to combine the factors of f
    F = combinekfinlf(f,l,limits);
    for (k = 1; k<= size(F);k++)
    {//for all possibilities of combinations of the factors of f
      result = result + merge_cf(F[k],g,limits);
      result = result + merge_icf(F[k],g,limits);
    }//for all possibilities of combinations of the factors of f
  }//all possibilities to combine the factors of f
  return(result);
}//Procedure mergence


//==================================================
//Checks, whether a list of factors doesn't exceed the given limits
static proc limitcheck(list g, intvec limits)
{//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;
  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,"==> 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[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,"==> 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;
  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[\theta], 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;
  intvec ivm11 = intvec(-1,1);
  dbprint(p,"==> Calculate one homogeneous factorization using homogfacFirstWeyl");
  //Compute again a homogeneous factorization
  one_hom_fac = homogfacFirstWeyl(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 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,"==> 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]");
  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,"==> 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,"==> 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,"==> 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],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,"==> 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 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)
{//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)
{//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 in 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;
  dbprint(p,"==> 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,"==> Successful");
  dbprint(p,"==> 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
  dbprint(p,"==> Successful");
  list result = list();
  int i;int j; int k; int l; //counter
  if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why
  {
    dbprint(p,"==> 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,"==> Successful");
    dbprint(p, "==> factorization of the polynomial with the routine sfacwa");
    list resulttemp = sfacwa(transf(h));
    dbprint(p,"==> Done");
    setring(r);
    map transfback = NTR, var(2),var(1);
    result = transfback(resulttemp);
  }
  else
  {
    dbprint(p, "==> factorization of the polynomial with the routine sfacwa");
    result = sfacwa(h);
    dbprint(p,"==> Done");
  }
  dbprint(p,"==> recursively check factors for irreducibility");
  list recursivetemp;
  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
          for(k=1; k<=size(recursivetemp);k++)
          {//insert factorized factors
            if(size(recursivetemp[k])>2)
            {//nontrivial
              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
        }//we have a nontrivial factorization
      }//Factorize every factor
    }//Nontrivial factorization
  }//recursively factorize factors
  dbprint(p,"==> Done");
  if (size(result)==0)
  {//only the trivial factorization could be found
    result = list(list(1,h));
  }//only the trivial factorization could be found
  //now, refine the possible redundant list
  return( delete_dublicates_noteval(result) );
}//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);
}

//This is the main program
static proc sfacwa(poly h)
"USAGE: sfacwa(h); h is a polynomial in the first Weyl algebra
RETURN: list
PURPOSE: Computes a factorization of a polynomial h in the first Weyl algebra
THEORY: @code{sfacwa} returns a list with some factorizations of the given
        polynomial. The possibilities of the factorization of the highest
        homogeneous part and those of the lowest will be merged. If during this
        procedure a factorization of the polynomial occurs, it will be added to
        the output list. For a more detailed description visit
        @url{http://www.math.rwth-aachen.de/\~Albert.Heinle}
SEE ALSO: homogfacFirstWeyl_all, homogfacFirstWeyl
"{//proc sfacwa
  int p=printlevel-voice+2; // for dbprint
  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);
  if(homogwithorder(h,ivm11))
  {
    dbprint(p,"==> Given polynomial is -1,1 homogeneous. Starting homog. fac. and ret. its result");
    return(homogfacFirstWeyl_all(h));
  }
  def r = basering;
  map invo = basering,-var(1),var(2);
  int i; int j; int k;
  dbprint(p,"==> Computing the degree-limits of the factorization");
  intvec limits = deg(h,iv11) ,deg(h,iv10),deg(h,iv01);
  def prod;
  //end finding the limits
  dbprint(p,"==> Computing the maximal and the minimal homogeneous part of the give polynomial");
  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 result;
  list temp;
  list homogtemp;
  def invhath;
  def invideal;
  list M; list hatM;
  dbprint(p,"==> Factorize the maximal and the minimal homog. part of the given polynomial");
  list f1 = homogfacFirstWeyl_all(maxh);
  list f2 = homogfacFirstWeyl_all(minh);
  int is_equal;
  poly hath;
  dbprint(p,"==> Merging the factorizations");
  for (i = 1; i<=size(f1);i++)
  {//Merge all combinations
    for (j = 1; j<=size(f2); j++)
    {
      M = M + mergence(f1[i],f2[j],limits);
    }
  }//Merge all combinations
  dbprint(p,"==> Filtering invalid combinations");
  for (i = 1 ; i<= size(M); i++)
  {//filter valid combinations
    if (product(M[i]) == h)
    {//We have one factorization
      result = result + list(M[i]);
      M = delete(M,i);
      continue;
    }//We have one factorization
    else
    {
      if (deg(h,ivm11)<=deg(h-product(M[i]),ivm11))
      {
        M = delete(M,i);
        continue;
      }
      if (deg(h,iv1m1)<=deg(h-product(M[i]),iv1m1))
      {
        M = delete(M,i);
        continue;
      }
    }
  }//filter valid combinations
  dbprint(p,"==> Deleting doublicates in the resulting list");
  M = delete_dublicates_eval(M);
  dbprint(p,"==> Iterating over all possible Combinations");
  while(size(M)>0)
  {//work on the elements of M
    hatM = list();
    for(i = 1; i<=size(M); i++)
    {//iterate over all elements of M
      hath = h-product(M[i]);
      temp = list();
      //First check for common inhomogeneous factors between hath and h
      invhath = involution(hath,invo);
      invideal = std(involution(ideal(M[i][1]),invo));
      if (involution(NF(invhath,invideal) ,invo)==0)
      {//hath and h have a common factor on the left
        j = 1;
        f1 = M[i];
        if (j+1<=size(f1))
        {//Checking for more than one common factor
          invideal = std(involution(ideal(product(f1[1..(j+1)])),invo));
          while(involution(NF(invhath,invideal),invo)==0)
          {
            if (j+1<size(f1))
            {
              j++;
            }
            else
            {
              break;
            }
          }
        }//Checking for more than one common factor
        invideal = lift(involution(product(f1[1..j]),invo),invhath);
        f2 = list(f1[1..j])+list(involution(invideal[1,1],invo));
        temp = temp + merge_cf(f2,f1,limits);
      }//hath and h have a common factor on the left
      if (reduce(hath, std(ideal(M[i][size(M[i])])))==0)
      {//hath and h have a common factor on the right
        j = size(M[i]);
        f1 = M[i];
        if (j-1>0)
        {//Checking for more than one factor
          while(reduce(hath,std(ideal(product(f1[(j-1)..size(f1)]))))==0)
          {
            if (j-1>1)
            {
              j--;
            }
            else
            {
              break;
            }
          }
        }//Checking for more than one factor
        f2 = list(lift(product(f1[j..size(f1)]),hath)[1,1])+list(f1[j..size(f1)]);
        temp = temp + merge_cf(f2,M[i],limits);
      }//hath and h have a common factor on the right
      //and now the homogeneous
      maxh = jet(hath,deg(hath,ivm11),ivm11)-jet(hath,deg(hath,ivm11)-1,ivm11);
      minh = jet(hath,deg(hath,iv1m1),iv1m1)-jet(hath,deg(hath,iv1m1)-1,iv1m1);
      f1 = homogfacFirstWeyl_all(maxh);
      f2 = homogfacFirstWeyl_all(minh);
      for (j = 1; j<=size(f1);j++)
      {
        for (k=1; k<=size(f2);k++)
        {
          homogtemp = mergence(f1[j],f2[k],limits);
        }
      }
      for (j = 1; j<= size(homogtemp); j++)
      {
        temp = temp + mergence(homogtemp[j],M[i],limits);
      }
      for (j = 1; j<=size(temp); j++)
      {//filtering invalid entries in temp
        if(product(temp[j])==h)
        {//This is already a result
          result = result + list(temp[j]);
          temp = delete(temp,j);
          continue;
        }//This is already a result
        if (deg(hath,ivm11)<=deg(hath-product(temp[j]),ivm11))
        {
          temp = delete(temp,j);
          continue;
        }
      }//filtering invalid entries in temp
      hatM = hatM + temp;
    }//iterate over all elements of M
    M = hatM;
    for (i = 1; i<=size(M);i++)
    {//checking for complete factorizations
      if (h == product(M[i]))
      {
        result = result + list(M[i]);
        M = delete(M,i);
        continue;
      }
    }//checking for complete factorizations
    M = delete_dublicates_eval(M);
  }//work on the elements of M
  //In the case, that there is none, write a constant factor before the factor of interest.
  dbprint(p,"==> Done");
  for (i = 1 ; i<=size(result);i++)
  {//add a constant factor
    if (deg(result[i][1],iv11)!=0)
    {
      result[i] = insert(result[i],1);
    }
  }//add a constant factor
  dbprint(p,"==> Deleting doublicates in the output list");
  result = delete_dublicates_noteval(result);
  return(result);
}//proc sfacwa


//==================================================
/*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;
  dbprint(p,"==> Checking the input");
  if (size(l)==0)
  {//The empty list is given
    dbprint(p,"==> Given list was empty");
    return(list());
  }//The empty list is given
  if (size(#)>2)
  {//We want max. two optional arguments
    dbprint(p,"==> More than two optional arguments");
    return(list());
  }//We want max. two optional arguments
  dbprint(p,"==> Done");
  list result;
  int i; int j;
  if (size(#)==0)
  {//No optional argument is given
    dbprint(p,"==> 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,"==> 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,"==> 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, X, D)
"USAGE:  facSubWeyl(h,x,y); 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;
  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,"==> 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,"==> 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;
  int shiftcounter;
  //Compute again a homogeneous factorization
  dbprint(p,"==> 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,"==> 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,"==> 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 in 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;
  dbprint(p,"==> 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,"==> Done");
  list result = list();
  int i;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,"==> 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,"==> Factorize the given polynomial with the subroutine facshift");
  list resulttemp = facshift(transf(h));
  dbprint(p,"==> Successful");
  setring(r);
  result = transfback(resulttemp);
  dbprint(p,"==> Recursively check the found factors for reducibility.");
  list recursivetemp;
  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 = facFirstShift(result[i][j]);
        if(size(recursivetemp)>1)
        {//we have a nontrivial factorization
          for(k=1; k<=size(recursivetemp);k++)
          {//insert factorized factors
            if(size(recursivetemp[k])>2)
            {//nontrivial
              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
        }//we have a nontrivial factorization
      }//Factorize every factor
    }//Nontrivial factorization
  }//recursively factorize factors
  //now, refine the possible redundant list
  dbprint(p,"==> Done");
  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 facshift(poly h)
"USAGE: facshift(h); h is a polynomial in the first Shift algebra
RETURN: list
PURPOSE: Computes a factorization of a polynomial h in the first Shift algebra
THEORY: @code{facshift} returns a list with some factorizations of the given
        polynomial. The possibilities of the factorization of the highest
        homogeneous part and those of the lowest will be merged. If during this
        procedure a factorization of the polynomial occurs, it will be added to
        the output list. For a more detailled description visit
        @url{http://www.math.rwth-aachen.de/\~Albert.Heinle}
SEE ALSO: homogfacFirstShift_all, homogfacFirstShift
"{//proc facshift
  int p=printlevel-voice+2; // for dbprint
  intvec iv01 = intvec(0,1);
  intvec iv11 = intvec(1,1);
  intvec iv10 = intvec(1,0);
  intvec iv0m1 = intvec(0,-1);
  if(homogwithorder(h,iv01))
  {
    dbprint(p,"==> Given polynomial was 1,0-homogeneous. Calling subroutine for homog. polys.");
    return(homogfacFirstShift_all(h));
  }
  def r = basering;
  map invo = basering,-var(1),-var(2);
  int i; int j; int k;
  intvec limits = deg(h,iv11) ,deg(h,iv10),deg(h,iv01);
  def prod;
  //end finding the limits
  poly maxh = jet(h,deg(h,iv01),iv01)-jet(h,deg(h,iv01)-1,iv01);
  poly minh = jet(h,deg(h,iv0m1),iv0m1)-jet(h,deg(h,iv0m1)-1,iv0m1);
  list result;
  list temp;
  list homogtemp;
  list M; list hatM;
  def invhath;
  def invideal;
  dbprint(p,"==> Factorizing the highest and the lowest homogeneous part");
  list f1 = homogfacFirstShift_all(maxh);
  list f2 = homogfacFirstShift_all(minh);
  int is_equal;
  poly hath;
  dbprint(p,"==> Merging their factors");
  for (i = 1; i<=size(f1);i++)
  {//Merge all combinations
    for (j = 1; j<=size(f2); j++)
    {
      M = M + mergence(f1[i],f2[j],limits);
    }
  }//Merge all combinations
  dbprint(p,"==> Filtering invalid combinations");
  for (i = 1 ; i<= size(M); i++)
  {//filter valid combinations
    if (product(M[i]) == h)
    {//We have one factorization
      result = result + list(M[i]);
      M = delete(M,i);
      continue;
    }//We have one factorization
    else
    {
      if (deg(h,iv01)<=deg(h-product(M[i]),iv01))
      {
        M = delete(M,i);
        continue;
      }
      if (deg(h,iv0m1)<=deg(h-product(M[i]),iv0m1))
      {
        M = delete(M,i);
        continue;
      }
    }
  }//filter valid combinations
  dbprint(p,"==> Iterating over all elements in M");
  M = delete_dublicates_eval(M);
  while(size(M)>0)
  {//work on the elements of M
    hatM = list();
    for(i = 1; i<=size(M); i++)
    {//iterate over all elements of M
      hath = h-product(M[i]);
      temp = list();
      //First check for common inhomogeneous factors between hath and h
      invhath = involution(hath,invo);
      invideal = std(involution(ideal(M[i][1]),invo));
      if (involution(NF(invhath, invideal ),invo)==0)
      {//hath and h have a common factor on the left
        j = 1;
        f1 = M[i];
        if (j+1<=size(f1))
        {//Checking for more than one common factor
          invideal = std(involution(ideal(product(f1[1..(j+1)])),invo));
          while(involution(NF(invhath,invideal),invo)==0)
          {
            if (j+1<size(f1))
            {
              j++;
            }
            else
            {
              break;
            }
          }
        }//Checking for more than one common factor
        if (deg(product(f1[1..j]),iv11)!=0)
        {
          invideal = lift(involution(product(f1[1..j]),invo),invhath);
          f2 = list(f1[1..j])+list(involution(invideal[1,1],invo));
        }
        else
        {
          invideal = lift(product(f1[1..j]),involution(hath,invo));
          f2 = list(f1[1..j])+list(involution(invideal[1,1],invo));
        }
        temp = temp + merge_cf(f2,f1,limits);
      }//hath and h have a common factor on the left
      if (reduce(hath, std(ideal(M[i][size(M[i])])))==0)
      {//hath and h have a common factor on the right
        j = size(M[i]);
        f1 = M[i];
        if (j-1>0)
        {//Checking for more than one factor
          while(reduce(hath,std(ideal(product(f1[(j-1)..size(f1)]))))==0)
          {
            if (j-1>1)
            {
              j--;
            }
            else
            {
              break;
            }
          }
        }//Checking for more than one factor
        f2 = list(lift(product(f1[j..size(f1)]),hath)[1,1])+list(f1[j..size(f1)]);
        temp = temp + merge_cf(f2,M[i],limits);
      }//hath and h have a common factor on the right
      //and now the homogeneous
      maxh = jet(hath,deg(hath,iv01),iv01)-jet(hath,deg(hath,iv01)-1,iv01);
      minh = jet(hath,deg(hath,iv0m1),iv0m1)-jet(hath,deg(hath,iv0m1)-1,iv0m1);
      f1 = homogfacFirstShift_all(maxh);
      f2 = homogfacFirstShift_all(minh);
      for (j = 1; j<=size(f1);j++)
      {
        for (k=1; k<=size(f2);k++)
        {
          homogtemp = mergence(f1[j],f2[k],limits);
        }
      }
      for (j = 1; j<= size(homogtemp); j++)
      {
        temp = temp + mergence(homogtemp[j],M[i],limits);
      }
      for (j = 1; j<=size(temp); j++)
      {//filtering invalid entries in temp
        if(product(temp[j])==h)
        {//This is already a result
          result = result + list(temp[j]);
          temp = delete(temp,j);
          continue;
        }//This is already a result
        if (deg(hath,iv01)<=deg(hath-product(temp[j]),iv01))
        {
          temp = delete(temp,j);
          continue;
        }
      }//filtering invalid entries in temp
      hatM = hatM + temp;
    }//iterate over all elements of M
    M = hatM;
    for (i = 1; i<=size(M);i++)
    {//checking for complete factorizations
      if (h == product(M[i]))
      {
        result = result + list(M[i]);
        M = delete(M,i);
        continue;
      }
    }//checking for complete factorizations
    M = delete_dublicates_eval(M);
  }//work on the elements of M
  //In the case, that there is none, write a constant factor before the factor of interest.
  for (i = 1 ; i<=size(result);i++)
  {//add a constant factor
    if (deg(result[i][1],iv11)!=0)
    {
      result[i] = insert(result[i],1);
    }
  }//add a constant factor
  dbprint(p,"==> Deleting double entries in the resulting list.");
  result = delete_dublicates_noteval(result);
  return(result);
}//proc facshift

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);
}

/*
  Example polynomials where one can find factorizations: K<x,y |yx=xy+1>
  (x^2+y)*(x^2+y);
  (x^2+x)*(x^2+y);
  (x^3+x+1)*(x^4+y*x+2);
  (x^2*y+y)*(y+x*y);
  y^3+x*y^3+2*y^2+2*(x+1)*y^2+y+(x+2)*y; //Example 5 Grigoriev-Schwarz.
  (y+1)*(y+1)*(y+x*y); //Landau Example projected to the first dimension.
*/


/* very hard things from Martin Lee:
// ex1, ex2
ring s = 0,(z,x),Ws(-1,1);
def S = nc_algebra(1,1); setring S;
poly a = 10z5x4+26z4x5+47z5x2-97z4x3; //Abgebrochen nach einer Stunde; yes, it takes long
def l= facFirstWeyl (a); l;
kill l;
poly b = -5328z8x5-5328z7x6+720z9x2+720z8x3-16976z7x4-38880z6x5-5184z7x3-5184z6x4-3774z5x5+2080z8x+5760z7x2-6144z6x3-59616z5x4+3108z3x6-4098z6x2-25704z5x3-21186z4x4+8640z6x-17916z4x3+22680z2x5+2040z5x-4848z4x2-9792z3x3+3024z2x4-10704z3x2-3519z2x3+34776zx4+12096zx3+2898x4-5040z2x+8064x3+6048x2; //Abgebrochen nach 1.5 Stunden; seems to be very complicated
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!

*/