Return to taka_poly.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

poly_solve_linear(Eq,V) uses the lexicographic order reverse(V) and poly_solve_linear(Eq,V | reverse=1)  uses the lexicographic order V is used.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_poly.rr,v 1.40 2022/07/14 10:29:29 takayama Exp $ */
#include "tags.h"

/*  tests
 F=[x^2+y^2+z^2-1,x*y+y*z+z*x-1,x*y*z-1];
 poly_grobner_basis(F);
 poly_grobner_basis(F | order=[10,3,0], v=[x,y,z]);
 poly_grobner_basis(F | order=[10,3,0], v=[x,y,z], str=1);
 poly_initial(F); 
 poly_initial(F | gb=1); 
 poly_initial(F | order=[1,1,0]);
 poly_in(F);
 F2=[x^2+y^2-1,x*y-1]; 
 poly_hilbert_polynomial(F2);
 poly_in(F2 | v=[x], gb=1);
 poly_initial_coefficients(F2 | v=[x],gb=1);

 poly_in_w([x^2+y^2-1,x*y-x] | v=[x,y],weight=[1,0]);
*/

def taka_poly_solve_poly1(F,V) {
  if (type(F) == 3) { /* Rational */
    G = dn(F); 
    if (deg(G,V) != 0) {
      return quote(root_of(F));
    }
  }
  if (type(getopt(num))>0) {
    F=base_replace(F,[[V,x]]);
    Ans=pari(roots,F);
    if (Ans==0) return([0]);
    else return(vtol(Ans));
  }
  Ans = [ ] ;
  G = nm(F);
  F = fctr(G);
  N = length(F);
  for (I=0; I<N; I++) {
    if (deg(F[I][0],V) == 0) {
    }else if (deg(F[I][0],V) == 1) {
       Sol = red( -coef(F[I][0],0,V)/coef(F[I][0],1,V) );
       for (J=0; J<F[I][1]; J++) {
         Ans = append(Ans,[ Sol ]);
       }
    }else {
       /* Sol = quote(root_of( F[I][0] )); */
       Sol = F[I][0];
       for (J=0; J<F[I][1]; J++) {
         Ans = append(Ans,[ Sol ]);
       }
    }
  }
  return Ans;
}

def taka_poly_solve_linear(Z,X) {
  if (type(getopt(p))>0) P=getopt(p);
  else P=0;
  if (type(getopt(reverse))>0) MyReverse=getopt(reverse);
  else MyReverse=0;

  Size_Z = length(Z);
  Sol = [ ];
  F = [];
  for (I = 0; I < Size_Z; I++) {
    if (Z[I] != 0) {
      if (type(Z[I]) == 3) {
        F = append(F,[ nm(red(Z[I])) ]);
      }else{
        F = append(F, [Z[I]]);
      }
    }
  }

  if (F == []) {
    return Sol;
  }

  if (P) F=map(ptozp,F);
  if (MyReverse) {
    G = nd_gr(F, X, P, 2);
  }else{
    G = nd_gr(F, reverse(X),P, 2);
  }

  if (G[0] == 1 || G[0] == -1) {
    return [];
  }

  Sol = newvect(length(X));
  N = 0;
  if (MyReverse) {
    for (I = 0; I < length(G); I++) {
      for (J=0; J < length(X) ; J++) {
        Nm = nm(G[I]); Dn = dn(G[I]);
        if (coef(Nm, 1, X[J]) != 0) {
          if (P) {
            Sol[J] = [X[J], X[J] - Nm*inv(coef(Nm, 1, X[J]),P)];
          }else{
            Sol[J] = [X[J], X[J] - Nm/coef(Nm, 1, X[J])];
          }
          N++;
          break;
        }
      }
    }
  }else{
    for (I = 0; I < length(G); I++) {
      for (J = length(X) - 1; J >= 0; J--) {
        Nm = nm(G[I]); Dn = dn(G[I]);
        if (coef(Nm, 1, X[J]) != 0) {
          if (P) {
            Sol[J] = [X[J], X[J] - Nm*inv(coef(Nm, 1, X[J]),P)];
          }else{
            Sol[J] = [X[J], X[J] - Nm/coef(Nm, 1, X[J])];
          }
          N++;
          break;
        }
      }
    }
  }
  Sol2 = newvect(N);
  I = 0;
  for (J=0; J<size(Sol)[0]; J++) {
    if (Sol[J] == 0) {
    }else{
       Sol2[I] = Sol[J]; I++;
    }
  }
  return vtol(Sol2);
}

