Changeset e2fd5d in git for Singular/LIB/sagbi.lib


Ignore:
Timestamp:
Dec 24, 2010, 9:49:04 PM (13 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
b90cc8b474180d8d55d1ffc520a7e8ae0a273c30
Parents:
9e269a0a45fc3fa26da7b381bc6593a5b4b7df4d
Message:
*levandov: canonicalform which is fraction free coeffs introduced, more on printlevel, minor bugfixes

git-svn-id: file:///usr/local/Singular/svn/trunk@13780 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/sagbi.lib

    r9e269a0 re2fd5d  
    1111SAGBI stands for 'subalgebra bases analogous to Groebner bases for ideals'.
    1212SAGBI bases provide important tools for working with finitely presented
    13 subalgebras of a polynomial ring. Note that in contrast to Groebner
     13subalgebras of a polynomial ring. Note, that in contrast to Groebner
    1414bases, SAGBI bases may be infinite.
    1515
    1616REFERENCES:
    17 Ana Bravo: Some Facts About Canonical Subalgebra Bases,
     17Ana Bravo: Some Facts About Canonical Subalgebra Bases, 
    1818MSRI Publications  51, p. 247-254
    1919
     
    136136   * 5. Compute the new algebraic relations
    137137   */
     138  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
     139  dbprint(ppl,"//Spoly-1- initialisation and precomputation"); 
    138140  def br=basering;
    139141  ideal varsBasering=maxideal(1);
     
    166168    kernNew[i-nvars(r)+nvars(br)]=leadTermsAlgebra[i]-listOfVariables[i+nvars(br)];
    167169  }
    168   //--------------- calulate kernel of Phi depending on method choosen ---------------
    169 
     170  //--------------- calculate kernel of Phi depending on method choosen ---------------
     171  dbprint(ppl,"//Spoly-2- Groebner basis computation");
    170172  attrib(kernOld,"isSB",1);
    171173  ideal kern=stdKernPhi(kernNew,kernOld,leadTermsAlgebra,method);
     174  dbprint(ppl-2,"//Spoly-2-1- ideal kern",kern);       
    172175  //-------------------------- calulate algebraic relations -----------------------
     176  dbprint(ppl,"//Spoly-3- computing new algebraic relations"); 
    173177  ideal algebraicRelations=nselect(kern,1..nvars(br));
    174178  attrib(algebraicRelationsOld,"isSB",1);
    175179  ideal algebraicRelationsNew=reduce(algebraicRelations,algebraicRelationsOld);
     180  /* canonicalizing: */
     181  algebraicRelationsNew=canonicalform(algebraicRelationsNew);
     182  dbprint(ppl-2,"//Spoly-3-1- ideal of new algebraic relations",algebraicRelationsNew);
    176183  /*        algebraicRelationsOld is a groebner basis by construction (as variable
    177184   *        ordering is
     
    193200   * K[y_1,...,y_m]->K[x_1,...,x_n], x(i)->x(i), @y(i)->leadterm(algebra[i])
    194201   */
     202  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
     203  dbprint(ppl,"//Spoly-1- initialisation and precomputation"); 
    195204  def br=basering;
    196205  int m=ncols(algebra);
     
    203212    leadCoefficients[i]=leadcoef(algebra[i]);
    204213  }
     214  dbprint(ppl-2,"//Spoly-1-1- Vector of leading coefficients",leadCoefficients);
    205215  int k=1;
    206216  for (i=1;i<=n;i++)
     
    215225  //of the leadings monomials in algebra.
    216226  intmat A[n][m]=intmat(tempVec,n,m);
     227  dbprint(ppl-2,"//Spoly-1-2- Matrix A",A);
    217228  //Create the preimage ring K[@y(1),...,@y(m)], where m=ncols(algebra).
    218229  string variableName=uniqueVariableName("@y");
     
    227238  setring rNew;
    228239  //Use toric_ideal to compute the kernel
     240   dbprint(ppl,"//Spoly-2- call of toric_ideal");
    229241  ideal algebraicRelations=toric_ideal(A,"ect");
    230242  //Suitable substitution
     243        dbprint(ppl,"//Spoly-3- substitutions");       
    231244  ideal leadCoefficients=fetch(br,leadCoefficients);
    232245  for (i=1; i<=m; i++)
     
    237250    }
    238251  }
     252        dbprint(ppl-2,"//Spoly-3-1- algebraic relations",algebraicRelations);   
    239253  export algebraicRelations;
    240254  return(rNew);
     
    251265   * p=reduce(leadF,kern) in K[@y(1),...,@y(m)] <=> leadF in K[lead(algebra)]
    252266   * Otherwise some precomputation has to be done, outlined below.
     267   * When using sagbiReduce,sagbi and sagbiPart the integer parRed will always be zero. Only the procedure
     268   * algebraicDependence causes this procedure to be called with parRed<>0. The only difference when parRed<>0
     269   * is that the reduction algorithms returns the non-zero constants it attains (instead of just returning zero as the
     270         * correct remainder), as they will be expressions in parameters for an algebraic dependence.
    253271   */
     272  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
     273  dbprint(ppl,"//Red-1- initialisation and precomputation");   
    254274  def br=basering;
    255275  int numVarsBasering=nvars(br);
     
    259279  if (numVarsBasering==nvars(r))
    260280  {
     281                dbprint(ppl-1,"//Red-1-1- Groebner basis computation");
    261282    /* Case that ring r is the same ring as the basering. Using proc extendRing,
    262283     * stdKernPhi
     
    281302    }
    282303    kern=stdKernPhi(kern,0,leadTermsAlgebra,method);
     304    dbprint(ppl-2,"//Red-1-1-1- Ideal kern",kern);
    283305  }
    284306  setring r;
     
    290312  intvec tempExp;
    291313  //-------------algebraic reduction for each polynomial F[i] ------------------------
     314  dbprint(ppl,"//Red-2- reduction, polynomial by polynomial");
    292315  for (i=1; i<=ncols(F);i++)
    293316  {
     317                dbprint(ppl-1,"//Red-2-"+string(i)+"- starting with new polynomial");
     318                dbprint(ppl-2,"//Red-2-"+string(i)+"-1- Polynomial before reduction",F[i]);
    294319    normalform=0;
    295320    while (F[i]!=0)
     
    300325      //K is always contained in the subalgebra,
    301326      //thus the remainder is zero in this case
    302         if (parRed) { break; }
    303         else { F[i]=0; break; }
     327        if (parRed)
     328        {
     329                                //If parRed<>0 save non-zero constants the reduction algorithms attains.
     330                                        break;
     331                                }
     332        else
     333        {
     334                                        F[i]=0;
     335                                        break;
     336                                }
    304337      }
    305338      //note: as the ordering in br and r might not be compatible
     
    338371      F[i] = normalform;
    339372    }
     373    dbprint(ppl-2,"//Red-2-"+string(i)+"-2- Polynomial after reduction",F[i]);
    340374  }
    341375  return(F);
     
    389423 *   and iterations specifies the
    390424 *   number of iterations. A degree boundary is not used here.
    391  * Note that parRed is used for testing a special modification
    392  * and can be ignored (assume parRed==0).
     425 * When this method is called via the procedures sagbi and sagbiPart the integer parRed
     426 * will always be zero. Only the procedure algebraicDependence calls this procedure with
     427 * parRed<>0. The only difference when parRed<>0 is that the reduction algorithms returns
     428 * the non-zero constants it attains (instead of just returning zero as the correct
     429 * remainder), as they will be expressions in parameters for an algebraic dependence.
     430 * These constants are saved in the ideal reducedParameters.
    393431 */
    394432{
     433  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
     434  dbprint(ppl,"// -0- initialisation and precomputation");     
    395435  def br=basering;
    396436  int i=1;
    397437
    398438  ideal reducedParameters;
    399   int numReducedParameters=1;
     439  int numReducedParameters=1; //number of elements plus one in reducedParameters
    400440  int j;
    401   if (parRed==0)
     441  if (parRed==0) //if parRed<>0 the algebra does not contain monomials and normalisation should be avoided
    402442  {
    403443    algebra=reduceByMonomials(algebra);
    404444    algebra=simplify(simplify(algebra,3),4);
    405445  }
    406   int step=1;
    407   if (iterations==-1) //case: infintitly many iterations
    408   {
    409     step=0;
    410     iterations=1;
    411   }
     446  // canonicalizing the gen's:
     447  algebra = canonicalform(algebra);
    412448  ideal P=1;
    413449  //note: P is initialized this way, so that the while loop is entered.
     
    417453  ideal spolynomialsNew;
    418454  def r=br;
    419   while (size(P)>0 && i<=iterations)
    420   {
    421     def rNew=spolynomialsGB(algebra,r,method);
     455  while (size(P)>0)
     456  {
     457    dbprint(ppl,"// -"+string(i)+"- interation of SAGBI construction algorithm");
     458    dbprint(ppl-1,"// -"+string(i)+"-1- Computing algebraic relations");
     459    def rNew=spolynomialsGB(algebra,r,method); /* canonicalizing inside! */
    422460    kill r;
    423461    def r=rNew;
    424462    kill rNew;
    425463    phi=r,varsBasering,algebra;
     464    dbprint(ppl-1,"// -"+string(i)+"-2- Substituting into algebraic relations");
    426465    spolynomialsNew=simplify(phi(algebraicRelationsNew),6);
    427466    //By construction spolynomialsNew only contains the spolynomials,
    428467    //that have not already
    429468    //been calculated in the steps before.
     469    dbprint(ppl-1,"// -"+string(i)+"-3- SAGBI reduction");
     470    dbprint(ppl-2,"// -"+string(i)+"-3-1- new S-polynomials before reduction",spolynomialsNew);
    430471    P=reductionGB(spolynomialsNew,algebra,r,tailreduction,method,parRed);
    431472    if (parRed)
     
    452493      P=simplify(P,6);
    453494    }
     495    /* canonicalize ! */
     496    P = canonicalform(P);
     497    dbprint(ppl-2,"// -"+string(i)+"-3-1- new S-polynomials after reduction",P);
    454498    algebra=algebra,P;
    455499    //Note that elements and order of elements must in algebra must not be changed,
    456500    //otherwise the already calculated
    457501    //ideal in r will give wrong results. Thus it is important to use a komma here.
    458     i=i+step;
    459   }
    460   if (step==1)
     502    i=i+1;
     503    if (iterations!=-1 && i>iterations) //When iterations==-1 the number of iterations is unlimited
     504    {
     505      break;
     506    }
     507  }
     508  if (iterations!=-1)
    461509  { //case that sagbiPart called this procedure
    462510    if (size(P)==0)
     
    480528    algebra=algebra,reducedParameters;
    481529  }
    482   algebra=simplify(algebra,6);
     530  algebra = simplify(algebra,6);
     531  algebra = canonicalform(algebra);
    483532  return(algebra);
    484533}
     
    735784  a=sagbiConstruction(algebra,-1,tailreduction,method,0);
    736785  a=simplify(a,7);
    737   a=interreduced(a);
     786  //  a=interreduced(a);
    738787  return(a);
    739788}
     
    803852  ideal a;
    804853  a=sagbiConstruction(algebra,iterations,tailreduction,method,0);
    805   a=simplify(a,7);
    806   a=interreduced(a);
     854  a=simplify(a,6);
     855  //  a=interreduced(a);
    807856  return(a);
    808857}
     
    821870
    822871
    823 
    824 // VL: finished the documentation
    825 // TO compare with algDependent from algebra_lib
    826 // TODO : remove constants from algebraic dependencies
    827872proc algebraicDependence(ideal I,int iterations)
    828873"USAGE: algebraicDependence(I,it); I an an ideal, it is an integer
    829874RETURN: ring
    830875ASSUME: basering is not a qring
    831 PURPOSE: In @code{it} iterations, compute algebraic dependencies between elements of I
     876PURPOSE: Returns a ring containing the ideal @code{algDep}, which contains possibly
     877@*              some algebraic dependencies of the elements of I obtained through @code{it}
     878@*              iterations of the SAGBI construction algorithms. See the example on how
     879@*              to access these objects.
    832880EXAMPLE: example algebraicDependence; shows an example"
    833881{
    834   assumeQring();
     882        assumeQring();
     883        int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
     884  dbprint(ppl,"//AlgDep-1- initialisation and precomputation");
    835885  def br=basering;
    836886  int i;
    837 
     887  I=simplify(I,2); //avoid that I contains zeros
     888
     889        //Create two polynomial rings, which both are extensions of the current basering.
     890  //The first ring will contain the additional paramteres @c(1),...,@c(m), the second one
     891  //will contain the additional variables @c(1),...,@c(m), where m=ncols(I).
    838892  string parameterName=uniqueVariableName("@c");
    839893  list l = ringlist(basering);
     
    843897    parList[i]=string(parameterName,"(",i,")");
    844898  }
    845   l[1]=list(l[1],parList,list(list("dp",1:ncols(I))));
     899  l[1]=list(l[1],parList,list(list("dp",1:ncols(I)))); //add @c(i) to the ring as paramteres
    846900  ideal temp=0;
    847901  l[1][4]=temp;
     
    862916  for (i=1; i<=ncols(I);i++)
    863917  {
    864     l2[2][i+nvars(br)]=string(extendVarName,"(",i,")");
     918    l2[2][i+nvars(br)]=string(extendVarName,"(",i,")"); //add @c(i) to the rings as variables
    865919  }
    866920  l2[3][size(l2[3])+1]=l2[3][size(l2[3])];
     
    899953  //         }
    900954
    901 
     955        //Compute a partial SAGBI basis of the algebra generated by I[1]-@c(1),...,I[m]-@c(m),
     956  //where the @c(n) are parameters
    902957  ideal I=fetch(br,I);
    903958  ideal algebra;
     
    906961    algebra[i]=I[i]-par(i);
    907962  }
     963  dbprint(ppl,"//AlgDep-2- call of SAGBI construction algorithm");     
    908964  algebra=sagbiConstruction(algebra, iterations,0,0,1);
     965  dbprint(ppl,"//AlgDep-3- postprocessing of results");
    909966  int j=1;
     967  //If K[x_1,...,x_n] was the basering, then algebra is in K(@c(1),...,@c(m))[x_1,...x_n]. We intersect
     968  //elements in algebra with K(@c(1),..,@c(n)) to get algDep. Note that @c(i) can only appear in the numerator,
     969  //as the SAGBI construction algorithms just multiplies and substracts polynomials. So actually we have
     970  //algDep=algebra intersect K[@c(1),...,@c(m)]
    910971  ideal algDep;
    911972  for (i=1; i<= ncols(algebra); i++)
    912973  {
    913     if (leadmonom(algebra[i])==1)
     974    if (leadmonom(algebra[i])==1) //leadmonom(algebra[i])==1 iff algebra[i] in K[@c(1),...,@c(m)]
    914975    {
    915976      algDep[j]=algebra[i];
     
    917978    }
    918979  }
     980  //Transfer algebraic dependencies to ring where @c(i) are not parameters, but now variables.
    919981  setring extendVarRing;
    920982  ideal algDep=imap(parameterRing,algDep);
    921983  ideal algebra=imap(parameterRing,algebra);
     984  //Now get rid of constants in K that may have been added to algDep.
     985  for (i=1; i<=ncols(algDep); i++)
     986  {
     987                if(leadmonom(algDep[i])==1)
     988                {
     989                                algDep[i]=0;
     990                }
     991        }
     992        algDep=simplify(algDep,2);
    922993  export algDep,algebra;
    923   //print(algDep);
    924994  setring br;
    925995  return(extendVarRing);
     
    9461016static proc interreduced(ideal I)
    9471017{
     1018  /* performs subalgebra interreduction of a set of subalgebra generators */
     1019  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
     1020  dbprint(ppl,"//Interred-1- starting interreduction");
    9481021  ideal J,B;
    9491022  int i,j,k;
     
    9511024  for(k=1;k<=ncols(I);k++)
    9521025  {
     1026    dbprint(ppl-1,"//Interred-1-"+string(k)+"- reducing next poly");
    9531027    f=I[k];
    9541028    I[k]=0;
     
    9571031  }
    9581032  I=simplify(I,2);
     1033  dbprint(ppl,"//Interred-2- interreduction completed");       
    9591034  return(I);
    9601035}
     
    11481223}
    11491224
     1225static proc canonicalform(ideal I)
     1226{
     1227  /* placeholder for the canonical form of a set of gen's */
     1228  /* for the time being we agree on content(p)=1; that is coeffs with no fractions */
     1229  int i; ideal J=I;
     1230  for(i=ncols(I); i>=1; i--)
     1231  {
     1232    J[i] = canonicalform_poly(I[i]);
     1233  }
     1234  return(J);
     1235}
     1236
     1237static proc canonicalform_poly(poly p)
     1238{
     1239  /* placeholder for the canonical form of a poly */
     1240  /* for the time being we agree on content(p)=1; that is coeffs with no fractions */
     1241  number n = content(p);
     1242  return( p/content(p) );
     1243}
     1244
    11501245/*
    11511246  ring r= 0,(x,y),dp;
Note: See TracChangeset for help on using the changeset viewer.