Changeset 0266ac in git


Ignore:
Timestamp:
Jun 20, 2006, 5:59:42 PM (17 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
Children:
e6cbed94a014891c3db90b82149bb853254bb0d9
Parents:
88b2dd92d922dd2e8a9f9759e7edde2f1e1d2829
Message:
*hannes: bug fixes for 1.117


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/primdec.lib

    r88b2dd r0266ac  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: primdec.lib,v 1.119 2006-06-20 11:55:10 Singular Exp $";
     2version="$Id: primdec.lib,v 1.120 2006-06-20 15:59:42 Singular Exp $";
    33category="Commutative Algebra";
    44info="
    … …  
    10221022  }
    10231023
    1024   if(voice>=8){ primary=extF(primary) }
     1024  if(voice>=8){primary=extF(primary);};
    10251025
    10261026//test whether all ideals in the decomposition are primary and
    … …  
    20642064   if(size(#) > 0){
    20652065     int valid;
    2066      for(j = 1; j <= size(#); j++)
    2067      {
     2066     for(j = 1; j <= size(#); j++){
    20682067       valid = 0;
    2069        if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
    2070        {
    2071          if (#[j] == 0) {facstdOption = "noFacstd"; valid = 1;} // facstd is not used.
    2072          if (#[j] == 1) {facstdOption = "facstd";   valid = 1;} // facstd is used.
    2073        }
    2074        if(typeof(#[j]) == "string")
    2075        {
    2076          if(#[j] == "GTZ" || #[j] == "SL")
    2077          {
     2068       if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) {
     2069         if (#[j] == 0) {facstdOption = "noFacstd"; valid = 1;};    // If #[j] == 0, facstd is not used.
     2070         if (#[j] == 1) {facstdOption = "facstd";   valid = 1;};    // If #[j] == 1, facstd is used.
     2071       };
     2072       if(typeof(#[j]) == "string"){
     2073         if(#[j] == "GTZ" || #[j] == "SL") {
    20782074           algorithm = #[j];
    20792075           valid = 1;
    2080          }
    2081          if(#[j] == "noFacstd" || #[j] == "facstd")
    2082          {
     2076         };
     2077         if(#[j] == "noFacstd" || #[j] == "facstd") {
    20832078           facstdOption = #[j];
    20842079           valid = 1;
    2085          }
    2086        }
    2087        if(valid == 0)
    2088        {
     2080         };
     2081       };
     2082       if(valid == 0) {
    20892083         dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
    2090        }
    2091      }
    2092    }
     2084       };
     2085     };
     2086   };
    20932087
    20942088   list q = simplifyIdeal(i);
    … …  
    21422136            attrib(pr[k], "isSB", 1);
    21432137            // Lines changed
    2144             if (algorithm == "GTZ")
    2145             {
     2138            if (algorithm == "GTZ") {
    21462139              qr = decomp(pr[k], 2);
    2147             }
    2148             else
    2149             {
     2140          } else {
    21502141              qr = newMinAssSL(pr[k]);
    2151             }
     2142          };
    21522143            for(j = 1; j <= size(qr) / 2; j++)
    21532144            {
    … …  
    21632154
    21642155   // Lines changed
    2165    if ((facstdOption == "noFacstd") || (dim(i) == 0))
    2166    {
    2167      if (algorithm == "GTZ")
    2168      {
    2169         re[1] = decomp(i, 2);
    2170      }
    2171      else
    2172      {
    2173        re[1] = newMinAssSL(i);
    2174      }
    2175      re = union(re);
    2176      option(set, op);
    2177      return(phi(re));
    2178    }
     2156   if ((facstdOption == "noFacstd") || (dim(i) == 0)) {
     2157      if (algorithm == "GTZ") {
     2158         re[1] = decomp(i, 2);
     2159      } else {
     2160         re[1] = newMinAssSL(i);
     2161    };
     2162      re = union(re);
     2163      option(set, op);
     2164      return(phi(re));
     2165   };
    21792166
    21802167   q = facstd(i);
    … …  
    21962183      }
    21972184      kill gnir;
    2198    }
     2185   };
    21992186*/
    22002187
    … …  
    22082195      dbprint(printlevel - voice, "We compute the decomp of component", j);
    22092196      // Lines changed
    2210       if (algorithm == "GTZ")
    2211       {
    2212         re[j] = decomp(q[j], 2);
    2213       }
    2214       else
    2215       {
    2216         re[j] = newMinAssSL(q[j]);
    2217       }
     2197      if (algorithm == "GTZ") {
     2198         re[j] = decomp(q[j], 2);
     2199      } else {
     2200         re[j] = newMinAssSL(q[j]);
     2201    };
    22182202      // Debug
    22192203      dbprint(printlevel - voice, "Number of components obtained for this component:", size(re[j]) / 2);
    22202204      dbprint(printlevel - voice, "re[j]:", re[j]);
    2221    }
     2205   };
    22222206   re = union(re);
    22232207
    … …  
    55185502   if(size(#) > 0){
    55195503     int valid;
    5520      for(j = 1; j <= size(#); j++)
    5521      {
     5504     for(j = 1; j <= size(#); j++){
    55225505       valid = 0;
    5523        if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
    5524        {
    5525          if (#[j] == 1) {facstdOption = "noFacstd"; valid = 1;} // facstd is not used.
    5526          if (#[j] == 0) {facstdOption = "facstd";   valid = 1;} // facstd is used.
    5527        }
    5528        if(typeof(#[j]) == "string")
    5529        {
    5530          if((#[j] == "GTZ") || (#[j] == "SL"))
    5531          {
     5506       if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) {
     5507         if (#[j] == 1) {facstdOption = "noFacstd"; valid = 1;};    // If #[j] == 1, facstd is not used.
     5508         if (#[j] == 0) {facstdOption = "facstd";   valid = 1;};    // If #[j] == 0, facstd is used.
     5509       };
     5510       if(typeof(#[j]) == "string"){
     5511         if((#[j] == "GTZ") || (#[j] == "SL")) {
    55325512           algorithm = #[j];
    55335513           valid = 1;
    5534          }
    5535          if((#[j] == "noFacstd") || (#[j] == "facstd"))
    5536          {
     5514         };
     5515         if((#[j] == "noFacstd") || (#[j] == "facstd")) {
    55375516           facstdOption = #[j];
    55385517           valid = 1;
    5539          }
    5540        }
    5541        if(valid == 0)
    5542        {
     5518         };
     5519       };
     5520       if(valid == 0) {
    55435521         dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
    5544        }
    5545      }
    5546    }
     5522       };
     5523     };
     5524   };
    55475525
    55485526   if(ord_test(basering)!=1)
    … …  
    55555533   {
    55565534      return(algeDeco(i,2));
    5557    }
     5535   };
    55585536
    55595537   list result;
    … …  
    56615639
    56625640   // Set input parameters
    5663    algorithm = "SL";     // Default: SL algorithm
    5664    il = 0;               // Default: Full radical (not only equidim part)
    5665    if (homog(i) == 1)
    5666    {                     // Default: facStd is used, except if the ideal is homogeneous.
    5667      useFac = 0;
    5668    }
    5669    else
    5670    {
    5671      useFac = 1;
    5672    }
    5673    if(size(#) > 0)
    5674    {
    5675      int valid;
    5676      for(j = 1; j <= size(#); j++)
    5677      {
    5678        valid = 0;
    5679        if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
    5680        {
    5681          il = #[j];        // If il == 1, equiRadical is computed
    5682          valid = 1;
    5683        }
    5684        if(typeof(#[j]) == "string")
    5685        {
    5686          if(#[j] == "KL")
    5687          {
    5688            algorithm = "KL";
    5689            valid = 1
    5690          }
    5691          if(#[j] == "KLdp")
    5692          {
    5693            algorithm = "KLdp";
    5694            valid = 1
    5695          }
    5696          if(#[j] == "SL")
    5697          {
    5698            algorithm = "SL";
    5699            valid = 1
    5700          }
    5701          if(#[j] == "noFacstd")
    5702          {
    5703            useFac = 0;
    5704            valid = 1
    5705          }
    5706          if(#[j] == "facstd")
    5707          {
    5708            useFac = 1;
    5709            valid = 1
    5710          }
    5711          if(#[j] == "equiRad")
    5712          {
    5713            il = 1;
    5714            valid = 1
    5715          }
    5716          if(#[j] == "fullRad")
    5717          {
    5718            il = 0;
    5719            valid = 1
    5720          }
    5721        }
    5722        if(valid == 0)
    5723        {
    5724          dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
    5725        }
    5726      }
    5727    }
     5641   algorithm = "SL";                                 // Default: SL algorithm
     5642   il = 0;                                                        // Default: Full radical (not only equidim part)
     5643   if (homog(i) == 1) {                                // Default: facStd is used, except if the ideal is homogeneous.
     5644      useFac = 0;
     5645   } else {
     5646             useFac = 1;
     5647   };
     5648   if(size(#) > 0){
     5649           int valid;
     5650           for(j = 1; j <= size(#); j++){
     5651                   valid = 0;
     5652                   if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) {
     5653                           il = #[j];                        // If il == 1, equiRadical is computed
     5654                           valid = 1;
     5655                   };
     5656                   if(typeof(#[j]) == "string"){
     5657                           if(#[j] == "KL") {
     5658                                   algorithm = "KL";
     5659                                   valid = 1};
     5660                           if(#[j] == "KLdp") {
     5661                                   algorithm = "KLdp";
     5662                                   valid = 1};
     5663                           if(#[j] == "SL") {
     5664                                   algorithm = "SL";
     5665                                   valid = 1};
     5666                           if(#[j] == "noFacstd") {
     5667                                   useFac = 0;
     5668                                   valid = 1};
     5669                           if(#[j] == "facstd") {
     5670                                   useFac = 1;
     5671                                      valid = 1};
     5672                           if(#[j] == "equiRad") {
     5673                                   il = 1;
     5674                                   valid = 1};
     5675                           if(#[j] == "fullRad") {
     5676                                   il = 0;
     5677                                   valid = 1};
     5678                   };
     5679                   if(valid == 0) {
     5680                           dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
     5681                   };
     5682           };
     5683   };
    57285684
    57295685   if(size(i) == 0){return(ideal(0));}
    … …  
    57515707   {
    57525708     pr = facstd(i);
    5753    }
    5754    else
    5755    {
    5756      pr = i
    5757    }
     5709   } else {
     5710         pr = i
     5711   };
    57585712   option(set, op);
    57595713   int s = size(pr);
    5760    if(useFac == 1)
    5761    {
     5714   if(useFac == 1) {
    57625715      dbprint(printlevel - voice, "Number of components returned by facstd: ", s);
    5763    }
     5716   };
    57645717   for(j = 1; j <= s; j++)
    57655718   {
    5766      attrib(pr[s + 1 - j], "isSB", 1);
    5767      if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il))
    5768      {
    5769        // SL Debug messages
    5770        dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
    5771        dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));
    5772 
    5773        if(algorithm == "KL")
    5774        {
    5775          rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
    5776        }
    5777        if(algorithm == "KLdp")
    5778        {
    5779          rad = intersect(rad, newRadicalKL(pr[s + 1 - j], rad, il));
    5780        }
    5781        if(algorithm == "SL")
    5782        {
    5783          rad = intersect(rad, newRadicalSL(pr[s + 1 - j], il));
    5784        }
    5785      }
    5786      else
    5787      {
    5788        // SL Debug
    5789        dbprint(printlevel-voice, "The radical of this component is not needed.");
    5790        dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 1)));
    5791        dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
    5792        dbprint(printlevel-voice, "il", il);
    5793      }
     5719      attrib(pr[s + 1 - j], "isSB", 1);
     5720      if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il))
     5721      {
     5722         // SL Debug messages
     5723         dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
     5724         dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));
     5725
     5726             if(algorithm == "KL") {
     5727                  rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
     5728              }
     5729             if(algorithm == "KLdp") {
     5730                  rad = intersect(rad, newRadicalKL(pr[s + 1 - j], rad, il));
     5731              }
     5732              if(algorithm == "SL") {
     5733                      rad = intersect(rad, newRadicalSL(pr[s + 1 - j], il));
     5734              };
     5735      } else
     5736      {
     5737              // SL Debug
     5738          dbprint(printlevel-voice, "The radical of this component is not needed.");
     5739          dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 1)));
     5740          dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
     5741          dbprint(printlevel-voice, "il", il);
     5742      };
    57945743   }
    57955744   return(interred(phi(rad)));
    … …  
    58165765// and radicalKL removed.
    58175766//
    5818 static proc newRadicalKL(ideal I, ideal ser, list #)
    5819 {
    5820 // I             The ideal for which the radical is computed
    5821 // ideal ser     Same as in radicalKL (used to reduce components already obtained)
    5822 // list #        If #[1] = 1, equiradical is computed.
    5823 
    5824   // I needs to be a Groebner basis.
    5825   if (attrib(I, "isSB") != 1)
    5826   {
    5827     I = groebner(I);
    5828   }
    5829 
    5830   ideal rad;                       // The radical
    5831   int allIndep = 1;                // All max independent sets are used
    5832 
    5833   list result = newRadicalReduction(I, ser, allIndep, #);
    5834   int done = result[3];
    5835   rad = result[1];
    5836   if (done == 0)
    5837   {
    5838     rad = intersect(rad, newRadicalKL(result[2], ideal(1), #));
    5839   }
    5840   return(rad);
    5841 }
     5767static proc newRadicalKL(ideal I, ideal ser, list #) {
     5768// I                                The ideal for which the radical is computed
     5769// ideal ser                Same as in radicalKL (used to reduce components already obtained)
     5770// list #                        If #[1] = 1, equiradical is computed.
     5771
     5772        // I needs to be a Groebner basis.
     5773        if (attrib(I, "isSB") != 1) {
     5774                I = groebner(I);
     5775        };
     5776
     5777        ideal rad;                                // The radical
     5778        int allIndep = 1;                // All max independent sets are used
     5779
     5780        list result = newRadicalReduction(I, ser, allIndep, #);
     5781        int done = result[3];
     5782        rad = result[1];
     5783        if (done == 0) {
     5784                rad = intersect(rad, newRadicalKL(result[2], ideal(1), #));
     5785        };
     5786        return(rad);
     5787};
     5788
    58425789
    58435790///////////////////////////////////////////////////////////////////////////////
    … …  
    58545801// obtained in intermediate steps.
    58555802{
    5856   ideal rad = 1;
    5857   ideal equiRad = 1;
    5858   list primes;
    5859   int k;          // Counter
    5860   int il;         // If il = 1, only the equiradical is required.
    5861   int iDim;       // The dimension of I
    5862   int stop = 0;   // Checks if the radical has been obtained
    5863 
    5864   if (attrib(I, "isSB") != 1)
    5865   {
    5866     I = groebner(I);
    5867   }
    5868   iDim = dim(I);
    5869 
    5870   // Checks if only equiradical is required
    5871   if (size(#) > 0)
    5872   {
    5873     il = #[1];
    5874   }
    5875 
    5876   while(stop == 0)
    5877   {
    5878     dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals.");
    5879     primes = newRadicalSLIteration(I, rad);                         // A list of primes or intersections of primes, not included in P
    5880     dbprint (printlevel - voice, "// Output of Iteration Step:");
    5881     dbprint (printlevel - voice, primes);
    5882     if (size(primes) > 0)
    5883     {
    5884       dbprint (printlevel - voice, "// We intersect P with the ideal just obtained.");
    5885       for(k = 1; k <= size(primes); k++)
    5886       {
    5887         rad = intersect(rad, primes[k]);
    5888         if (il == 1)
    5889         {
    5890           if (attrib(primes[k], "isSB") != 1)
    5891           {
    5892             primes[k] = groebner(primes[k]);
    5893           }
    5894           if (iDim == dim(primes[k]))
    5895           {
    5896             equiRad = intersect(equiRad, primes[k]);
    5897           }
    5898         }
    5899       }
    5900     }
    5901     else
    5902     {
    5903       stop = 1;
    5904     }
    5905   }
    5906   if (il == 0)
    5907   {
    5908     return(rad);
    5909   }
    5910   else
    5911   {
    5912     return(equiRad);
    5913   }
    5914 }
     5803        ideal rad = 1;
     5804        ideal equiRad = 1;
     5805        list primes;
     5806        int k;                        // Counter
     5807        int il;                 // If il = 1, only the equiradical is required.
     5808        int iDim;                // The dimension of I
     5809        int stop = 0;   // Checks if the radical has been obtained
     5810
     5811        if (attrib(I, "isSB") != 1) {
     5812                I = groebner(I);
     5813        };
     5814        iDim = dim(I);
     5815
     5816    // Checks if only equiradical is required
     5817    if (size(#) > 0) {
     5818       il = #[1];
     5819    };
     5820
     5821        while(stop == 0) {
     5822                dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals.");
     5823                primes = newRadicalSLIteration(I, rad);                         // A list of primes or intersections of primes, not included in P
     5824                dbprint (printlevel - voice, "// Output of Iteration Step:");
     5825                dbprint (printlevel - voice, primes);
     5826                if (size(primes) > 0) {
     5827                        dbprint (printlevel - voice, "// We intersect P with the ideal just obtained.");
     5828                        for(k = 1; k <= size(primes); k++) {
     5829                                rad = intersect(rad, primes[k]);
     5830                                if (il == 1){
     5831                                        if (attrib(primes[k], "isSB") != 1) {
     5832                                                primes[k] = groebner(primes[k]);
     5833                                        };
     5834                                        if (iDim == dim(primes[k])) {
     5835                                                equiRad = intersect(equiRad, primes[k]);
     5836                                        };
     5837                                };
     5838                        };
     5839                } else {
     5840                        stop = 1;
     5841                };
     5842        };
     5843        if (il == 0) {
     5844                return(rad);
     5845        } else {
     5846                return(equiRad);
     5847        };
     5848};
    59155849
    59165850//////////////////////////////////////////////////////////////////////////
    … …  
    59275861//       This is not used in the new algorithm. It is part of KL algorithm
    59285862
    5929 static proc newRadicalReduction(ideal I, ideal ser, int allIndep, list #)
    5930 {
     5863static proc newRadicalReduction(ideal I, ideal ser, int allIndep, list #) {
    59315864// allMaximal                1 -> Indicates that the reduction to the zerodim case
    59325865//                       must be done for all indep set of the leading terms ideal
    … …  
    59365869//                  It is used to set the value of done.)
    59375870
    5938   attrib(I, "isSB", 1);   // I needs to be a reduced standard basis
    5939   list indep, fett;
    5940   intvec @w, @hilb, op;
    5941   int @wr, @n, @m, lauf, di;
    5942   ideal fac, @h, collectrad, lsau;
    5943   poly @q;
    5944   string @va, quotring;
    5945 
    5946   def @P = basering;
    5947   int jdim = dim(I);               // Computes the dimension of I
    5948   int  homo = homog(I);            // Finds out if I is homogeneous
    5949   ideal rad = ideal(1);            // The unit ideal
    5950   ideal te = ser;
    5951   if(size(#) > 0)
    5952   {
    5953     @wr = #[1];
    5954   }
    5955   if(homo == 1)
    5956   {
    5957     for(@n = 1; @n <= nvars(basering); @n++)
    5958     {
    5959       @w[@n] = ord(var(@n));
    5960     }
    5961     @hilb = hilb(I, 1, @w);
    5962   }
    5963 
    5964   // SL 2006.04.11 1 Debug messages
    5965   dbprint(printlevel-voice, "//Computes the radical of the ideal:", I);
    5966   // SL 2006.04.11 2 Debug messages
     5871   attrib(I, "isSB", 1);   // I needs to be a reduced standard basis
     5872   list indep, fett;
     5873   intvec @w, @hilb, op;
     5874   int @wr, @n, @m, lauf, di;
     5875   ideal fac, @h, collectrad, lsau;
     5876   poly @q;
     5877   string @va, quotring;
     5878
     5879   def @P = basering;
     5880   int jdim = dim(I);               // Computes the dimension of I
     5881   int  homo = homog(I);            // Finds out if I is homogeneous
     5882   ideal rad = ideal(1);            // The unit ideal
     5883   ideal te = ser;
     5884   if(size(#) > 0)
     5885   {
     5886      @wr = #[1];
     5887   }
     5888   if(homo == 1)
     5889   {
     5890      for(@n = 1; @n <= nvars(basering); @n++)
     5891      {
     5892         @w[@n] = ord(var(@n));
     5893      }
     5894      @hilb = hilb(I, 1, @w);
     5895   }
     5896
     5897   // SL 2006.04.11 1 Debug messages
     5898   dbprint(printlevel-voice, "//Computes the radical of the ideal:", I);
     5899   // SL 2006.04.11 2 Debug messages
    59675900
    59685901  //---------------------------------------------------------------------------
    … …  
    59705903  //---------------------------------------------------------------------------
    59715904
    5972   if (jdim==-1)
    5973   {
    5974     return(ideal(1), ideal(1), 1);
    5975   }
     5905   if (jdim==-1)
     5906   {
     5907      return(ideal(1), ideal(1), 1);
     5908   }
    59765909
    59775910  //---------------------------------------------------------------------------
    … …  
    59795912  //---------------------------------------------------------------------------
    59805913
    5981   if (jdim==0)
    5982   {
    5983     return(zeroRad(I), ideal(1), 1);
    5984   }
    5985 
    5986   //-------------------------------------------------------------------------
    5987   //search for a maximal independent set indep,i.e.
    5988   //look for subring such that the intersection with the ideal is zero
    5989   //j intersected with K[var(indep[3]+1),...,var(nvar)] is zero,
    5990   //indep[1] is the new varstring, indep[2] the string for the block-ordering
    5991   //-------------------------------------------------------------------------
    5992 
    5993   // SL 2006-04-24 1         If allIndep = 0, then it only computes one maximal
    5994   //                                        independent set.
    5995   //                                        This looks better for the new algorithm but not for KL
    5996   //                                        algorithm
    5997   list parameters = allIndep;
    5998   indep = newMaxIndependSet(I, parameters);
    5999   // SL 2006-04-24 2
    6000 
    6001   for(@m = 1; @m <= size(indep); @m++)
    6002   {
    6003     if((indep[@m][1] == varstr(basering)) && (@m == 1))
    6004     //this is the good case, nothing to do, just to have the same notations
    6005     //change the ring
    6006     {
    6007       execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
     5914   if (jdim==0)
     5915   {
     5916      return(zeroRad(I), ideal(1), 1);
     5917   }
     5918
     5919   //-------------------------------------------------------------------------
     5920   //search for a maximal independent set indep,i.e.
     5921   //look for subring such that the intersection with the ideal is zero
     5922   //j intersected with K[var(indep[3]+1),...,var(nvar)] is zero,
     5923   //indep[1] is the new varstring, indep[2] the string for the block-ordering
     5924   //-------------------------------------------------------------------------
     5925
     5926   // SL 2006-04-24 1         If allIndep = 0, then it only computes one maximal
     5927   //                                        independent set.
     5928   //                                        This looks better for the new algorithm but not for KL
     5929   //                                        algorithm
     5930   list parameters = allIndep;
     5931   indep = newMaxIndependSet(I, parameters);
     5932   // SL 2006-04-24 2
     5933
     5934   for(@m = 1; @m <= size(indep); @m++)
     5935   {
     5936      if((indep[@m][1] == varstr(basering)) && (@m == 1))
     5937      //this is the good case, nothing to do, just to have the same notations
     5938      //change the ring
     5939      {
     5940         execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
    60085941                              +ordstr(basering)+");");
    6009       ideal @j = fetch(@P, I);
    6010       attrib(@j, "isSB", 1);
    6011     }
    6012     else
    6013     {
    6014       @va = string(maxideal(1));
    6015 
    6016       execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),("
     5942         ideal @j = fetch(@P, I);
     5943         attrib(@j, "isSB", 1);
     5944      }
     5945      else
     5946      {
     5947         @va = string(maxideal(1));
     5948
     5949         execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),("
    60175950                              + indep[@m][2] + ");");
    6018       execute("map phi = @P," + @va + ";");
    6019       if(homo == 1)
    6020       {
    6021         ideal @j = std(phi(I), @hilb, @w);
     5951         execute("map phi = @P," + @va + ";");
     5952         if(homo == 1)
     5953         {
     5954            ideal @j = std(phi(I), @hilb, @w);
     5955         }
     5956         else
     5957         {
     5958            ideal @j = groebner(phi(I));
     5959         }
     5960      }
     5961      if((deg(@j[1]) == 0) || (dim(@j) < jdim))
     5962      {
     5963         setring @P;
     5964         break;
     5965      }
     5966      for (lauf = 1; lauf <= size(@j); lauf++)
     5967      {
     5968         fett[lauf] = size(@j[lauf]);
     5969      }
     5970      //------------------------------------------------------------------------
     5971      // We have now the following situation:
     5972      // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     5973      // to this quotientring, j is there still a standardbasis, the
     5974      // leading coefficients of the polynomials there (polynomials in
     5975      // K[var(nnp+1),..,var(nva)]) are collected in the list h,
     5976      // we need their LCM, gh, because of the following:
     5977      // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
     5978      // intersected with K[var(1),...,var(nva)] is (j:gh^n)
     5979      // on the other hand j = ((j, gh^n) intersected with (j : gh^n))
     5980
     5981      //------------------------------------------------------------------------
     5982      // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
     5983      // and the map phi:K[var(1),...,var(nva)] ----->
     5984      // K(var(nnpr+1),..,var(nva))[..the rest..]
     5985      //------------------------------------------------------------------------
     5986      quotring = newPrepareQuotientRing(nvars(basering) - indep[@m][3]);
     5987      //------------------------------------------------------------------------
     5988      // We pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     5989      //------------------------------------------------------------------------
     5990
     5991      execute(quotring);
     5992
     5993      // @j considered in the quotientring
     5994      ideal @j = imap(gnir1, @j);
     5995
     5996      kill gnir1;
     5997
     5998      // j is a standardbasis in the quotientring but usually not minimal
     5999      // here it becomes minimal
     6000
     6001      @j = clearSB(@j, fett);
     6002
     6003      // We need later LCM(h[1],...) = gh for saturation
     6004      ideal @h;
     6005      if(deg(@j[1]) > 0)
     6006      {
     6007         for(@n = 1; @n <= size(@j); @n++)
     6008         {
     6009            @h[@n] = leadcoef(@j[@n]);
     6010         }
     6011         op = option(get);
     6012         option(redSB);
     6013         @j = interred(@j);  //to obtain a reduced standardbasis
     6014         attrib(@j, "isSB", 1);
     6015         option(set, op);
     6016
     6017         // SL 1 Debug messages
     6018         dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j)));
     6019         ideal zero_rad = zeroRad(@j);
     6020         dbprint(printlevel - voice, "zero_rad passed");
     6021         // SL 2
    60226022      }
    60236023      else
    60246024      {
    6025         ideal @j = groebner(phi(I));
    6026       }
    6027     }
    6028     if((deg(@j[1]) == 0) || (dim(@j) < jdim))
    6029     {
    6030       setring @P;
    6031       break;
    6032     }
    6033     for (lauf = 1; lauf <= size(@j); lauf++)
    6034     {
    6035       fett[lauf] = size(@j[lauf]);
    6036     }
    6037     //------------------------------------------------------------------------
    6038     // We have now the following situation:
    6039     // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
    6040     // to this quotientring, j is there still a standardbasis, the
    6041     // leading coefficients of the polynomials there (polynomials in
    6042     // K[var(nnp+1),..,var(nva)]) are collected in the list h,
    6043     // we need their LCM, gh, because of the following:
    6044     // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
    6045     // intersected with K[var(1),...,var(nva)] is (j:gh^n)
    6046     // on the other hand j = ((j, gh^n) intersected with (j : gh^n))
    6047 
    6048     //------------------------------------------------------------------------
    6049     // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
    6050     // and the map phi:K[var(1),...,var(nva)] ----->
    6051     // K(var(nnpr+1),..,var(nva))[..the rest..]
    6052     //------------------------------------------------------------------------
    6053     quotring = newPrepareQuotientRing(nvars(basering) - indep[@m][3]);
    6054     //------------------------------------------------------------------------
    6055     // We pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
    6056     //------------------------------------------------------------------------
    6057 
    6058     execute(quotring);
    6059 
    6060     // @j considered in the quotientring
    6061     ideal @j = imap(gnir1, @j);
    6062 
    6063     // j is a standardbasis in the quotientring but usually not minimal
    6064     // here it becomes minimal
    6065 
    6066     @j = clearSB(@j, fett);
    6067 
    6068     // We need later LCM(h[1],...) = gh for saturation
    6069     ideal @h;
    6070     if(deg(@j[1]) > 0)
    6071     {
    6072       for(@n = 1; @n <= size(@j); @n++)
    6073       {
    6074         @h[@n] = leadcoef(@j[@n]);
    6075       }
    6076       op = option(get);
    6077       option(redSB);
    6078       @j = interred(@j);  //to obtain a reduced standardbasis
    6079       attrib(@j, "isSB", 1);
    6080       option(set, op);
    6081 
    6082       // SL 1 Debug messages
    6083       dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j)));
    6084       ideal zero_rad = zeroRad(@j);
    6085       dbprint(printlevel - voice, "zero_rad passed");
    6086       // SL 2
    6087     }
    6088     else
    6089     {
    6090       ideal zero_rad = ideal(1);
    6091     }
    6092 
    6093     // We need the intersection of the ideals in the list quprimary with the
    6094     // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
    6095     // but fi polynomials, then the intersection of q with the polynomialring
    6096     // is the saturation of the ideal generated by f1,...,fr with respect to
    6097     // h which is the lcm of the leading coefficients of the fi considered in
     6025         ideal zero_rad = ideal(1);
     6026      }
     6027
     6028      // We need the intersection of the ideals in the list quprimary with the
     6029      // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     6030      // but fi polynomials, then the intersection of q with the polynomialring
     6031      // is the saturation of the ideal generated by f1,...,fr with respect to
     6032      // h which is the lcm of the leading coefficients of the fi considered in
    60986033      // the quotientring: this is coded in saturn
    60996034
    … …  
    61866121   } else {
    61876122      done = 0;
    6188    }
     6123   };
    61896124   // SL 2006.04.21 2
    61906125
    … …  
    61936128   return(result);
    61946129   // SL 2006.04.21 2
    6195 }
     6130};
    61966131
    61976132
    … …  
    62156150                good = 1 - rad_con(P[k], I);
    62166151                k++;
    6217         }
     6152        };
    62186153        k--;
    62196154        if (good == 0) {
    … …  
    62216156                list emptyList = list();
    62226157                return (emptyList);
    6223         }
     6158        };
    62246159        dbprint(printlevel - voice, "// That one was good!");
    62256160        dbprint(printlevel - voice, "// We saturate I with respect to this element.");
    … …  
    62296164                dbprint(printlevel - voice, "// The polynomial is 1, the saturation in not actually computed.");
    62306165                ideal J = I;
    6231         }
     6166        };
    62326167
    62336168        // We now call proc radicalNew;
    … …  
    62426177
    62436178        return(result[1]);
    6244 }
     6179};
    62456180
    62466181
    … …  
    62786213   } else {
    62796214           allMaximal = 1;
    6280    }
     6215   };
    62816216
    62826217   int nMax;
    … …  
    62856220   } else {
    62866221          nMax = 1;
    6287    }
     6222   };
    62886223
    62896224   for(n = 1; n <= nMax; n++)
    … …  
    66566591       indepOption = #[count];
    66576592        count++;
    6658      }
    6659    }
     6593     };
     6594   };
    66606595   if(typeof(#[count]) == "string") {
    66616596     if ((#[count] == "intersect") or (#[count] == "noIntersect")){
    66626597       intersectOption = #[count];
    66636598        count++;
    6664      }
    6665    }
     6599     };
     6600   };
    66666601     if((typeof(#[count]) == "int") or (typeof(#[count]) == "number"))
    66676602     {
    … …  
    66716606          if(@wr==3){abspri = 1; @wr = 0;}
    66726607          count++;
    6673        }
    6674      }
     6608       };
     6609     };
    66756610     if(size(#)>count)
    66766611     {
    … …  
    66786613          peek=#[count + 1];
    66796614          ser=#[count + 2];
    6680      }
    6681   }
     6615     };
     6616  };
    66826617  if(abspri)
    66836618  {
    … …  
    67146649       } else {
    67156650         return(l);
    6716        }
     6651       };
    67176652        }
    67186653        if (intersectOption == "intersect") {
    … …  
    67206655       } else {
    67216656         return(primary);
    6722        }
     6657       };
    67236658
    67246659     }
    … …  
    67426677    } else {
    67436678       return(primary);
    6744    }
     6679   };
    67456680  }
    67466681
    … …  
    68116746     } else {
    68126747       return(primary);
    6813      }
     6748     };
    68146749    }
    68156750    if(size(fried)>0)
    … …  
    68516786       } else {
    68526787         list pr = result;
    6853        }
     6788       };
    68546789
    68556790       setring gnir;
    … …  
    68596794         @j=pr[@k]+fried;
    68606795         pr[@k]=@j;
    6861        }
     6796       };
    68626797       if (intersectOption == "intersect") {
    68636798      ideal intersection = imap(@deirf, intersection);
    68646799      @j = intersection + fried;
    68656800      intersection = @j;
    6866      }
     6801     };
    68676802     setring @P;
    68686803     if (intersectOption == "intersect") {
    … …  
    68706805     } else {
    68716806       return(imap(gnir,pr));
    6872      }
     6807     };
    68736808    }
    68746809  }
    … …  
    68856820   } else {
    68866821     return(primary);
    6887    }
     6822   };
    68886823  }
    68896824
    … …  
    69136848        if (intersectOption == "intersect") {
    69146849          generator = generator * fac[1][@k];
    6915         }
    6916      }
     6850        };
     6851     };
    69176852   if (intersectOption == "intersect") {
    69186853       gIntersection = generator;
    6919    }
     6854   };
    69206855
    69216856     setring @P;
    … …  
    69236858   if (intersectOption == "intersect") {
    69246859       ideal intersection = fetch(gnir,gIntersection);
    6925      }
     6860     };
    69266861
    69276862//HIER
    … …  
    69406875       for(ab=1;ab<=size(primary);ab++) {
    69416876          intersection = intersect(intersection, primary[ab][2]);
    6942         }
    6943      }
     6877        };
     6878     };
    69446879     if (intersectOption == "intersect") {
    69456880     return(list(primary, intersection));
    69466881     } else {
    69476882     return(primary);
    6948      }
     6883     };
    69496884  }
    69506885
    … …  
    69826917   } else {
    69836918     return(primary);
    6984    }
     6919   };
    69856920
    69866921  }
    … …  
    70697004      absprimary = absprimary + result[2];
    70707005      abskeep = abskeep + result[3];
    7071       }
     7006      };
    70727007    @h = result[5];
    70737008    ser = result[4];
    … …  
    71287063        ser = imap(@Phelp, ser);
    71297064        @j = imap(@Phelp, jwork);
    7130      }
     7065     };
    71317066  }
    71327067
    … …  
    74777412   } else {
    74787413     return(primary);
    7479    }
     7414   };
    74807415}
    74817416example
    … …  
    76357570
    76367571       int zeroMinAss = @wr;
    7637        if (@wr == 2) {zeroMinAss = 1;}
     7572       if (@wr == 2) {zeroMinAss = 1;};
    76387573        list uprimary= newZero_decomp(@j, ser, zeroMinAss);
    76397574
    … …  
    77467681
    77477682     return(quprimary, absprimary, abskeep, ser, @h);
    7748 }
     7683};
    77497684
    77507685
    … …  
    78347769      for (k = 1; k <= size(pd); k++) {
    78357770        P = intersect(P, pd[k]);
    7836       }
     7771      };
    78377772      */
    78387773      P = intersect(P, pd[2]);
    … …  
    78407775      dbprint(printlevel - voice, "// Intersection finished.");
    78417776    } else {
    7842       stop = 1;
    7843     }
    7844   }
     7777      stop = 1
     7778    };
     7779  };
    78457780  //return(insert(primaryDec, P));
    78467781  // Returns only the primary components, not the radical.
    78477782  return(primaryDec);
    7848 }
     7783};
    78497784
    78507785///////////////////////////////////////////////////////////////////////////////
    … …  
    78637798    good = 1 - rad_con(P[k], I);
    78647799    k++;
    7865   }
     7800  };
    78667801  k--;
    78677802  if (good == 0) {
    … …  
    78697804    dbprint (printlevel - voice, "// No element was found, P = sqrt(I).");
    78707805    return (list(primaryDec, ideal(0)));
    7871   }
     7806  };
    78727807  // Debug
    78737808  dbprint (printlevel - voice, "// We found h = ", P[k]);
    … …  
    79207855   } else {
    79217856     indepOption = "allIndep";
    7922    }
     7857   };
    79237858
    79247859   int nMax;
    … …  
    79277862   } else {
    79287863    nMax = 1;
    7929    }
     7864   };
    79307865
    79317866   for(n = 1; n <= nMax; n++)
    … …  
    80337968      if (#[1] == "nest") {
    80347969        nestLevel = #[2];
    8035       }
     7970      };
    80367971      # = list();
    8037     }
    8038   }
     7972    };
     7973  };
    80397974
    80407975  if(vdim(j)==deg(j[1]))
    … …  
    81928127  }
    81938128
    8194   if(nestLevel > 1){primary=extF(primary)}
     8129  if(nestLevel > 1){primary=extF(primary);};
    81958130
    81968131//test whether all ideals in the decomposition are primary and
Note: See TracChangeset for help on using the changeset viewer.