def taka_poly_factor(S) {
  if (type(S) == LIST || type(S) == MATRIX || type(S) == VECTOR) {
    return( map(taka_poly_factor,S) );
  }
  if (type(S) == RPOLYNOMIAL) {
     S1 = new_poly_factored_polynomial();
     S1->F = fctr(S);
     return S1;
  }
  if (type(S) == RATIONAL) {
     NN = nm(S); DD = dn(S);
     NN = taka_poly_factor(NN);
     DD = taka_poly_factor(DD);
     S = new_poly_factored_rational();
     S->Numerator = NN;
     S->Denominator = DD;
     return S;
  }
  if (type(S) == QUOTE) {
    A = quote_to_obj(S);
    B = taka_poly_factor(A);
    return B; /* object is not translated to a quote. */
  }
  /* It has not yet been implemented. */
  return S;
}

/* 
  printing method for POLY_FACTORED_FORM 
*/
def taka_poly_input_form_poly_factored_polynomial(S) {
  A="new__poly_factored_polynomial(";
  A += print_input_form(S->F);
  A += ")";
  return A;
}
def taka_poly_input_form_poly_factored_rational(S) {
  A="new__poly_factored_rational(";
  A += print_input_form(S->Numerator);
  A += ",";
  A += print_input_form(S->Denominator);
  A += ")";
  return A;
}

def taka_poly_terminal_form_poly_factored_polynomial(S) {
  A="{";
  A += "F->";
  A += print_terminal_form(S->F);
  A += "}";
  return A;
}
def taka_poly_terminal_form_poly_factored_rational(S) {
  A="{";
  A += "Numerator->";
  A += print_terminal_form(S->Numerator);
  A += ", Denominator->";
  A += print_terminal_form(S->Denominator);
  A += "}";
  return A;
}

def taka_poly_tex_form_poly_factored_polynomial(S,Tb) {
  return taka_tex_form_quote( quote_factored_form_to_quote(S),Tb);
}

/* Printing and input method
   for POLY_RING
*/
def taka_poly_input_form_poly_ring(S) {
  A  = "new__poly_ring(";
  A += taka_input_form(S->Variables)+",";
  A += taka_input_form(S->Order)+",";
  A += taka_input_form(S->K)+",";
  A += taka_input_form(S->Weyl);
  A += ")";
  return A;
}
def taka_poly_terminal_form_poly_ring(S) {
  A  = "{";
  A += "Variables->"+taka_terminal_form(S->Variables)+",";
  A += "Order->"+taka_terminal_form(S->Order)+",";
  A += "K->"+taka_terminal_form(S->K)+",";
  A += "Weyl->"+taka_terminal_form(S->Weyl)+"}";
  return A;
}

def taka_poly_tex_form_poly_ring(S,Tb) {
  taka_tex_form(S->K,Tb);
  taka_tex_form(S->Variables,Tb);
  write_to_tb(" \\mbox{\ where the order is defined by \ } ",Tb);
  taka_tex_form(S->Order,Tb);
}

/* Printing and input method
   for POLY_POLYNOMIAL
*/
def taka_poly_input_form_poly_polynomial(S) {
  A  = "new__poly_polynomial(";
  A += taka_input_form(S->Ring)+",";
  A += taka_input_form(S->F)+")";
  return A;
}
def taka_poly_terminal_form_poly_polynomial(S) {
  A  = "{";
  /* A += "Ring->"+taka_terminal_form(S->Ring)+","; */
  R = S->Ring;
  P = quote_sort_polynomial(S->F,R->Variables,R->Order);
  A += "F->"+taka_terminal_form(P)+"}";
  return A;
}

def taka_poly_tex_form_poly_polynomial(S,Tb) {
  R = S->Ring;
  P = quote_sort_polynomial(S->F,R->Variables,R->Order);
  taka_tex_form(P,Tb);
}

