////////////////////////////////////////////////////////////////////
version="$Id$";
category="Symbolic-numerical solving";
info="
LIBRARY: ffsolve.lib        multivariate equation solving over finite fields
AUTHOR: Gergo Gyula Borus, borisz@borisz.net
KEYWORDS: multivariate equations; finite field

PROCEDURES:
ffsolve();         finite field solving using heuristically chosen method
PEsolve();         solve system of multivariate equations over finite field
simplesolver();    solver using modified exhausting search
GBsolve();         multivariate solver using Groebner-basis
XLsolve();         multivariate polynomial solver using linearization
ZZsolve();         solve system of multivariate equations over finite field
";

LIB "presolve.lib";
LIB "general.lib";
LIB "ring.lib";
LIB "standard.lib";
LIB "matrix.lib";

////////////////////////////////////////////////////////////////////
proc ffsolve(ideal equations, list #)
"USAGE:         ffsolve(I[, L]); I ideal, L list of strings
RETURN:         list L, the common roots of I as ideal
ASSUME:         basering is a finite field of type (p^n,a)
"
{
  list solutions, lSolvers, tempsols;
  int i,j, k,n, R, found;
  ideal factors, linfacs;
  poly lp;
  // check assumptions
  if(npars(basering)>1){
    ERROR("Basering must have at most one parameter");
  }
  if(char(basering)==0){
    ERROR("Basering must have finite characteristic");
  }

  if(size(#)){
    if(size(#)==1 and typeof(#[1])=="list"){
      lSolvers = #[1];
    }else{
      lSolvers = #;
    }
  }else{
    if(deg(equations) == 2){
      lSolvers = "XLsolve", "PEsolve", "simplesolver", "GBsolve", "ZZsolve";
    }else{
      lSolvers = "PEsolve", "simplesolver", "GBsolve", "ZZsolve", "XLsolve";
    }
    if(deg(equations) == 1){
      lSolvers = "GBsolve";
    }
  }
  n = size(lSolvers);
  R = random(1, n*(3*n+1) div 2);
  for(i=1;i<n+1;i++){
    if(R<=(2*n+1-i)){
      string solver = lSolvers[i];
    }else{
      R=R-(2*n+1-i);
    }
  }

  if(nvars(basering)==1){
    return(facstd(equations));
  }

  // search for the first univariate polynomial
  found = 0;
  for(i=1; i<=ncols(equations); i++){
    if(univariate(equations[i])>0){
      factors=factorize(equations[i],1);
      for(j=1; j<=ncols(factors); j++){
        if(deg(factors[j])==1){
          linfacs[size(linfacs)+1] = factors[j];
        }
      }
      if(deg(linfacs[1])>0){
        found=1;
        break;
      }
    }
  }
  //   if there is, collect its the linear factors
  if(found){
    // substitute the root and call recursively
    ideal neweqs, invmapideal, ti;
    map invmap;
    for(k=1; k<=ncols(linfacs); k++){
      lp = linfacs[k];
      neweqs = reduce(equations, lp);

      intvec varexp = leadexp(lp);
      def original_ring = basering;
      def newRing = clonering(nvars(original_ring)-1);
      setring newRing;
      ideal mappingIdeal;
      j=1;
      for(i=1; i<=size(varexp); i++){
        if(varexp[i]){
          mappingIdeal[i] = 0;
        }else{
          mappingIdeal[i] = var(j);
          j++;
        }
      }
      map recmap = original_ring, mappingIdeal;
      list tsols = ffsolve(recmap(neweqs), lSolvers);
      if(size(tsols)==0){
        tsols = list(ideal(1));
      }
      setring original_ring;
      j=1;
      for(i=1;i<=size(varexp);i++){
        if(varexp[i]==0){
          invmapideal[j] = var(i);
          j++;
        }
      }
      invmap = newRing, invmapideal;
      tempsols = invmap(tsols);

      // combine the solutions
      for(j=1; j<=size(tempsols); j++){
        ti =  std(tempsols[j]+lp);
        if(deg(ti)>0){
          solutions = insert(solutions,ti);
        }
      }
    }
  }else{
    execute("solutions="+solver+"(equations);") ;
  }
  return(solutions);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = (2,a),x(1..3),lp;
  minpoly=a2+a+1;
  ideal I;
  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
  ffsolve(I);
}
////////////////////////////////////////////////////////////////////
proc PEsolve(ideal L, list #)
"USAGE:         PEsolve(I[, i]); I ideal, i optional integer
                solve I (system of multivariate equations) over a
                finite field using an equvalence property when i is
                not given or set to 2, otherwise if i is set to 0
                then check whether common roots exists
RETURN:         list if optional parameter is not given or set to 2,
                integer if optional is set to 0
ASSUME:         basering is a finite field of type (p^n,a)
NOTE:           When the optional parameter is set to 0, speoff only
                checks if I has common roots, then return 1, otherwise
                return 0.

"
{
  int mode, i,j;
  list results, rs, start;
  poly g;
  // check assumptions
  if(npars(basering)>1){
    ERROR("Basering must have at most one parameter");
  }
  if(char(basering)==0){
    ERROR("Basering must have finite characteristic");
  }

  if( size(#) > 0 ){
    mode = #[1];
  }else{
    mode = 2;
  }
  L = simplify(L,15);
  g = productOfEqs( L );

  if(g == 0){
    if(mode==0){
      return(0);
    }
    return( list() );
  }
  if(g == 1){
    list vectors = every_vector();
    for(j=1; j<=size(vectors); j++){
      ideal res;
      for(i=1; i<=nvars(basering); i++){
        res[i] = var(i)-vectors[j][i];
      }
      results[size(results)+1] = std(res);
    }
    return( results );
  }

  if( mode == 0 ){
    return( 1 );
  }else{
    for(i=1; i<=nvars(basering); i++){
      start[i] = 0:order_of_extension();
    }

    if( mode == 1){
      results[size(results)+1] = melyseg(g, start);
    }else{
      while(1){
        start = melyseg(g, start);
        if( size(start) > 0 ){
          ideal res;
          for(i=1; i<=nvars(basering); i++){
            res[i] = var(i)-vec2elm(start[i]);
          }
          results[size(results)+1] = std(res);
          start = increment(start);
        }else{
          break;
        }
      }
    }
  }
  return(results);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = (2,a),x(1..3),lp;
  minpoly=a2+a+1;
  ideal I;
  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
  PEsolve(I);
}
////////////////////////////////////////////////////////////////////
proc simplesolver(ideal E)
"USAGE:         simplesolver(I); I ideal
                solve I (system of multivariate equations) over a
                finite field by exhausting search
RETURN:         list L, the common roots of I as ideal
ASSUME:         basering is a finite field of type (p^n,a)
"
{
  int i,j,k,t, correct;
  list solutions = list(std(ideal()));
  list partial_solutions;
  ideal partial_system, curr_sol, curr_sys, factors;
  poly univar_poly;
  E = E+defaultIdeal();
  // check assumptions
  if(npars(basering)>1){
    ERROR("Basering must have at most one parameter");
  }
  if(char(basering)==0){
    ERROR("Basering must have finite characteristic");
  }
  for(k=1; k<=nvars(basering); k++){
    partial_solutions = list();
    for(i=1; i<=size(solutions); i++){
      partial_system = reduce(E, solutions[i]);
      for(j=1; j<=ncols(partial_system); j++){
        if(univariate(partial_system[j])>0){
          univar_poly = partial_system[j];
          break;
        }
      }
      factors = factorize(univar_poly,1);
      for(j=1; j<=ncols(factors); j++){
        if(deg(factors[j])==1){
          curr_sol = std(solutions[i]+ideal(factors[j]));
          curr_sys = reduce(E, curr_sol);
          correct = 1;
          for(t=1; t<=ncols(curr_sys); t++){
            if(deg(curr_sys[t])==0){
              correct = 0;
              break;
            }
          }
          if(correct){
            partial_solutions = insert(partial_solutions, curr_sol);
          }
        }
      }
    }
    solutions = partial_solutions;
  }
  return(solutions);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = (2,a),x(1..3),lp;
  minpoly=a2+a+1;
  ideal I;
  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
  simplesolver(I);
}
////////////////////////////////////////////////////////////////////
proc GBsolve(ideal equation_system)
"USAGE:         GBsolve(I); I ideal
                solve I (system of multivariate equations) over an
                extension of Z/p by Groebner basis methods
RETURN:         list L, the common roots of I as ideal
ASSUME:         basering is a finite field of type (p^n,a)
"
{
  int i,j, prop, newelement, number_new_vars;
  ideal ls;
  list results, slvbl, linsol, ctrl, new_sols, varinfo;
  ideal I, linear_solution, unsolved_part, univar_part, multivar_part, unsolved_vars;
  intvec unsolved_var_nums;
  string new_vars;
  // check assumptions
  if(npars(basering)>1){
    ERROR("Basering must have at most one parameter");
  }
  if(char(basering)==0){
    ERROR("Basering must have finite characteristic");
  }

  def original_ring = basering;
  if(npars(basering)==1){
    int prime_coeff_field=0;
    string minpolystr = "minpoly="+
    get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
  }else{
    int prime_coeff_field=1;
  }

  option(redSB);

  equation_system = simplify(equation_system,15);

  ideal standard_basis = std(equation_system);
  list basis_factors = facstd(standard_basis);
  if( basis_factors[1][1] == 1){
    return(results)
  };

  for(i=1; i<= size(basis_factors); i++){
    prop = 0;
    for(j=1; j<=size(basis_factors[i]); j++){
      if( univariate(basis_factors[i][j])>0 and deg(basis_factors[i][j])>1){
        prop =1;
        break;
      }
    }
    if(prop == 0){
      ls = solvelinearpart( basis_factors[i] );
      if(ncols(ls) == nvars(basering) ){
        ctrl, newelement = add_if_new(ctrl, ls);
        if(newelement){
          results = insert(results, ls);
        }
      }else{
        slvbl = insert(slvbl, list(basis_factors[i],ls) );
      }
    }
  }
  if(size(slvbl)<>0){
    for(int E = 1; E<= size(slvbl); E++){
      I = slvbl[E][1];
      linear_solution = slvbl[E][2];
      attrib(I,"isSB",1);
      unsolved_part = reduce(I,linear_solution);
      univar_part = ideal();
      multivar_part = ideal();
      for(i=1; i<=ncols(I); i++){
        if(univariate(I[i])>0){
          univar_part = univar_part+I[i];
        }else{
          multivar_part = multivar_part+I[i];
        }
      }
      varinfo = findvars(univar_part,1);
      unsolved_vars = varinfo[3];
      unsolved_var_nums = varinfo[4];
      number_new_vars = ncols(unsolved_vars);

      new_vars = "@y(1.."+string(number_new_vars)+")";
      def R_new = changevar(new_vars, original_ring);
      setring R_new;
      if( !prime_coeff_field ){
        execute(minpolystr);
      }

      ideal mapping_ideal;
      for(i=1; i<=size(unsolved_var_nums); i++){
        mapping_ideal[unsolved_var_nums[i]] = var(i);
      }

      map F = original_ring, mapping_ideal;
      ideal I_new = F( multivar_part );

      list sol_new;
      int unsolvable = 0;
      sol_new = simplesolver(I_new);
      if( size(sol_new) == 0){
        unsolvable = 1;
      }

      setring original_ring;

      if(unsolvable){
        list sol_old = list();
      }else{
        map G = R_new, unsolved_vars;
        new_sols = G(sol_new);
        for(i=1; i<=size(new_sols); i++){
          ideal sol = new_sols[i]+linear_solution;
          sol = std(sol);
          ctrl, newelement = add_if_new(ctrl, sol);
          if(newelement){
            results = insert(results, sol);
          }
          kill sol;
        }
      }
      kill G;
      kill R_new;
    }
  }
  return( results  );
}
example
{
  "EXAMPLE:";echo=2;
  ring R = (2,a),x(1..3),lp;
  minpoly=a2+a+1;
  ideal I;
  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
  GBsolve(I);
}
////////////////////////////////////////////////////////////////////
proc XLsolve(ideal I, list #)
"USAGE:         XLsolve(I[, d]); I ideal, d optional integer
                solve I (system of multivariate polynomials) with a
                variant of the linearization technique, multiplying
                the polynomials with monomials of degree at most d
                (default is 2)
RETURN:         list L of the common roots of I as ideals
ASSUME:         basering is a finite field of type (p^n,a)"
{
  int i,j,k, D;
  int SD = deg(I);
  list solutions;
  if(size(#)){
    if(typeof(#[1])=="int"){ D = #[1]; }
  }else{
    D = 2;
  }
  list lMonomialsForMultiplying = monomialsOfDegreeAtMost(D+SD);

  int m = ncols(I);
  list extended_system;
  list mm;
  for(k=1; k<=size(lMonomialsForMultiplying)-SD; k++){
    mm = lMonomialsForMultiplying[k];
    for(i=1; i<=m; i++){
      for(j=1; j<=size(mm); j++){
        extended_system[size(extended_system)+1] = reduce(I[i]*mm[j], defaultIdeal());
      }
    }
  }
  ideal new_system = I;
  for(i=1; i<=size(extended_system); i++){
    new_system[m+i] = extended_system[i];
  }
  ideal reduced_system = linearReduce( new_system, lMonomialsForMultiplying);

  solutions = simplesolver(reduced_system);

  return(solutions);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = (2,a),x(1..3),lp;
  minpoly=a2+a+1;
  ideal I;
  I[1]=(a)*x(1)^2+x(2)^2+(a+1);
  I[2]=(a)*x(1)^2+(a)*x(1)*x(3)+(a)*x(2)^2+1;
  I[3]=(a)*x(1)*x(3)+1;
  I[4]=x(1)^2+x(1)*x(3)+(a);
  XLsolve(I, 3);
}

////////////////////////////////////////////////////////////////////
proc ZZsolve(ideal I)
"USAGE:         ZZsolve(I); I ideal
solve I (system of multivariate equations) over a
finite field by mapping the polynomials to a single
univariate polynomial over extension of the basering
RETURN:         list, the common roots of I as ideal
ASSUME:         basering is a finite field of type (p^n,a)
"
{
  int i, j, nv, numeqs,r,l,e;
  def original_ring = basering;
  // check assumptions
  if(npars(basering)>1){
    ERROR("Basering must have at most one parameter");
  }
  if(char(basering)==0){
    ERROR("Basering must have finite characteristic");
  }

  nv = nvars(original_ring);
  numeqs = ncols(I);
  l = numeqs % nv;
  if( l == 0){
    r = numeqs div nv;
  }else{
    r = (numeqs div nv) +1;
  }


  list list_of_equations;
  for(i=1; i<=r; i++){
    list_of_equations[i] = ideal();
  }
  for(i=0; i<numeqs; i++){
    list_of_equations[(i div nv)+1][(i % nv) +1] = I[i+1];
  }

  ring ring_for_matrix = (char(original_ring),@y),(x(1..nv),@X,@c(1..nv)(1..nv)),lp;
  execute("minpoly="+Z_get_minpoly(size(original_ring)^nv, parstr(1))+";");

  ideal IV;
  for(i=1; i<=nv; i++){
    IV[i] = var(i);
  }

  matrix M_C[nv][nv];
  for(i=1;i<=nrows(M_C); i++){
    for(j=1; j<=ncols(M_C); j++){
      M_C[i,j] = @c(i)(j);
    }
  }

  poly X = Z_phi(IV);
  ideal IX_power_poly;
  ideal IX_power_var;
  for(i=1; i<=nv; i++){
    e = (size(original_ring)^(i-1));
    IX_power_poly[i] = X^e;
    IX_power_var[i] = @X^e;
  }
  IX_power_poly = reduce(IX_power_poly, Z_default_ideal(nv, size(original_ring)));

  def M = matrix(IX_power_poly,1,nv)*M_C - matrix(IV,1,nv);

  ideal IC;
  for(i=1; i<=ncols(M); i++){
    for(j=1; j<=ncols(IV); j++){
      IC[(i-1)*ncols(M)+j] = coeffs(M[1,i],IV[j])[2,1];
    }
  }

  ideal IC_solultion = std(Presolve::solvelinearpart(IC));


  matrix M_C_sol[nv][nv];
  for(i=1;i<=nrows(M_C_sol); i++){
    for(j=1; j<=ncols(M_C_sol); j++){
      M_C_sol[i,j] = reduce(@c(i)(j), std(IC_solultion));
    }
  }
  ideal I_subs;
  I_subs = ideal(matrix(IX_power_var,1,nv)*M_C_sol);

  setring original_ring;
  string var_str = varstr(original_ring)+",@X,@y";
  string minpoly_str = "minpoly="+string(minpoly)+";";
  def ring_for_substitution = Ring::changevar(var_str, original_ring);

  setring ring_for_substitution;
  execute(minpoly_str);

  ideal I_subs = imap(ring_for_matrix, I_subs);
  ideal I = imap(original_ring, I);
  list list_of_equations = imap(original_ring, list_of_equations);

  list list_of_F;
  for(i=1; i<=r; i++){
    list_of_F[i] = Z_phi( list_of_equations[i] );
  }

  for(i=1; i<=nv; i++){

    for(j=1; j<=r; j++){
      list_of_F[j] = subst( list_of_F[j], var(i), I_subs[i] );
    }
  }
  int s = size(original_ring);
  if(npars(original_ring)==1){
    for(j=1; j<=r; j++){
      list_of_F[j] = subst(list_of_F[j], par(1), (@y^( (s^nv-1) div (s-1) )));
    }
  }

  ring temp_ring = (char(original_ring),@y),@X,lp;
  list list_of_F = imap(ring_for_substitution, list_of_F);

  ring ring_for_factorization = (char(original_ring),@y),X,lp;
  execute("minpoly="+Z_get_minpoly(size(original_ring)^nv, parstr(1))+";");
  map rho = temp_ring,X;
  list list_of_F = rho(list_of_F);
  poly G = 0;
  for(i=1; i<=r; i++){
    G = gcd(G, list_of_F[i]);
  }
  if(G==1){
    return(list());
  }

  list factors = Presolve::linearpart(factorize(G,1));

  ideal check;
  for(i=1; i<=nv; i++){
    check[i] = X^(size(original_ring)^(i-1));
  }
  list fsols;

  matrix sc;
  list sl;
  def sM;
  matrix M_for_sol = fetch(ring_for_matrix, M_C_sol);
  for(i=1; i<=size(factors[1]); i++){
    sc = matrix(reduce(check, std(factors[1][i])), 1,nv  );

    sl = list();
    sM = sc*M_for_sol;
    for(j=1; j<=ncols(sM); j++){
      sl[j] = sM[1,j];
    }
    fsols[i] = sl;
  }
  if(size(fsols)==0){
    return(list());
  }
  setring ring_for_substitution;
  list ssols = imap(ring_for_factorization, fsols);
  if(npars(original_ring)==1){
    execute("poly P="+Z_get_minpoly(size(original_ring)^nv, "@y"));
    poly RP = gcd(P,  (@y^( (s^nv-1) div (s-1) ))-a);
    for(i=1; i<=size(ssols); i++){
      for(j=1; j<=size(ssols[i]); j++){
        ssols[i][j] = reduce( ssols[i][j], std(RP));
      }
    }
  }
  setring original_ring;
  list solutions = imap(ring_for_substitution, ssols);
  list final_solutions;
  ideal ps;
  for(i=1; i<=size(solutions); i++){
    ps = ideal();
    for(j=1; j<=nvars(original_ring); j++){
      ps[j] = var(j)-solutions[i][j];
    }
    final_solutions = insert(final_solutions, std(ps));
  }
  return(final_solutions);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = (2,a),x(1..3),lp;
  minpoly=a2+a+1;
  ideal I;
  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
  ZZsolve(I);
}
////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
static proc linearReduce(ideal I, list mons)
{
  system("--no-warn", 1);
  int LRtime = rtimer;
  int i,j ;
  int prime_field = 1;
  list solutions, monomials;
  for(i=1; i<=size(mons); i++){
    monomials = reorderMonomials(mons[i])+monomials;
  }
  int number_of_monomials = size(monomials);

  def original_ring = basering;
  if(npars(basering)==1){
    prime_field=0;
    string minpolystr = "minpoly="
    +get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
  }
  string old_vars = varstr(original_ring);
  string new_vars = "@y(1.."+string( number_of_monomials )+")";

  def ring_for_var_change = changevar( old_vars+","+new_vars, original_ring);
  setring ring_for_var_change;
  if( prime_field == 0){
    execute(minpolystr);
  }

  list monomials = imap(original_ring, monomials);
  ideal I = imap(original_ring, I);
  ideal C;
  intvec weights=1:nvars(original_ring);

  for(i=1; i<= number_of_monomials; i++){
    C[i] = monomials[i] - @y(i);
    weights = weights,deg(monomials[i])+1;
  }
  ideal linear_eqs = I;
  for(i=1; i<=ncols(C); i++){
    linear_eqs = reduce(linear_eqs, C[i]);
  }

  def ring_for_elimination = changevar( new_vars, ring_for_var_change);
  setring ring_for_elimination;
  if( prime_field == 0){
    execute(minpolystr);
  }

  ideal I = imap(ring_for_var_change, linear_eqs);
  ideal lin_sol = solvelinearpart(I);
  def ring_for_back_change = changeord( list(list("wp",weights),list("C",0:1)), ring_for_var_change);

  setring ring_for_back_change;
  if( prime_field == 0){
    execute(minpolystr);
  }

  ideal lin_sol = imap(ring_for_elimination, lin_sol);
  ideal C = imap(ring_for_var_change, C);
  ideal J = lin_sol;
  for(i=1; i<=ncols(C); i++){
    J = reduce(J, C[i]);
  }
  setring original_ring;
  ideal J = imap(ring_for_back_change, J);
  return(J);
}

static proc monomialsOfDegreeAtMost(int k)
{
  int Mtime = rtimer;
  list monomials, monoms, monoms_one, lower_monoms;
  int n = nvars(basering);
  int t,i,l,j,s;
  for(i=1; i<=n; i++){
    monoms_one[i] = var(i);
  }
  monomials = list(monoms_one);
  if(1 < k){
    for(t=2; t<=k; t++){
      lower_monoms = monomials[t-1];
      monoms = list();
      s = 1;
      for(i=1; i<=n; i++){
        for(j=s; j<=size(lower_monoms); j++){
          monoms = monoms+list(lower_monoms[j]*var(i));
        }
        s = s + int(binomial(n+t-2-i, t-2));
      }
      monomials[t] = monoms;
    }
  }
  return(monomials);
}

static proc reorderMonomials(list monomials)
{
  list univar_monoms, mixed_monoms;

  for(int j=1; j<=size(monomials); j++){
    if( univariate(monomials[j])>0 ){
      univar_monoms = univar_monoms + list(monomials[j]);
    }else{
      mixed_monoms = mixed_monoms + list(monomials[j]);
    }
  }

  return(univar_monoms + mixed_monoms);
}

static proc melyseg(poly g, list start)
{
  list gsub = g;
  int i = 1;

  while( start[1][1] <> char(basering) ){
    gsub[i+1] = subst( gsub[i], var(i), vec2elm(start[i]));
    if( gsub[i+1] == 0 ){
      list new = increment(start,i);
      for(int l=1; l<=size(start); l++){
        if(start[l]<>new[l]){
          i = l;
          break;
        }
      }
      start = new;
    }else{
      if(i == nvars(basering)){
        return(start);
      }else{
        i++;
      }
    }
  }
  return(list());
}

static proc productOfEqs(ideal I)
{
  system("--no-warn", 1);
  ideal eqs = sort_ideal(I);
  int i,q;
  poly g = 1;
  q = size(basering);
  ideal I = defaultIdeal();

  for(i=1; i<=size(eqs); i++){
    if(g==0){return(g);}
    g = reduce(g*(eqs[i]^(q-1)-1), I);
  }
  return( g );
}

static proc clonering(list #)
{
  def original_ring = basering;
  int n = nvars(original_ring);
  int prime_field=npars(basering);
  if(prime_field){
    string minpolystr = "minpoly="+
    get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
  }

  if(size(#)){
    int newvars = #[1];

  }else{
    int newvars = nvars(original_ring);
  }
  string newvarstr = "v(1.."+string(newvars)+")";
  def newring = changevar(newvarstr, original_ring);
  setring newring;
  if( prime_field ){
    execute(minpolystr);
  }
  return(newring);
}

static proc defaultIdeal()
{
  ideal I;
  for(int i=1; i<=nvars(basering); i++){
    I[i] = var(i)^size(basering)-var(i);
  }
  return( std(I) );
}

static proc order_of_extension()
{
  int oe=1;
  list rl = ringlist(basering);
  if( size(rl[1]) <> 1){
    oe = deg( subst(minpoly,par(1),var(1)) );
  }
  return(oe);
}

static proc vec2elm(intvec v)
{
  number g = 1;
  if(npars(basering) == 1){ g=par(1); }
  number e=0;
  int oe = size(v);
  for(int i=1; i<=oe; i++){
    e = e+v[i]*g^(oe-i);
  }
  return(e);
}

static proc increment(list l, list #)
{
  int c, i, j, oe;
  oe = order_of_extension();
  c = char(basering);

  if( size(#) == 1 ){
    i = #[1];
  }else{
    i = size(l);
  }

  l[i] = nextVec(l[i]);
  while( l[i][1] == c && i>1 ){
    l[i] = 0:oe;
    i--;
    l[i] = nextVec(l[i]);
  }
  if( i < size(l) ){
    for(j=i+1; j<=size(l); j++){
      l[j] = 0:oe;
    }
  }
  return(l);
}

static proc nextVec(intvec l)
{
  int c, i, j;
  i = size(l);
  c = char(basering);
  l[i] = l[i] + 1;
  while( l[i] == c && i>1 ){
    l[i] = 0;
    i--;
    l[i] = l[i] + 1;
  }
  return(l);
}

static proc every_vector()
{
  list element, list_of_elements;

  for(int i=1; i<=nvars(basering); i++){
    element[i] = 0:order_of_extension();
  }

  while(size(list_of_elements) < size(basering)^nvars(basering)){
    list_of_elements = list_of_elements + list(element);
    element = increment(element);
  }
  for(int i=1; i<=size(list_of_elements); i++){
    for(int j=1; j<=size(list_of_elements[i]); j++){
      list_of_elements[i][j] = vec2elm(list_of_elements[i][j]);
    }
  }
  return(list_of_elements);
}

static proc num2int(number a)
{
  int N=0;
  if(order_of_extension() == 1){
    N = int(a);
    if(N<0){
      N = N + char(basering);
    }
  }else{
    ideal C = coeffs(subst(a,par(1),var(1)),var(1));
    for(int i=1; i<=ncols(C); i++){
      int c = int(C[i]);
      if(c<0){ c = c + char(basering); }
      N = N + c*char(basering)^(i-1);
    }
  }
  return(N);
}

static proc get_minpoly_str(int size_of_ring, string parname)
{
  def original_ring = basering;
  ring new_ring = (size_of_ring, A),x,lp;
  string S = string(minpoly);
  string SMP;
  if(S=="0"){
    SMP = SMP+parname;
  }else{
    for(int i=1; i<=size(S); i++){
      if(S[i]=="A"){
        SMP = SMP+parname;
      }else{
        SMP = SMP+S[i];
      }
    }
  }
  return(SMP);
}

static proc sort_ideal(ideal I)
{
  ideal OI;
  int i,j,M;
  poly P;
  M = ncols(I);
  OI = I;
  for(i=2; i<=M; i++){
    j=i;
    while(size(OI[j-1])>size(OI[j])){
      P = OI[j-1];
      OI[j-1] = OI[j];
      OI[j] = P;
      j--;
      if(j==1){ break; }
    }
  }
  return(OI);
}

static proc add_if_new(list L, ideal I)
{
  int i, newelement;
  poly P;

  I=std(I);
  for(i=1; i<=nvars(basering); i++){
    P = P +  reduce(var(i),I)*var(1)^(i-1);
  }
  newelement=1;
  for(i=1; i<=size(L); i++){
    if(L[i]==P){
      newelement=0;
      break;
    }
  }
  if(newelement){
    L = insert(L, P);
  }
  return(L,newelement);
}

static proc Z_get_minpoly(int size_of_ring, string parname)
{
  def original_ring = basering;
  ring new_ring = (size_of_ring, A),x,lp;
  string S = string(minpoly);
  string SMP;
  if(S=="0"){
    SMP = SMP+parname;
  }else{
    for(int i=1; i<=size(S); i++){
      if(S[i]=="A"){
        SMP = SMP+parname;
      }else{
        SMP = SMP+S[i];
      }
    }
  }
  return(SMP);
}

static proc Z_phi(ideal I)
{
  poly f;
  for(int i=1; i<= ncols(I); i++){
    f = f+I[i]*@y^(i-1);
  }
  return(f);
}

static proc Z_default_ideal(int number_of_variables, int q)
{
  ideal DI;
  for(int i=1; i<=number_of_variables; i++){
    DI[i] = var(i)^q-var(i);
  }
  return(std(DI));
}