Changeset 5ebd453 in git


Ignore:
Timestamp:
May 25, 2009, 6:29:04 PM (15 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
8512090f711dc406701bd340e3fd7e96f79a294c
Parents:
5b2685e36d3336e4237ec537cdbf20e8dd4f57e3
Message:
*santiago: new version


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/normal.lib

    r5b2685e r5ebd453  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: normal.lib,v 1.53 2009-04-15 11:15:56 seelisch Exp $";
     2version="$Id: normal.lib,v 1.54 2009-05-25 16:29:04 Singular Exp $";
    33category="Commutative Algebra";
    44info="
     
    166166  if ( find(method,"withDelta") or find(method,"wd") or find(method,"withdelta"))
    167167  {
    168     if((decomp == 0) or (decomp == 3)){
     168    if((decomp == 0) or (decomp == 3))
     169    {
    169170      withDelta = 1;
    170     } else {
     171    }
     172    else
     173    {
    171174      "WARNING: the delta invariants cannot be computed with an equidimensional";
    172175      "         decomposition.";
     
    191194                          // R is the ring where computations will be done
    192195
    193   if((useRing  == 1) and (isGlobal == 1)){
     196  if((useRing  == 1) and (isGlobal == 1))
     197  {
    194198    def globR = basering;
    195   } else {
     199  }
     200  else
     201  {
    196202    // We change to dp ordering.
    197203    list rl = ringlist(origR);
     
    206212//------------------------ trivial checkings ------------------------
    207213  id = groebner(id);
    208   if((size(id) == 0) or (id[1] == 1)){
     214  if((size(id) == 0) or (id[1] == 1))
     215  {
    209216    // The original ring R/I was normal. Nothing to do.
    210217    // We define anyway a new ring, equal to R, to be able to return it.
     
    218225    export normap;
    219226    setring origR;
    220     if(withDelta){
     227    if(withDelta)
     228    {
    221229      result = list(list(ideal(1), poly(1), ROut), 0);
    222230    } else {
     
    228236//------------------------ preliminary decomposition-----------------------
    229237  list prim;
    230   if(decomp == 2){
     238  if(decomp == 2)
     239  {
    231240    dbprint(dbg, "// Computing the equidimensional decomposition...");
    232241    prim = equidim(id);
    233242  }
    234   if((decomp == 0) or (decomp == 1)){
     243  if((decomp == 0) or (decomp == 1))
     244  {
    235245    prim = id;
    236246  }
    237   if(decomp == 3){
     247  if(decomp == 3)
     248  {
    238249    dbprint(dbg, "// Computing the minimal associated primes...");
    239250    if( noFac )
     
    253264//----------------- back to the original ring if required ------------------
    254265// if ring was not global and useRing is on, we go back to the original ring
    255   if((useRing == 1) and (isGlobal != 1)){
     266  if((useRing == 1) and (isGlobal != 1))
     267  {
    256268    setring origR;
    257269    def R = basering;
    258270    list prim = fetch(globR, prim);
    259   } else {
     271  }
     272  else
     273  {
    260274    def R = basering;
    261     if(useRing == 0){
     275    if(useRing == 0)
     276    {
    262277      ideal U;
    263278      poly c;
     
    283298  {
    284299    if(dbg>=2){pause();}
    285     if(dbg>=1){
     300    if(dbg>=1)
     301    {
    286302      "// start computation of component",i;
    287303      "   --------------------------------";
    288304    }
    289     if(groebner(prim[i])[1] != 1){
    290       if(dbg>=2){
     305    if(groebner(prim[i])[1] != 1)
     306    {
     307      if(dbg>=2)
     308      {
    291309        "We compute the normalization in the ring"; basering;
    292310      }
     
    294312      norComp = normalM(prim[i], decomp, withDelta);
    295313      printlevel = printlevel - 1;
    296       for(j = 1; j <= size(norComp); j++){
     314      for(j = 1; j <= size(norComp); j++)
     315      {
    297316        newR = norComp[j][3];
    298317        newRList = ringlist(newR);
    299318        U = norComp[j][1];
    300319        c = norComp[j][2];
    301         if(withDelta){
     320        if(withDelta)
     321        {
    302322          delt = norComp[j][4];
    303           if((delt >= 0) and (deltI >= 0)){
     323          if((delt >= 0) and (deltI >= 0))
     324          {
    304325            deltI = deltI + delt;
    305           } else {
     326          }
     327          else
     328          {
    306329            deltI = -1;
    307330          }
    308331        }
    309332        // -- incorporate result for this component to the list of results ---
    310         if(useRing == 0){
     333        if(useRing == 0)
     334        {
    311335          // We go back to the original ring.
    312336          setring origR;
     
    315339          newRListO = imap(R, newRList);
    316340          // We change the ordering in the new ring.
    317           if(nvars(newR) > nvars(origR)){
     341          if(nvars(newR) > nvars(origR))
     342          {
    318343            newRListO[3]=insert(origOrd, newRListO[3][1]);
    319           } else {
     344          }
     345          else
     346          {
    320347            newRListO[3] = origOrd;
    321348          }
     
    329356          totalComps++;
    330357          result[totalComps] = list(U, c, newROrigOrd);
    331           if(withDelta){
     358          if(withDelta)
     359          {
    332360            result[totalComps] = insert(result[totalComps], delt, 3);
    333361          }
    334362          setring R;
    335         } else {
     363        }
     364        else
     365        {
    336366          setring R;
    337367          totalComps++;
     
    343373
    344374// -------------------------- delta computation ----------------------------
    345   if(withDelta == 1){
     375  if(withDelta == 1)
     376  {
    346377    // Intersection multiplicities of list prim, sp=size(prim).
    347378    if ( dbg >= 1 )
     
    371402  list MG;      // Module generators
    372403  intvec DV;    // Vector of delta's of each component
    373   for(i = 1; i <= size(result); i++){
     404  for(i = 1; i <= size(result); i++)
     405  {
    374406    RL[i] = result[i][3];
    375407    MG[i] = lineUpLast(result[i][1], result[i][2]);
    376     if(withDelta){
     408    if(withDelta)
     409    {
    377410      DV[i] = result[i][4];
    378411    }
     
    380413  list resultNew;
    381414
    382   if(withDelta){
     415  if(withDelta)
     416  {
    383417    resultNew = list(RL, MG, list(DV, deltI));
    384418  } else {
     
    391425  {
    392426    "";
    393     if(!withDelta){
     427    if(!withDelta)
     428    {
    394429      "// 'normal' created a list, say nor, of two elements.";
    395430    } else {
     
    410445    "// 1/ci * Ui, with Ui the i-th ideal of nor[2] and ci the last element";
    411446    "// Ui[size(Ui)] of Ui.";
    412     if(withDelta){
     447    if(withDelta)
     448    {
    413449      "// * nor[3] is a list of an intvec of size r, the delta invariants ";
    414450      "// of the r components, and an integer, the total delta invariant ";
     
    484520   poly p = Li[4];
    485521   int noRed = 0;
    486    if(size(Li) > 4){
    487      if(Li[5] == 1){
     522   if(size(Li) > 4)
     523   {
     524     if(Li[5] == 1)
     525     {
    488526       noRed = 1;
    489527     }
     
    732770
    733771
    734   if(noRed == 0){
     772  if(noRed == 0)
     773  {
    735774    L = substpart(endid,endphi,homo,rw);
    736775    def lastRing=L[1];
     
    32483287proc normalP(ideal id,list #)
    32493288"USAGE:   normalP(id [,choose]); id a radical ideal, choose a comma separated
    3250          list of optional strings: \"withRing\" or \"isPrim\" or \"noFac\".@*
    3251          If choose = \"noFac\", factorization is avoided during the computation
    3252          of the minimal associated primes; choose = \"isPrim\" disables the
    3253          computation of the minimal associated primes (should only be used
    3254          if the ideal is known to be prime).@*
     3289         list of optional strings \"withRing\", \"isPrim\", \"noFac\",
     3290         \"noRed\", where@*
     3291         - \"noFac\", factorization is avoided during the computation
     3292         of the minimal associated primes.@*
     3293         - \"isPrim\", disables the computation of the minimal associated
     3294         primes (should only be used if the ideal is known to be prime).@*
     3295         - \"withRing\", the ring structure of the normalization is
     3296         computed. The number of variables in the new ring is reduced as much
     3297         as possible.@*
     3298         - \"noRed\", when computing the ring structure no reduction on the
     3299         number of  variables is done, it creates one new variable for each
     3300         of the new generators of the integral closure in the quotient field.@*
    32553301ASSUME:  The characteristic of the ground field must be positive. The ideal
    32563302         id is assumed to be radical. However, if choose does not contain
     
    33013347"
    33023348{
    3303    int i,j,y, wring, isPrim, noFac, sr, del, co;;
     3349   int i,j,y, sr, del, co;
    33043350   intvec deli;
    33053351   list resu, Resu, prim, Gens, mstdid;
    33063352   ideal gens;
     3353
     3354   // Default options
     3355   int wring = 0;           // The ring structure is not computed.
     3356   int noRed = 0;           // No reduction is done in the ring structure
     3357   int isPrim = 0;          // Ideal is not assumed to be prime
     3358   int noFac = 0;           // Use facstd when computing the decomposition
     3359
     3360
    33073361   y = printlevel-voice+2;
    33083362
     
    33373391     wring=1;
    33383392   }
     3393   if ( find(method,"noRed") or find(method,"nored") )
     3394   {
     3395     noRed=1;
     3396   }
    33393397   if ( find(method,"isPrim") or find(method,"isprim") )
    33403398   {
    33413399     isPrim=1;
    33423400   }
    3343    if ( find(method,"noFac") )
     3401   if ( find(method,"noFac") or find(method,"nofac"))
    33443402   {
    33453403     noFac=1;
     
    34043462         {
    34053463         //------ 2-nd main subprocedure: compute ring structure -----------
    3406             resu = list(computeRing(II,prim[i])) + resu;
     3464           if(noRed == 0){
     3465             resu = list(computeRing(II,prim[i])) + resu;
     3466           } else {
     3467             resu = list(computeRing(II,prim[i], "noRed")) + resu;
     3468           }
    34073469         }
    34083470         printlevel = printlevel-1;
     
    36423704      list g3 = gnirlist[3];             //contains blocks of orderings
    36433705      int n3 = size(g3);
     3706
    36443707   //----------------- first identify module ordering ------------------
    36453708      if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
     
    36563719      }
    36573720   //---- delete variables which were substituted and weights  --------
    3658 //gnirlist;"";
    36593721      intvec V;
    36603722      int n(0);
     
    36633725      for ( ii=1; ii<=n3-1; ii++ )
    36643726      {
    3665          V = V,g3[ii][2];           //copy weights for ordering in each block
    3666          if ( ii==1 )               //into one intvector
    3667          {
    3668             V = V[2..size(V)];
    3669          }
     3727        // If order is a matrix ordering, it is replaced by dp ordering.
     3728        // TODO: replace it only when some of the original
     3729        //       variables are eliminated.
     3730        if(g3[ii][1] == "M"){
     3731          g3[ii][1] = "dp";
     3732          g3[ii][2] = (1..sqroot(size(g3[ii][2])))*0+1;
     3733        }
     3734        V = V,g3[ii][2];           //copy weights for ordering in each block
     3735        if ( ii==1 )               //into one intvector
     3736        {
     3737           V = V[2..size(V)];
     3738        }
    36703739        // int n(ii) = size(g3[ii][2]);
    36713740        int n(ii) = size(V);
     
    36743743        for ( jj = n(ii-1)+1; jj<=n(ii); jj++)
    36753744        {
    3676 //"ii, jj", ii,jj;
    36773745          if( Le[4][jj] !=0 )     //jj=index of var which was not substituted
    36783746          {
     
    36803748            newg2[kk] = g2[jj];   //not substituted var from varlist
    36813749            V(ii)=V(ii),V[jj];    //weight of not substituted variable
    3682 //"V(",ii,")", V(ii);
    36833750          }
    36843751        }
     
    40094076  def R = basering;
    40104077
    4011   // ---------------- computation of the singular locus ---------------------
    4012   // We compute the radical of the ideal of minors modulo the original ideal.
    4013   // This is done only in the first step.
    4014 
     4078//------------------------ Groebner bases and dimension of I-----------------
    40154079  if(isGlobal == 1){
    40164080    list IM = mstd(I);
     
    40254089  int d = dim(I);
    40264090
     4091  // ---------------- computation of the singular locus ---------------------
     4092  // We compute the radical of the ideal of minors modulo the original ideal.
     4093  // This is done only in the first step.
    40274094  qring Q = I;   // We work in the quotient by the groebner base of the ideal I
    40284095  option("redSB");
     
    40334100
    40344101  dbprint(dbg, "Computing the jacobian ideal...");
     4102
    40354103  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
    40364104  J = groebner(J);
     4105
     4106  //------------------ We check if the singular locus is empty -------------
     4107  if(J[1] == 1){
     4108    // The original ring R/I was normal. Nothing to do.
     4109    // We define anyway a new ring, equal to R, to be able to return it.
     4110    setring R;
     4111    list lR = ringlist(R);
     4112    def ROut = ring(lR);
     4113    setring ROut;
     4114    ideal norid = fetch(R, I);
     4115    ideal normap = maxideal(1);
     4116    export norid;
     4117    export normap;
     4118    setring R;
     4119    if(withDelta){
     4120      list output = ideal(1), poly(1), ROut, 0;
     4121    } else {
     4122      list output = ideal(1), poly(1), ROut;
     4123    }
     4124    return(list(output));
     4125  }
     4126
    40374127
    40384128  // -------------------- election of the conductor -------------------------
     
    41484238      // I = Id1 \cap Id2
    41494239      printlevel = printlevel + 1;
     4240
    41504241      list nor1 = normalM(Id1, decomp, withDelta)[1];
    41514242      list nor2 = normalM(Id2, decomp, withDelta)[1];
     
    42104301  ideal oldU = 1;
    42114302
    4212   if(dbg >= 2){
     4303  if(dbg >= 2)
     4304  {
    42134305    "The quotient is"; U;
    42144306  }
     
    42174309  // We check if the equality in Grauert - Remmert criterion is satisfied.
    42184310  isNormal = checkInclusions(D*oldU, U);
    4219   if(isNormal == 0){
    4220     if(dbg >= 1){
     4311  if(isNormal == 0)
     4312  {
     4313    if(dbg >= 1)
     4314    {
    42214315      "In this step, we have the ring 1/c * U, with c =", c;
    42224316      "and U = "; U;
    42234317    }
    4224   } else {
     4318  }
     4319  else
     4320  {
    42254321    // The original ring R/I was normal. Nothing to do.
    42264322    // We define anyway a new ring, equal to R, to be able to return it.
     
    42294325    def ROut = ring(lR);
    42304326    setring ROut;
    4231     ideal norid = 0;
     4327    ideal norid = fetch(R, I);
    42324328    ideal normap = maxideal(1);
    42334329    export norid;
    42344330    export normap;
    42354331    setring R;
    4236     if(withDelta){
     4332    if(withDelta)
     4333    {
    42374334      list output = ideal(1), poly(1), ROut, 0;
    4238     } else {
     4335    }
     4336    else
     4337    {
    42394338      list output = ideal(1), poly(1), ROut;
    42404339    }
     
    42434342
    42444343  // ----- computation of the chain of ideals A1 c A2 c ... c An ------------
    4245   while(isNormal == 0){
     4344  while(isNormal == 0)
     4345  {
    42464346    step++;
    4247     if(dbg >= 1){
     4347    if(dbg >= 1)
     4348    {
    42484349      "";
    42494350      "Step ", step, " begins.";
     
    42704371    cJMod = getGenerators(cJ, U, c);
    42714372
    4272     if(dbg >= 2){
     4373    if(dbg >= 2)
     4374    {
    42734375      "The test ideal in this step is ";
    42744376      cJMod;
     
    42874389    isNormal = checkInclusions(D*oldU, U);
    42884390
    4289     if(isNormal == 1){
     4391    if(isNormal == 1)
     4392    {
    42904393      // We go one step back. In the last step we didnt get antyhing new,
    42914394      // we just verified that the ring was already normal.
     
    42934396      dbprint(dbg, "");
    42944397      U = oldU;
    4295     } else {
     4398    }
     4399    else
     4400    {
    42964401      // ------------- preparation for next iteration ----------------------
    42974402      // We have to go on.
    42984403      // The new denominator is chosen.
    42994404      c = D * c;
    4300       if(deg(c) > deg(condu)){
     4405      if(deg(c) > deg(condu))
     4406      {
    43014407        U = changeDenominatorQ(U, c, condu);
    43024408        c = condu;
    43034409      }
    4304       if(dbg >= 1){
     4410      if(dbg >= 1)
     4411      {
    43054412        "In this step, we have the ring 1/c * U, with c =", c;
    43064413        "and U = ";
     
    43124419
    43134420  // ------------------------- delta computation ----------------------------
    4314   if(withDelta){
     4421  if(withDelta)
     4422  {
    43154423    ideal UD = groebner(U);
    43164424    delt = vdim(std(modulo(UD, c)));
     
    43404448  int i;
    43414449  ideal newU;
    4342   for (i = 1; i <= ncols(U); i++){
    4343     if(U[i] != c){
    4344       if(size(newU) == 0){
     4450  for (i = 1; i <= ncols(U); i++)
     4451  {
     4452    if(U[i] != c)
     4453    {
     4454      if(size(newU) == 0)
     4455      {
    43454456        newU = U[i];
    4346       } else {
     4457      }
     4458      else
     4459      {
    43474460        newU = newU, U[i];
    43484461      }
    43494462    }
    43504463  }
    4351   if(size(newU) == 0){
     4464  if(size(newU) == 0)
     4465  {
    43524466    newU = c;
    4353   } else {
     4467  }
     4468  else
     4469  {
    43544470    newU = newU, c;
    43554471  }
     
    43644480  int i;
    43654481  ideal newU = c;
    4366   for (i = 1; i <= ncols(U); i++){
    4367     if(U[i] != c){
     4482  for (i = 1; i <= ncols(U); i++)
     4483  {
     4484    if(U[i] != c)
     4485    {
    43684486      newU = newU, U[i];
    43694487    }
     
    43744492///////////////////////////////////////////////////////////////////////////////
    43754493
    4376 static proc getSmallest(ideal J){
     4494static proc getSmallest(ideal J)
     4495{
    43774496// Computes the polynomial of smallest degree of J.
    43784497// If there are more than one, it chooses the one with smallest number
     
    43824501  int d = deg(p);
    43834502  int di;
    4384   for(i = 2; i <= ncols(J); i++){
    4385     if(J[i] != 0){
     4503  for(i = 2; i <= ncols(J); i++)
     4504  {
     4505    if(J[i] != 0)
     4506    {
    43864507      di = deg(J[i]);
    4387       if(di < d){
     4508      if(di < d)
     4509      {
    43884510        p = J[i];
    43894511        d = di;
    4390       } else {
    4391         if(di == d){
    4392           if(size(J[i]) < size(p)){
     4512      }
     4513      else
     4514      {
     4515        if(di == d)
     4516        {
     4517          if(size(J[i]) < size(p))
     4518          {
    43934519            p = J[i];
    43944520          }
     
    44024528///////////////////////////////////////////////////////////////////////////////
    44034529
    4404 static proc getGenerators(ideal J, ideal U, poly c){
     4530static proc getGenerators(ideal J, ideal U, poly c)
     4531{
    44054532// Computes the generators of J as an ideal in the original ring,
    44064533// where J is given by generators in the new ring.
     
    44144541
    44154542  if(dbg>1){"Checking for new generators...";}
    4416   for(i = 1; i <= ncols(J); i++){
    4417     for(j = 1; j <= ncols(U); j++){
     4543  for(i = 1; i <= ncols(J); i++)
     4544  {
     4545    for(j = 1; j <= ncols(U); j++)
     4546    {
    44184547      p = lift(c, J[i]*U[j])[1,1];
    44194548      p = reduce(p, JGr);
    4420       if(p != 0){
    4421         if(dbg>1){
     4549      if(p != 0)
     4550      {
     4551        if(dbg>1)
     4552        {
    44224553          "New polynoial added:", p;
    44234554          if(dbg>4){pause();}
     
    44344565///////////////////////////////////////////////////////////////////////////////
    44354566
    4436 static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D){
     4567static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D)
     4568{
    44374569// Internal procedure, used in normalM.
    44384570// Computes the test ideal in the new ring.
     
    44614593  int isGlobal = ord_test(origEre);      // Checks if the original ring has
    44624594                                         // global ordering.
    4463   if(isGlobal != 1){
     4595  if(isGlobal != 1)
     4596  {
    44644597    list rl = ringlist(origEre);
    44654598    list origOrd = rl[3];
     
    44704603    setring ere;
    44714604    ideal norid = imap(origEre, norid);
    4472   } else {
     4605  }
     4606  else
     4607  {
    44734608    def ere = origEre;
    44744609  }
     
    44884623  J = groebner(J);
    44894624  //J = interred(J);
    4490   if(dbg > 1){
     4625  if(dbg > 1)
     4626  {
    44914627    "The radical in the generated ring is";
    44924628    J;
     
    45174653  // ---------------------- step 1 of the mapping ---------------------------
    45184654  intvec degs;
    4519   for(i = 1; i<=ncols(J); i++){
     4655  for(i = 1; i<=ncols(J); i++)
     4656  {
    45204657    degs[i] = degSubring(J[i], @v);
    45214658  }
    4522   if(dbg > 1){
     4659  if(dbg > 1)
     4660  {
    45234661    "The degrees with respect to the new variables are";
    45244662    degs;
     
    45324670  // ---------------------- step 3 of the mapping ---------------------------
    45334671  ideal z;                    // The variables of the original ring in order.
    4534   for(i = 1; i<=nvars(R); i++){
     4672  for(i = 1; i<=nvars(R); i++)
     4673  {
    45354674    z[i] = var(i);
    45364675  }
    45374676
    45384677  map f = ere, U[2..ncols(U)], z[1..ncols(z)]; // The map to the original ring.
    4539   if(dbg > 1){
     4678  if(dbg > 1)
     4679  {
    45404680    "The map is ";
    45414681    f;
     
    45434683  }
    45444684
    4545   if(dbg > 1){
     4685  if(dbg > 1)
     4686  {
    45464687    "Computing the map...";
    45474688  }
    45484689
    45494690  J = f(mapJ);
    4550   if(dbg > 1){
     4691  if(dbg > 1)
     4692  {
    45514693    "The ideal J mapped back (before lifting) is";
    45524694    J;
     
    45584700  ideal J = imap(R, J);
    45594701  poly c = imap(R, c);
    4560   for(i = 1; i<=ncols(J); i++){
    4561     if(degs[i]>1){
     4702  for(i = 1; i<=ncols(J); i++)
     4703  {
     4704    if(degs[i]>1)
     4705    {
    45624706      J[i] = lift(c^(degs[i]-1), J[i])[1,1];
    4563     } else {
    4564       if(degs[i]==0){
     4707    }
     4708    else
     4709    {
     4710      if(degs[i]==0)
     4711      {
    45654712        J[i] = c*J[i];
    45664713      }
     
    45684715  }
    45694716
    4570   if(dbg > 1){
     4717  if(dbg > 1)
     4718  {
    45714719    "The ideal J lifted is";
    45724720    J;
     
    45854733///////////////////////////////////////////////////////////////////////////////
    45864734
    4587 static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I){
     4735static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I)
     4736{
    45884737// Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the
    45894738// ring is 1/c2 * U2.
     
    46034752///////////////////////////////////////////////////////////////////////////////
    46044753
    4605 static proc changeDenominatorQ(ideal U1, poly c1, poly c2){
     4754static proc changeDenominatorQ(ideal U1, poly c1, poly c2)
     4755{
    46064756// Given a ring in the form 1/c1 * U, it computes a new U2 st the ring
    46074757// is 1/c2 * U2.
     
    46194769///////////////////////////////////////////////////////////////////////////////
    46204770
    4621 static proc checkInclusions(ideal U1, ideal U2){
     4771static proc checkInclusions(ideal U1, ideal U2)
     4772{
    46224773// Checks if the identity A = Hom(J, J) of Grauert-Remmert criterion is
    46234774// satisfied.
     
    46364787  // The following check should always be satisfied.
    46374788  // This is only used for debugging.
    4638   if(dbg > 1){
     4789  if(dbg > 1)
     4790  {
    46394791    "and the inclusion A c Hom(J, J): (this should always be satisfied)";
    46404792    // This interred is used only because a bug in groebner!
     
    46454797      "Something went wrong... (this inclusion should always be satisfied)";
    46464798      ~;
    4647     } else {
     4799    }
     4800    else
     4801    {
    46484802      if(dbg>4){pause();}
    46494803    }
    46504804  }
    46514805
    4652   if(size(reduction1[1]) == 0){
     4806  if(size(reduction1[1]) == 0)
     4807  {
    46534808    // We are done! The ring computed in the last step was normal.
    46544809    return(1);
    4655   } else {
     4810  }
     4811  else
     4812  {
    46564813    return(0);
    46574814  }
     
    46604817///////////////////////////////////////////////////////////////////////////////
    46614818
    4662 static proc degSubring(poly p, intvec @v){
     4819static proc degSubring(poly p, intvec @v)
     4820{
    46634821// Computes the degree of a polynomial taking only some variables as variables
    46644822// and the others as parameters.
     
    46684826  int d = 0;  // The degree
    46694827  int e;      // Degree (auxiliar variable)
    4670   for(i = 1; i <= size(p); i++){
     4828  for(i = 1; i <= size(p); i++)
     4829  {
    46714830    e = sum(leadexp(p[i]), @v);
    46724831    if(e > d){d = e;}
     
    46774836///////////////////////////////////////////////////////////////////////////////
    46784837
    4679 static proc mapBackIdeal(ideal I, poly c, intvec @v){
     4838static proc mapBackIdeal(ideal I, poly c, intvec @v)
     4839{
    46804840// Modifies all polynomials in I so that a map x(i) -> y(i)/c can be
    46814841// carried out.
     
    46844844
    46854845  int i;  // counter
    4686   for(i = 1; i <= ncols(I); i++){
     4846  for(i = 1; i <= ncols(I); i++)
     4847  {
    46874848    I[i] = mapBackPoly(I[i], c, @v);
    46884849  }
     
    46924853///////////////////////////////////////////////////////////////////////////////
    46934854
    4694 static proc mapBackPoly(poly p, poly c, intvec @v){
     4855static proc mapBackPoly(poly p, poly c, intvec @v)
     4856{
    46954857// Multiplies each monomial of p by a power of c so that a map x(i) -> y(i)/c
    46964858// can be carried out.
     
    47014863  int d = degSubring(p, @v);
    47024864  poly g = 0;
    4703   for(i = 1; i <= size(p); i++){
     4865  for(i = 1; i <= size(p); i++)
     4866  {
    47044867    e = sum(leadexp(p[i]), @v);
    47054868    g = g + p[i] * c^(d-e);
     
    55675730   return(result);
    55685731}
    5569 
    55705732example
    55715733{ "EXAMPLE:";
     
    56175779  int n=size(L);
    56185780  for (int i=1;i<=n;i++)
     5781  {
     5782    if (defined(R(i)))
    56195783    {
    5620       if (defined(R(i))) {
    5621         string s="Fixed name R("+string(i)+") leads to conflict with existing "
    5622               +"object having this name";
    5623         ERROR(s);
    5624       }
    5625       def R(i)=L[i];
    5626       export R(i);
    5627     }
    5628 
     5784      string s="Fixed name R("+string(i)+") leads to conflict with existing "
     5785            +"object having this name";
     5786      ERROR(s);
     5787    }
     5788    def R(i)=L[i];
     5789    export R(i);
     5790  }
    56295791  return();
    56305792}
     
    56625824//               ("primCl" option must also be set).
    56635825
    5664 proc timeNormal(ideal I, list #){
     5826proc timeNormal(ideal I, list #)
     5827{
     5828  def r = basering;
     5829
    56655830  //--------------------------- define the method ---------------------------
    56665831  int isPrim, useRing;
     
    57005865
    57015866  int tt;
    5702   if(norM){
     5867  if(norM)
     5868  {
    57035869    tt = timer;
    5704     if(decomp == 0){
     5870    if(decomp == 0)
     5871    {
    57055872      "Running normal(useRing, isPrim)...";
    57065873      list a1 = normal(I, "useRing", "isPrim");
    57075874      "Time normal(useRing, isPrim): ", timer - tt;
    5708     } else {
     5875    }
     5876    else
     5877    {
    57095878      "Running normal(useRing)...";
    57105879      list a1 = normal(I, "useRing");
     
    57135882    "";
    57145883  }
    5715   if(norP){
     5884  if(norP)
     5885  {
    57165886    tt = timer;
    5717     if(decomp == 0){
     5887    if(decomp == 0)
     5888    {
    57185889      "Running normalP(isPrim)...";
    57195890      list a2 = normalP(I, "isPrim");
    57205891      "Time normalP(isPrim): ", timer - tt;
    5721     } else {
     5892    }
     5893    else
     5894    {
    57225895      "Running normalP()...";
    57235896      list a2 = normalP(I);
     
    57275900  }
    57285901
    5729   if(norC){
     5902  if(norC)
     5903  {
    57305904    tt = timer;
    5731     if(decomp == 0){
     5905    if(decomp == 0)
     5906    {
    57325907      "Running normalC(isPrim)...";
    57335908      list a3 = normalC(I, "isPrim");
    57345909      "Time normalC(isPrim): ", timer - tt;
    5735     } else {
     5910    }
     5911    else
     5912    {
    57365913      "Running normalC()...";
    57375914      list a3 = normalC(I);
     
    57415918  }
    57425919
    5743   if(primCl){
     5920  if(primCl)
     5921  {
    57445922    tt = timer;
    5745     if(decomp == 0){
     5923    if(decomp == 0)
     5924    {
    57465925      "Running normalC(withGens, isPrim)...";
    57475926      list a4 = normalC(I, "isPrim", "withGens");
    57485927      "Time normalC(withGens, isPrim): ", timer - tt;
    5749     } else {
     5928    }
     5929    else
     5930    {
    57505931      "Running normalC(withGens)...";
    57515932      list a4 = normalC(I, "withGens");
     
    57555936  }
    57565937
    5757   if(check111 and norM){
     5938  if(check111 and norM)
     5939  {
    57585940    "Checking output with norTest...";
    57595941    "WARNING: this checking only works if the original ideal was prime.";
     
    57625944  }
    57635945
    5764   if(checkP and norP and norM){
     5946  if(checkP and norP and norM)
     5947  {
    57655948    "Comparing with normalP output...";
    5766     if(size(a2) > 0){
     5949    if(size(a2) > 0)
     5950    {
    57675951      "WARNING: this checking only works if the original ideal was prime.";
    57685952      U1 = a1[2][1];
     
    57755959      ideal W = fetch(r, W);
    57765960      ch = 0;
    5777       if(size(reduce(U2, groebner(W))) == 0){
     5961      if(size(reduce(U2, groebner(W))) == 0)
     5962      {
    57785963        "U2 c U1";
    57795964        ch = 1;
    57805965      }
    5781       if(size(reduce(W, groebner(U2))) == 0){
     5966      if(size(reduce(W, groebner(U2))) == 0)
     5967      {
    57825968        "U1 c U2";
    57835969        ch = ch + 1;
    57845970      }
    5785       if(ch == 2){
     5971      if(ch == 2)
     5972      {
    57865973        "Output of normalP is equal.";
    5787       } else {
     5974      }
     5975      else
     5976      {
    57885977        "ERROR: Output of normalP is different.";
    57895978      }
    57905979      setring r;
    57915980      kill q;
    5792     } else {
     5981    }
     5982    else
     5983    {
    57935984      "normalP returned no output. Comparison is not possible.";
    57945985    }
     
    57965987  }
    57975988
    5798   if(checkPC and norM and primCl){
     5989  if(checkPC and norM and primCl)
     5990  {
    57995991    "Comparing with primeClosure output...";
    5800     if(size(a4) > 0){
     5992    if(size(a4) > 0)
     5993    {
    58015994      "WARNING: this checking only works if the original ideal was prime.";
    58025995      // primeClosure check
     
    58106003      ideal W = fetch(r, W);
    58116004      ch = 0;
    5812       if(size(reduce(U2, groebner(W))) == 0){
     6005      if(size(reduce(U2, groebner(W))) == 0)
     6006      {
    58136007        "U2 c U1";
    58146008        ch = 1;
    58156009      }
    5816       if(size(reduce(W, groebner(U2))) == 0){
     6010      if(size(reduce(W, groebner(U2))) == 0)
     6011      {
    58176012        "U1 c U2";
    58186013        ch = ch + 1;
    58196014      }
    5820       if(ch == 2){
     6015      if(ch == 2)
     6016      {
    58216017        "Output of normalC(withGens) is equal.";
    5822       } else {
     6018      }
     6019      else
     6020      {
    58236021        "ERROR: Output of normalC(withGens) is different.";
    58246022      }
    58256023      setring r;
    58266024      kill q;
    5827     } else {
     6025    }
     6026    else
     6027    {
    58286028      "normalC(withGens) returned no output. Comparison is not possible.";
    58296029    }
    58306030    "";
    58316031  }
     6032}
     6033
     6034///////////////////////////////////////////////////////////////////////////
     6035static proc sqroot(int n);
     6036{
     6037  int s = 1;
     6038  while(s*s < n)
     6039  {
     6040    s++;
     6041  }
     6042  return(s);
    58326043}
    58336044
Note: See TracChangeset for help on using the changeset viewer.