/* Printing and input method 
   for POLY_IDEAL
*/
def taka_poly_input_form_poly_ideal(S) {
  A  = "new__poly_ideal(";
  A += taka_input_form(S->Generators)+",";
  A += taka_input_form(S->Ring)+",";
  A += taka_input_form(S->Grobner);
  A += ")";
  return A;
}
def taka_poly_terminal_form_poly_ideal(S) {
  R = S->Ring;
  P = map(quote_sort_polynomial,S->Generators,R->Variables,R->Order);
  A  = "{";
  A += "Generators->";
  A += taka_terminal_form(P);
  A += "}";
  return A;
}

def taka_poly_tex_form_poly_ideal(S,Tb) {
  R = S->Ring;
  P = map(quote_sort_polynomial,S->Generators,R->Variables,R->Order);
  taka_tex_form(P,Tb);
}

/* Utility function to get variables. */
def taka_poly_get_variables(L) {
  if (type(L) == LIST) {
    N = length(L);
    V = [ ];
    for (I=0; I<N; I++) {
      V = base_set_union(V,vars(L[I]));
    }
    return V;
  }else if (type(L) == STRUCT) {
    if (struct_type(L) == POLY_IDEAL) {
      return L->Ring->Variables;
    }
  }else if (type(L) == RPOLYNOMIAL) {
    return vars(L);
  }else if (L == 0) {
    return [ ];
  }
  return "???";
}

/*
  taka_poly_elimination_ideal(I,V,VV|grobner_basis=yes);
 --> new interface
  taka_poly_elimination_ideal(I,VV | v=V, grobner_basis=yes);
*/
def taka_poly_elimination_ideal(I,VV) {
  /* option grobner */
  Gb = getopt(grobner_basis);
  if (type(Gb) == -1) {
    Gb = 0;
  }else{
    Gb = 1;
  }
  if (type(getopt(gb))>0) Gb=1; else Gb=0;
  /* option v */
  V = getopt(v);
  if (type(V) == -1) {
    V = taka_poly_get_variables(I);
  }
  Strategy=1;
  if (type(getopt(strategy)) >= 0); Strategy=getopt(strategy);

  II = I;
  if (type(I) == STRUCT) {
    if (struct_type(I) == POLY_IDEAL) {
      I = II->Generators;
      V = II->Ring->Variables;
      if (II->Grobner) {
         Gb = 1;
      }else {
         Gb = 0;
      }
    }else{
      error("Argument should be ideal.");
    }
  }

  if ((Strategy==0) || Gb) {
   return(taka_poly_elimination_ideal_bare(I,V,VV,Gb | option_list=getopt()));
  }else {
   return(noro_poly_elimination_ideal(I,V,VV | option_list=getopt()));
  }
}

def noro_poly_elimination_ideal(B,V,W) {
  /* compute GB of <<B>cap Q[W]> w.r.t. grlex */
  /* trace=0 => nd_gr */
  /* trace=1 => nd_gr_trace */
  /* trace=1, homo=1 => nd_gr_trace with homogenization */
	if ( type(Homo=getopt(homo)) == -1 ) Homo = 0;
	if ( type(Trace=getopt(trace)) == -1 ) Trace = 0;
	X = setminus(V,W);  /* defined in bfct */
	Ord = [[0,length(X)],[0,length(W)]];
	if ( Trace ) {
		G0 = nd_gr_trace(B,V,Homo,1,0);
		dp_ord(0); HC = map(dp_hc,map(dp_ptod,G0,V));
		G1 = ndbf.nd_gb_candidate(G0,append(X,W),Ord,Homo,HC,0);
	} else {
		G0 = nd_gr(B,V,0,0);
		G1 = nd_gr(G0,append(X,W),0,Ord);
	}
	G = ndbf.elimination(G1,W);
	return G;
}

def taka_poly_elimination_ideal_bare(I,V,VV,Gb) {
  if (!Gb) {
     V2 = base_set_minus(V,VV);
     V2 = append(V2,VV);
     /* Computing a lexicographic GB. */
     I = nd_gr(I,V2,0,2);
  }

  S = map(taka_poly_elimination_poly,I,V,VV);
  S = base_prune(0,S);
  if (type(II) == STRUCT) {
    SS = new__poly_ideal(new__poly_ring(VV,2,II->Ring->K,II->Weyl),S,0);
    return SS;
  }else{
    return S;
  }
}

def taka_poly_elimination_poly(F,V,VV) {
  U = vars(F);
  U = base_intersection(U,V);
  if (base_subsetq(U,VV)) {
     return F;
  }else{
     return 0;
  }
}

def taka_poly_grobner_basis(I) {
  Opt = getopt();
  Opt = base_rebuild_opt(Opt | remove_keys=["order","v","str"]);
  II = I;
  /* option order */
  Or = getopt(order);
  if (type(Or) == -1) {
    Or = 0;
  }
  /* option v */
  V = getopt(v);
  if (type(V) == -1) {
    if (type(I) == LIST && length(I) > 0 && type(I[0]) == QUOTE) {
       I = quote_to_obj(I); QuoteIn = 1;
       V = taka_poly_get_variables(I);
    }else{
       V = taka_poly_get_variables(I);
    }
  }
  Str = getopt(str);
  if (type(Str) == -1) Str=0;

  Gb = 0;
  if (type(I) == STRUCT) {
    if (struct_type(I) == POLY_IDEAL) {
      I = II->Generators;
      V = II->Ring->Variables;
      Or = II->Ring->Order;
      if (II->Grobner) {
         Gb = 1;
      }else {
         Gb = 0;
      }
      R = II->Ring;
    }else{
      error("Argument should be ideal.");
    }
  }else{
    R = new__poly_ring(V,Or,base_Q,0);
  }
  if (!Gb) {
    Or = taka_poly_build_order(Or);
    /* print(V); print(Or);*/
    G = nd_gr_trace(I,V,1,1,Or | option_list=Opt);
    if (QuoteIn) G = quote_to_quote(G);
  }else{
    G = II->Generators;
  }
  if (Str) SS = new__poly_ideal(R,G,1);
  else SS = G;
  return SS;
}
/* Or=[[1,1,0,0]] weight vector,
   Number
*/
def taka_poly_build_order(Or) {
  if (type(Or) <= 1) return(Or);
  /* Block order */
  if ((type(Or) == 4) && (rtostr(Or[0]) == "block")) return(cdr(Or));
  /* weight vectors. Matrix order */
  if (type(Or) == 6) { /* matrix */
    Or=matrix_matrix_to_list(Or);
  } else if (type(Or) == 5) { /* vector */
    Or=[matrix_matrix_to_list(Or)];
  } else if ((type(Or) == 4) && (type(Or[0]) != 4)) {
    Or = [Or];
  }
  N = length(Or[0]);
  M = newmat(N,N);
  for (J=0; J<N; J++) M[0][J] = 1;
  for (I=1; I<N; I++) M[I][N-I] = -1; 
  M = matrix_matrix_to_list(M);
  M2 = append(Or,M);
  M2 = newmat(length(M2),N,M2);
  return(M2);
}

/* old version 1 */
def taka_poly_grobner_basis_v1(I) {
  II = I;
  /* option order */
  Or = getopt(order);
  if (type(Or) == -1) {
    Or = 0;
  }
  /* option v */
  V = getopt(v);
  if (type(V) == -1) {
    if (type(I) == LIST && length(I) > 0 && type(I[0]) == QUOTE) {
       I = quote_to_obj(I); QuoteIn = 1;
       V = taka_poly_get_variables(I);
    }else{
       V = taka_poly_get_variables(I);
    }
  }
  Gb = 0;
  if (type(I) == STRUCT) {
    if (struct_type(I) == POLY_IDEAL) {
      I = II->Generators;
      V = II->Ring->Variables;
      Or = II->Ring->Order;
      if (II->Grobner) {
         Gb = 1;
      }else {
         Gb = 0;
      }
      R = II->Ring;
    }else{
      error("Argument should be ideal.");
    }
  }else{
    R = new__poly_ring(V,Or,base_Q,0);
  }
  if (!Gb) {
    Or = taka_poly_legacy_order(Or,V);
    /* print(V); print(Or); */
    G = dp_gr_main(I | v=V, order=Or);
    if (QuoteIn) G = quote_to_quote(G);
  }else{
    G = II->Generators;
  }

  SS = new__poly_ideal(R,G,1);
  return SS;
}

def taka_poly_legacy_order(Ord,V) {
  if (type(Ord) > 1) return Ord;
  if (Ord == 0) return [append([@grlex],V)];
  else if (Ord == 1) return [append([@glex],V)];
  else if (Ord == 2) return [append([@lex],V)];
  error("Unknown order.");
}

def taka_poly_hilbert_polynomial(I) {
  /* option grobner */
  S = getopt(s);
  if (type(S) == -1) {
    S = h;
  }
  /* option v */
  V = getopt(v);
  if (type(V) == -1) {
    V = taka_poly_get_variables(I);
  }
  Gb = 0;
  if (type(getopt(gb)) > 0) Gb=1;
  if (type(getopt(grobner_basis)) > 0) Gb=1;
  /* printf("Gb=%a\n",Gb); */
  if (type(I) == STRUCT) {
   if (struct_type(I) == POLY_IDEAL) {
     if ((I->Grobner) && (I->Ring->Order == 0)) {
       Gb = 1;
     }
     G = I->Generators;
   }
  }else {
     G = I;
  }
  if (!Gb) {
    C=poly_initial(G | v=V,order=0); 
  }else{
    C=poly_initial(G | v=V,order=0,gb=1);
  }
  if (!base_is_asir2018() || (type(getopt(sm1))>0)) {
    H = sm1.hilbert([C,V]);
    F = H-subst(H,h,h-1);
    H = subst(H,h,S);
    return [H,subst(F,h,S)];
  }
  A = dp_monomial_hilbert_poincare(C,V);
  return base_replace([not_implemeted_yet,A[3],A[2],A[0],A[1]],[[t,S],[n,S]]);
}

def taka_poly_initial_aux(I) {
  II=I;
  Gb=0; Or=0;
  if (type(getopt(gb)) > 0) Gb=1;
  if (type(getopt(grobner_basis)) > 0) Gb=1;
  if (type(I) == STRUCT) {
    if (struct_type(I) == POLY_IDEAL) {
      I = II->Generators;
      V = II->Ring->Variables;
      Or = II->Ring->Order;
      if (II->Grobner) {
         Gb = 1;
      }else {
         error("The ideal must be a Grobner basis.");
      }
      R = II->Ring;
    }else{
      error("Argument should be ideal.");
    }
  }else{
    /* weight or order is used */
    if (type(getopt(weight)) >= 0) Or=getopt(weight); 
    if (type(getopt(order)) >= 0) Or= getopt(order);
    if (type(getopt(v)) > 0) V=getopt(v);
    else V=taka_poly_get_variables(I);
  }
  Opt = base_rebuild_opt(getopt() | remove_keys=["gb","grobner_basis"]);
  if (!Gb) {
    G = poly_grobner_basis(I | option_list=Opt);
  }else G=I;

  dp_ord(taka_poly_build_order(Or));
  G2 = map(dp_ptod,G,V);
  return([G2,V,Or]);
}
def taka_poly_initial(I) {
  A = taka_poly_initial_aux(I | option_list=getopt());
  G2=A[0]; V=A[1];
  G3 = map(dp_ht,G2); 
  G4 = map(dp_dtop,G3,V);
  SS = G4;
  /* SS = new__poly_ideal(R,G4,1); */
  return SS;
}

def taka_poly_initial_coefficients(I) {
  A = taka_poly_initial_aux(I | option_list=getopt());
  G2=A[0]; V=A[1];
  G3 = map(dp_hc,G2); 
  SS = G3;
  return SS;
}

def taka_poly_initial_term(F) {
  V = getopt(v);
  Order = getopt(order);
  Weight = getopt(weight);
  if (type(F) == STRUCT) {
    if (struct_type(F) == POLY_POLYNOMIAL) {
      VV = F->Ring->Variables;
      OOrder = F->Ring->Order;
      F = F->F;
    }else{
      error("Invalid structure as an argument.");
    }
  }
  if (type(V) == -1) {
    if (VV == 0) {
      V = taka_poly_get_variables(F);  
    }else {
      V = VV;
    }
  }
  if (type(Order) == -1) {
    if (OOrder == 0) {
      Order = 0;
    }else{
      Order = OOrder;
    }
  }
  if (type(Weight) == -1) {
    Weight = 0;
  }
  if (Weight == 0) {
     dp_ord(Order);
     G = dp_ptod(F,V);
     T = dp_hm(G);

     return [dp_dtop(T,V)];
  }else {
     dp_ord(taka_poly_weight_vector(Weight,V));
     G = dp_ptod(F,V);
     C = 0;  T = 0;
     D = dp_ht(G);
     D = taka_poly_degree(D|weight=Weight,v=V);
     while (G != 0) {
       H = dp_ht(G);
       if (D > taka_poly_degree(H|weight=Weight,v=V)) break;
       T += dp_hm(G);
       G = G-dp_hm(G);
     }
     return [dp_dtop(T,V),D];
  }
}

def taka_poly_degree(F) {
  if (F == 0) return 0;
  V = getopt(v);
  Weight = getopt(weight);
  if (type(F) == STRUCT) {
    if (struct_type(F) == POLY_POLYNOMIAL) {
      VV = F->Ring->Variables;
      F = F->F;
    }else{
      error("Invalid structure as an argument.");
    }
  }
  if (type(V) == -1) {
    if (VV == 0) {
      if (type(F) != DPOLYNOMIAL) {
        V = taka_poly_get_variables(F);  
      }else{
        V = dp_etov(F);
        V = vtol(V); /* Dummy */
      }
    }else {
      V = VV;
    }
  }
  N = length(V);
  if (type(Weight) == -1) {
    Weight = newvect(N);
    for (I=0; I<N; I++) Weight[I] = 1;
  }else{
    Weight = taka_poly_weight_vector_parse(Weight,V);
  }

  if (F == 0) return 0;
  if (type(F) != DPOLYNOMIAL) {
    dp_ord(taka_poly_weight_vector(Weight,V));
    F = dp_ptod(F,V);
  }
  if (type(Weight) != VECTOR) {
    Weight = newvect(N,Weight);
  }
  D = matrix_inner_product(dp_etov(F),Weight);
  return D;
}

/*
  taka_poly_weight_vector(W,V) returns "order".
  cf. dp_ord().
  Example:
  taka_poly_weight_vector([1,0,0],[x,y,z]);
  taka_poly_weight_vector([[1,1,1],[1,-1,0]],[x,y,z]);
  taka_poly_weight_vector([x,1,y,1],[x,y,z]);
*/
def taka_poly_weight_vector(W,V) {
  if (type(W) == VECTOR || type(W) == MATRIX) {
     W = omatrix_mtol(W);
  }else if (type(W) == NUMBER || type(W) == 0) {
     return W;
  }    
  if (type(W[0]) == LIST) {
    return taka_poly_weight_vectors(W,V);
  }else{
    return taka_poly_weight_vector1(W,V);
  }
  error("Invalid argument for taka_poly_weight_vector");
}

def taka_poly_weight_vector1(W,V) {
  WW = taka_poly_weight_vector_parse(W,V);
  N = length(V);
  WW = newmat(N+2,N,[WW]);
  for (I=0; I<N; I++) WW[1][I] = 1;
  for (I=2; I<N+2; I++) {
    WW[I][N-(I-1)] = -1;
  }
  return WW;
}

def taka_poly_weight_vectors(W,V) {
  L = length(W);
  WW = [ ];
  for (I=0; I<L; I++) {
    WW = append(WW,[taka_poly_weight_vector_parse(W[I],V)]);
  }
  N = length(V);
  M = N+L;
  WW = newmat(M+1,N,WW);
  for (I=0; I<N; I++) WW[L][I] = 1;
  J = 1;
  for (I=L+1; I<M+1; I++) {
    WW[I][N-J] = -1; J++;
  }
  return WW;
}

/* w = [1,0,0,1], or [x,1,y,1] */
def taka_poly_weight_vector_parse(W,V) {
  N = length(V);
  M = length(W);
  Ans = newvect(N);
  Type = 0;
  for (I=0; I<M; I++) {
    if (type(W[I]) == STRING) {
      W[I] = eval_str(W[I]);
    }
    if (type(W[I]) == RPOLYNOMIAL) {
       Type = 1;
       K = base_position(W[I],V);
       if (K >=0 && K < N) {
          Ans[K] = W[I+1];
       }else{
         print("W=",0); print(W);
         print("V=",0); print(V);
         error("weight vector format error.");
       }
    }else if (!Type) {
       if (N != M) {
         error("weight vector format error.");
         print("W=",0); print(W);
         print("V=",0); print(V);
       }
       Ans[I] = W[I];
    }
  }
  return vtol(Ans);
}

/* V=[x,y], W=[1,1] --> [[x,1,y,1]] */
def taka_poly_weight_list(V,W) {
  L = [];
  N = length(V);
  for (I=0; I<N; I++) {
    L = cons(V[I],L);
    L = cons(W[I],L);
  }
  return [reverse(L)];
}
/* V=[x,y], W=[1,1] --> [[x,1],[y,1]] */
def taka_poly_weight_list2(V,W) {
  L = [];
  N = length(V);
  for (I=0; I<N; I++) {
    L = cons([V[I],W[I]],L);
  }
  return reverse(L);
}

/*
 The block matrix to find a GB in R=K(x)<dx>, |x|=N  
 or  R=K(x)[y], |x|=|y|=N
*/
def taka_poly_dblock(N) {
  M=newmat(2+N*2,2*N);
  for (I=N; I<2*N; I++) M[0][I]=1;
  for (I=2*N-1; I>=N; I--) M[1+(2*N-1)-I][I]= -1;
  for (I=0; I<N; I++) M[N+1][I]=1;
  for (I=N-1; I>=0; I--) M[N+1+1+(N-1)-I][I] = -1;
  return M;  
}

def taka_poly_gr_w(F,V,W) {
    G=poly_grobner_basis(F | v=V, order=[W],str=1);
    return(G->Generators);
}
def taka_poly_in_w(F) {
  if (type(F) == LIST) {
    A = taka_poly_initial_aux(F | option_list=getopt());
    V = A[1];
    W= A[2];
    G = map(dp_dtop,A[0],V);
    L=[];
    for (; G != []; G=cdr(G)) {
       L=cons(poly_initial_term(car(G) | v=V, weight=W)[0],L);
    }
    return reverse(L);
  }
  W = 0;
  if (getopt(v) > 0) V=getopt(v);
  else V = taka_poly_get_variables(F);
  if (type(getopt(weight)) >= 0) W = getopt(weight);
  if (type(getopt(order)) >= 0) W=getopt(order); 
  /* printf("W=%a, type(W)=%a\n",W,type(W)); */
  G=poly_initial_term(F | v=V, weight=W)[0];
  return(G);
}

def taka_poly_ord_w(F,V,W) {
  if (type(F) == LIST) {
    G=F;
    L=[];
    for (; G != []; G=cdr(G)) {
       L=cons(poly_initial_term(car(G) | v=V, weight=W)[1],L);
    }
    return reverse(L);
  }
  G=poly_initial_term(F | v=V, weight=W)[1];
  return(G);
}

def taka_poly_coefficient(F,Deg,V) {
  if  (type(F) < 2) return 0;
  if  (type(F) == 2) return coef(F,Deg,V);
  if  (type(F) == 3) {
    DD=dn(F); NN=nm(F);
    if (deg(DD,V) > 0) {
      printf("Denominator contains the variable %a\n",V);
      error("The denominator DD contains the specified variable V");
    }
    return red(taka_poly_coefficient(NN,Deg,V)/DD);
  }
  return map(taka_poly_coefficient,F,Deg,V);
}

def taka_poly_lcm(F) {
  G=1;
  N=length(F);
  for (I=0; I<N; I++) {
    if (F[I] != 0) G=red(G*F[I]/gcd(G,F[I]));
  }
  return(G);
}
def taka_poly_denominator(F) {
  if (type(F) == 0) return(1);
  if (type(F) < 4) return(dn(F));
  F = base_flatten(F);
  G=map(taka_poly_denominator,F);
  return(taka_poly_lcm(base_flatten(G)));
}
def taka_poly_dvar2(V,D) {
  if (type(V)==2) {
    return eval_str(D+rtostr(V));
  }
  if (type(V)>=4) {
    return map(taka_poly_dvar2,V,D);
  }
  error("Invalid argument");
}
def taka_poly_dvar(V) {
  if (type(getopt(d))<0) D="d";
  else D=rtostr(getopt(d));
  return taka_poly_dvar2(V,D);
}

Taka_poly_loaded = 1$
end$