Changeset 808a9f3 in git


Ignore:
Timestamp:
Jun 20, 2006, 1:55:10 PM (17 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
Children:
88b2dd92d922dd2e8a9f9759e7edde2f1e1d2829
Parents:
951a6327b56dcd11021810ab04eba401d1432b98
Message:
*hannes: bug fix to 1.117+format


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/primdec.lib

    r951a63 r808a9f3  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: primdec.lib,v 1.118 2006-06-19 15:08:32 Singular Exp $";
     2version="$Id: primdec.lib,v 1.119 2006-06-20 11:55:10 Singular Exp $";
    33category="Commutative Algebra";
    44info="
    … …  
    77@*        Wolfram Decker, decker@math.uni-sb.de         (SY)
    88@*        Hans Schoenemann, hannes@mathematik.uni-kl.de (SY)
     9@*        Santiago Laplagne, slaplagn@dm.uba.ar        (GTZ)
    910
    1011OVERVIEW:
    … …  
    10211022  }
    10221023
    1023   if(voice>=8){primary=extF(primary)};
     1024  if(voice>=8){ primary=extF(primary) }
    10241025
    10251026//test whether all ideals in the decomposition are primary and
    … …  
    20432044static proc minAssPrimes(ideal i, list #)
    20442045"USAGE:   minAssPrimes(i); i ideal
    2045          minAssPrimes(i,1); i ideal  (to use also the factorizing Groebner)
     2046      Optional parameters in list #: (can be entered in any order)
     2047      0, "facstd"   ->   uses facstd to first decompose the ideal
     2048      1, "noFacstd" ->  does not use facstd (default)
     2049      "SL" ->     the new algorithm is used (default)
     2050      "GTZ" ->     the old algorithm is used
    20462051RETURN:  list = the minimal associated prime ideals of i
    20472052EXAMPLE: example minAssPrimes; shows an example
    20482053"
    20492054{
    2050    def P=basering;
    2051    if(size(i)==0){return(list(ideal(0)));}
    2052    list q=simplifyIdeal(i);
    2053    list re=maxideal(1);
    2054    int j,a,k;
    2055    intvec op=option(get);
    2056    map phi=P,q[2];
    2057 
    2058    if(npars(P)==0){option(redSB);}
     2055   string algorithm;    // Algorithm to be used
     2056   string facstdOption;    // To uses proc facstd
     2057   int j;          // Counter
     2058   def P = basering;
     2059   if(size(i) == 0){return(list(ideal(0)));}
     2060
     2061   // Set input parameters
     2062   algorithm = "SL";         // Default: SL algorithm
     2063   facstdOption = "noFacstd";    // Default: facstd is not used
     2064   if(size(#) > 0){
     2065     int valid;
     2066     for(j = 1; j <= size(#); j++)
     2067     {
     2068       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         {
     2078           algorithm = #[j];
     2079           valid = 1;
     2080         }
     2081         if(#[j] == "noFacstd" || #[j] == "facstd")
     2082         {
     2083           facstdOption = #[j];
     2084           valid = 1;
     2085         }
     2086       }
     2087       if(valid == 0)
     2088       {
     2089         dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
     2090       }
     2091     }
     2092   }
     2093
     2094   list q = simplifyIdeal(i);
     2095   list re = maxideal(1);
     2096   int a, k;
     2097   intvec op = option(get);
     2098   map phi = P,q[2];
     2099
     2100   list result;
     2101
     2102   if(npars(P) == 0){option(redSB);}
    20592103
    20602104   if(attrib(i,"isSB")!=1)
    … …  
    20732117      }
    20742118   }
    2075    if(dim(i)==-1){return(ideal(1));}
    2076    if((dim(i)==0)&&(npars(P)==0))
    2077    {
    2078       int di=vdim(i);
     2119
     2120   if(dim(i) == -1){re = ideal(1); return(re);}
     2121   if((dim(i) == 0) && (npars(P) == 0))
     2122   {
     2123      int di = vdim(i);
    20792124      execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
    2080       ideal J=std(imap(P,i));
    2081       attrib(J,"isSB",1);
    2082       if(vdim(J)!=di)
    2083       {
    2084          J=fglm(P,i);
    2085       }
    2086       list pr=triangMH(J,2);
    2087       list qr,re;
    2088 
    2089       for(k=1;k<=size(pr);k++)
     2125      ideal J = std(imap(P,i));
     2126      attrib(J, "isSB", 1);
     2127      if(vdim(J) != di)
     2128      {
     2129         J = fglm(P, i);
     2130      }
     2131      list pr = triangMH(J, 2);
     2132      list qr, re;
     2133
     2134      for(k = 1; k <= size(pr); k++)
    20902135      {
    20912136         if(primT(pr[k]))
    20922137         {
    2093             re[size(re)+1]=pr[k];
     2138            re[size(re) + 1] = pr[k];
    20942139         }
    20952140         else
    20962141         {
    2097             attrib(pr[k],"isSB",1);
    2098             qr=decomp(pr[k],2);
    2099             for(j=1;j<=size(qr)/2;j++)
     2142            attrib(pr[k], "isSB", 1);
     2143            // Lines changed
     2144            if (algorithm == "GTZ")
    21002145            {
    2101                re[size(re)+1]=qr[2*j];
     2146              qr = decomp(pr[k], 2);
     2147            }
     2148            else
     2149            {
     2150              qr = newMinAssSL(pr[k]);
     2151            }
     2152            for(j = 1; j <= size(qr) / 2; j++)
     2153            {
     2154               re[size(re) + 1] = qr[2 * j];
    21022155            }
    21032156         }
    21042157      }
    21052158      setring P;
    2106       re=imap(gnir,re);
    2107       option(set,op);
     2159      re = imap(gnir, re);
     2160      option(set, op);
    21082161      return(phi(re));
    21092162   }
    21102163
    2111    if((size(#)==0)||(dim(i)==0))
    2112    {
    2113       re[1]=decomp(i,2);
    2114       re=union(re);
    2115       option(set,op);
    2116       return(phi(re));
    2117    }
    2118 
    2119    q=facstd(i);
     2164   // 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   }
     2179
     2180   q = facstd(i);
    21202181
    21212182/*
    2122    if((size(q)==1)&&(dim(i)>1))
     2183   if((size(q) == 1) && (dim(i) > 1))
    21232184   {
    21242185      execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
    2125 
    2126       list p=facstd(fetch(P,i));
    2127       if(size(p)>1)
    2128       {
    2129          a=1;
     2186      list p = facstd(fetch(P, i));
     2187      if(size(p) > 1)
     2188      {
     2189         a = 1;
    21302190         setring P;
    2131          q=fetch(gnir,p);
     2191         q = fetch(gnir,p);
    21322192      }
    21332193      else
    … …  
    21402200
    21412201   option(set,op);
    2142    for(j=1;j<=size(q);j++)
    2143    {
    2144       if(a==0){attrib(q[j],"isSB",1);}
    2145       re[j]=decomp(q[j],2);
    2146    }
    2147    re=union(re);
     2202   // Debug
     2203   dbprint(printlevel - voice, "Components returned by facstd", size(q), q);
     2204   for(j = 1; j <= size(q); j++)
     2205   {
     2206      if(a == 0){attrib(q[j], "isSB", 1);}
     2207      // Debug
     2208      dbprint(printlevel - voice, "We compute the decomp of component", j);
     2209      // 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      }
     2218      // Debug
     2219      dbprint(printlevel - voice, "Number of components obtained for this component:", size(re[j]) / 2);
     2220      dbprint(printlevel - voice, "re[j]:", re[j]);
     2221   }
     2222   re = union(re);
     2223
    21482224   return(phi(re));
    21492225}
    … …  
    54205496///////////////////////////////////////////////////////////////////////////////
    54215497proc minAssGTZ(ideal i,list #)
    5422 "USAGE:   minAssGTZ(i); i ideal
    5423           minAssGTZ(i,1); i ideal  does not use the factorizing Groebner
     5498"USAGE:    minAssGTZ(i); i ideal
     5499      Optional parameters in list #: (can be entered in any order)
     5500      0, "facstd"   ->   uses facstd to first decompose the ideal (default)
     5501      1, "noFacstd" ->  does not use facstd
     5502      "SL" ->     the new algorithm is used (default)
     5503      "GTZ" ->     the old algorithm is used
    54245504RETURN:  a list, the minimal associated prime ideals of i.
    54255505NOTE:    Designed for characteristic 0, works also in char k > 0 based
    … …  
    54285508"
    54295509{
     5510   int j;
     5511   string algorithm;
     5512   string facstdOption;
     5513   int useFac;
     5514
     5515   // Set input parameters
     5516   algorithm = "SL";         // Default: SL algorithm
     5517   facstdOption = "facstd";
     5518   if(size(#) > 0){
     5519     int valid;
     5520     for(j = 1; j <= size(#); j++)
     5521     {
     5522       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         {
     5532           algorithm = #[j];
     5533           valid = 1;
     5534         }
     5535         if((#[j] == "noFacstd") || (#[j] == "facstd"))
     5536         {
     5537           facstdOption = #[j];
     5538           valid = 1;
     5539         }
     5540       }
     5541       if(valid == 0)
     5542       {
     5543         dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
     5544       }
     5545     }
     5546   }
     5547
    54305548   if(ord_test(basering)!=1)
    54315549   {
    … …  
    54385556      return(algeDeco(i,2));
    54395557   }
    5440    if(size(#)==0)
    5441    {
    5442       return(minAssPrimes(i,1));
    5443    }
    5444    return(minAssPrimes(i));
     5558
     5559   list result;
     5560   result = minAssPrimes(i, facstdOption, algorithm);
     5561   return(result);
    54455562}
    54465563example
    … …  
    55445661
    55455662   // Set input parameters
    5546    algorithm = "SL";                                 // Default: SL algorithm
    5547    il = 0;                                                        // Default: Full radical (not only equidim part)
    5548    if (homog(i) == 1) {                                // Default: facStd is used, except if the ideal is homogeneous.
    5549       useFac = 0;
    5550    } else {
    5551              useFac = 1;
    5552    };
    5553    if(size(#) > 0){
    5554            int valid;
    5555            for(j = 1; j <= size(#); j++){
    5556                    valid = 0;
    5557                    if((typeof(#[j]) == "int") or (typeof(#[j]) == "number")) {
    5558                            il = #[j];                        // If il == 1, equiRadical is computed
    5559                            valid = 1;
    5560                    };
    5561                    if(typeof(#[j]) == "string"){
    5562                            if(#[j] == "KL") {
    5563                                    algorithm = "KL";
    5564                                    valid = 1};
    5565                            if(#[j] == "KLdp") {
    5566                                    algorithm = "KLdp";
    5567                                    valid = 1};
    5568                            if(#[j] == "SL") {
    5569                                    algorithm = "SL";
    5570                                    valid = 1};
    5571                            if(#[j] == "noFacstd") {
    5572                                    useFac = 0;
    5573                                    valid = 1};
    5574                            if(#[j] == "facstd") {
    5575                                    useFac = 1;
    5576                                       valid = 1};
    5577                            if(#[j] == "equiRad") {
    5578                                    il = 1;
    5579                                    valid = 1};
    5580                            if(#[j] == "fullRad") {
    5581                                    il = 0;
    5582                                    valid = 1};
    5583                    };
    5584                    if(valid == 0) {
    5585                            dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
    5586                    };
    5587            };
    5588    };
     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   }
    55895728
    55905729   if(size(i) == 0){return(ideal(0));}
    … …  
    56125751   {
    56135752     pr = facstd(i);
    5614    } else {
    5615          pr = i
    5616    };
     5753   }
     5754   else
     5755   {
     5756     pr = i
     5757   }
    56175758   option(set, op);
    56185759   int s = size(pr);
    5619    if(useFac == 1) {
     5760   if(useFac == 1)
     5761   {
    56205762      dbprint(printlevel - voice, "Number of components returned by facstd: ", s);
    5621    };
     5763   }
    56225764   for(j = 1; j <= s; j++)
    56235765   {
    5624       attrib(pr[s + 1 - j], "isSB", 1);
    5625       if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il))
    5626       {
    5627          // SL Debug messages
    5628          dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
    5629          dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));
    5630 
    5631              if(algorithm == "KL") {
    5632                   rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
    5633               }
    5634              if(algorithm == "KLdp") {
    5635                   rad = intersect(rad, newRadicalKL(pr[s + 1 - j], rad, il));
    5636               }
    5637               if(algorithm == "SL") {
    5638                       rad = intersect(rad, newRadicalSL(pr[s + 1 - j], il));
    5639               };
    5640       } else
    5641       {
    5642               // SL Debug
    5643           dbprint(printlevel-voice, "The radical of this component is not needed.");
    5644           dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))", size(reduce(rad, pr[s + 1 - j], 1)));
    5645           dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
    5646           dbprint(printlevel-voice, "il", il);
    5647       };
     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     }
    56485794   }
    56495795   return(interred(phi(rad)));
    … …  
    56705816// and radicalKL removed.
    56715817//
    5672 static proc newRadicalKL(ideal I, ideal ser, list #) {
    5673 // I                                The ideal for which the radical is computed
    5674 // ideal ser                Same as in radicalKL (used to reduce components already obtained)
    5675 // list #                        If #[1] = 1, equiradical is computed.
    5676 
    5677         // I needs to be a Groebner basis.
    5678         if (attrib(I, "isSB") != 1) {
    5679                 I = groebner(I);
    5680         };
    5681 
    5682         ideal rad;                                // The radical
    5683         int allIndep = 1;                // All max independent sets are used
    5684 
    5685         list result = newRadicalReduction(I, ser, allIndep, #);
    5686         int done = result[3];
    5687         rad = result[1];
    5688         if (done == 0) {
    5689                 rad = intersect(rad, newRadicalKL(result[2], ideal(1), #));
    5690         };
    5691         return(rad);
    5692 };
    5693 
     5818static 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}
    56945842
    56955843///////////////////////////////////////////////////////////////////////////////
    … …  
    57065854// obtained in intermediate steps.
    57075855{
    5708         ideal rad = 1;
    5709         ideal equiRad = 1;
    5710         list primes;
    5711         int k;                        // Counter
    5712         int il;                 // If il = 1, only the equiradical is required.
    5713         int iDim;                // The dimension of I
    5714         int stop = 0;   // Checks if the radical has been obtained
    5715 
    5716         if (attrib(I, "isSB") != 1) {
    5717                 I = groebner(I);
    5718         };
    5719         iDim = dim(I);
    5720 
    5721     // Checks if only equiradical is required
    5722     if (size(#) > 0) {
    5723        il = #[1];
    5724     };
    5725 
    5726         while(stop == 0) {
    5727                 dbprint (printlevel-voice, "// We call radLoopR to find new prime ideals.");
    5728                 primes = newRadicalSLIteration(I, rad);                         // A list of primes or intersections of primes, not included in P
    5729                 dbprint (printlevel - voice, "// Output of Iteration Step:");
    5730                 dbprint (printlevel - voice, primes);
    5731                 if (size(primes) > 0) {
    5732                         dbprint (printlevel - voice, "// We intersect P with the ideal just obtained.");
    5733                         for(k = 1; k <= size(primes); k++) {
    5734                                 rad = intersect(rad, primes[k]);
    5735                                 if (il == 1){
    5736                                         if (attrib(primes[k], "isSB") != 1) {
    5737                                                 primes[k] = groebner(primes[k]);
    5738                                         };
    5739                                         if (iDim == dim(primes[k])) {
    5740                                                 equiRad = intersect(equiRad, primes[k]);
    5741                                         };
    5742                                 };
    5743                         };
    5744                 } else {
    5745                         stop = 1
    5746                 };
    5747         };
    5748         if (il == 0) {
    5749                 return(rad);
    5750         } else {
    5751                 return(equiRad);
    5752         };
    5753 };
     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}
    57545915
    57555916//////////////////////////////////////////////////////////////////////////
    … …  
    57665927//       This is not used in the new algorithm. It is part of KL algorithm
    57675928
    5768 static proc newRadicalReduction(ideal I, ideal ser, int allIndep, list #) {
     5929static proc newRadicalReduction(ideal I, ideal ser, int allIndep, list #)
     5930{
    57695931// allMaximal                1 -> Indicates that the reduction to the zerodim case
    57705932//                       must be done for all indep set of the leading terms ideal
    … …  
    57745936//                  It is used to set the value of done.)
    57755937
    5776    attrib(I, "isSB", 1);   // I needs to be a reduced standard basis
    5777    list indep, fett;
    5778    intvec @w, @hilb, op;
    5779    int @wr, @n, @m, lauf, di;
    5780    ideal fac, @h, collectrad, lsau;
    5781    poly @q;
    5782    string @va, quotring;
    5783 
    5784    def @P = basering;
    5785    int jdim = dim(I);               // Computes the dimension of I
    5786    int  homo = homog(I);            // Finds out if I is homogeneous
    5787    ideal rad = ideal(1);            // The unit ideal
    5788    ideal te = ser;
    5789    if(size(#) > 0)
    5790    {
    5791       @wr = #[1];
    5792    }
    5793    if(homo == 1)
    5794    {
    5795       for(@n = 1; @n <= nvars(basering); @n++)
    5796       {
    5797          @w[@n] = ord(var(@n));
    5798       }
    5799       @hilb = hilb(I, 1, @w);
    5800    }
    5801 
    5802    // SL 2006.04.11 1 Debug messages
    5803    dbprint(printlevel-voice, "//Computes the radical of the ideal:", I);
    5804    // SL 2006.04.11 2 Debug messages
     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
    58055967
    58065968  //---------------------------------------------------------------------------
    … …  
    58085970  //---------------------------------------------------------------------------
    58095971
    5810    if (jdim==-1)
    5811    {
    5812       return(ideal(1), ideal(1), 1);
    5813    }
     5972  if (jdim==-1)
     5973  {
     5974    return(ideal(1), ideal(1), 1);
     5975  }
    58145976
    58155977  //---------------------------------------------------------------------------
    … …  
    58175979  //---------------------------------------------------------------------------
    58185980
    5819    if (jdim==0)
    5820    {
    5821       return(zeroRad(I), ideal(1), 1);
    5822    }
    5823 
    5824    //-------------------------------------------------------------------------
    5825    //search for a maximal independent set indep,i.e.
    5826    //look for subring such that the intersection with the ideal is zero
    5827    //j intersected with K[var(indep[3]+1),...,var(nvar)] is zero,
    5828    //indep[1] is the new varstring, indep[2] the string for the block-ordering
    5829    //-------------------------------------------------------------------------
    5830 
    5831    // SL 2006-04-24 1         If allIndep = 0, then it only computes one maximal
    5832    //                                        independent set.
    5833    //                                        This looks better for the new algorithm but not for KL
    5834    //                                        algorithm
    5835    list parameters = allIndep;
    5836    indep = newMaxIndependSet(I, parameters);
    5837    // SL 2006-04-24 2
    5838 
    5839    for(@m = 1; @m <= size(indep); @m++)
    5840    {
    5841       if((indep[@m][1] == varstr(basering)) && (@m == 1))
    5842       //this is the good case, nothing to do, just to have the same notations
    5843       //change the ring
    5844       {
    5845          execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
     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)+"),("
    58466008                              +ordstr(basering)+");");
    5847          ideal @j = fetch(@P, I);
    5848          attrib(@j, "isSB", 1);
     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] + "),("
     6017                              + indep[@m][2] + ");");
     6018      execute("map phi = @P," + @va + ";");
     6019      if(homo == 1)
     6020      {
     6021        ideal @j = std(phi(I), @hilb, @w);
    58496022      }
    58506023      else
    58516024      {
    5852          @va = string(maxideal(1));
    5853 
    5854          execute("ring gnir1 = (" + charstr(basering) + "), (" + indep[@m][1] + "),("
    5855                               + indep[@m][2] + ");");
    5856          execute("map phi = @P," + @va + ";");
    5857          if(homo == 1)
    5858          {
    5859             ideal @j = std(phi(I), @hilb, @w);
    5860          }
    5861          else
    5862          {
    5863             ideal @j = groebner(phi(I));
    5864          }
    5865       }
    5866       if((deg(@j[1]) == 0) || (dim(@j) < jdim))
    5867       {
    5868          setring @P;
    5869          break;
    5870       }
    5871       for (lauf = 1; lauf <= size(@j); lauf++)
    5872       {
    5873          fett[lauf] = size(@j[lauf]);
    5874       }
    5875       //------------------------------------------------------------------------
    5876       // We have now the following situation:
    5877       // j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
    5878       // to this quotientring, j is there still a standardbasis, the
    5879       // leading coefficients of the polynomials there (polynomials in
    5880       // K[var(nnp+1),..,var(nva)]) are collected in the list h,
    5881       // we need their LCM, gh, because of the following:
    5882       // let (j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..rest..]
    5883       // intersected with K[var(1),...,var(nva)] is (j:gh^n)
    5884       // on the other hand j = ((j, gh^n) intersected with (j : gh^n))
    5885 
    5886       //------------------------------------------------------------------------
    5887       // The arrangement for the quotientring K(var(nnp+1),..,var(nva))[..rest..]
    5888       // and the map phi:K[var(1),...,var(nva)] ----->
    5889       // K(var(nnpr+1),..,var(nva))[..the rest..]
    5890       //------------------------------------------------------------------------
    5891       quotring = newPrepareQuotientRing(nvars(basering) - indep[@m][3]);
    5892       //------------------------------------------------------------------------
    5893       // We pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
    5894       //------------------------------------------------------------------------
    5895 
    5896       execute(quotring);
    5897 
    5898       // @j considered in the quotientring
    5899       ideal @j = imap(gnir1, @j);
    5900 
    5901       kill gnir1;
    5902 
    5903       // j is a standardbasis in the quotientring but usually not minimal
    5904       // here it becomes minimal
    5905 
    5906       @j = clearSB(@j, fett);
    5907 
    5908       // We need later LCM(h[1],...) = gh for saturation
    5909       ideal @h;
    5910       if(deg(@j[1]) > 0)
    5911       {
    5912          for(@n = 1; @n <= size(@j); @n++)
    5913          {
    5914             @h[@n] = leadcoef(@j[@n]);
    5915          }
    5916          op = option(get);
    5917          option(redSB);
    5918          @j = interred(@j);  //to obtain a reduced standardbasis
    5919          attrib(@j, "isSB", 1);
    5920          option(set, op);
    5921 
    5922          // SL 1 Debug messages
    5923          dbprint(printlevel - voice, "zero_rad", basering, @j, dim(groebner(@j)));
    5924          ideal zero_rad = zeroRad(@j);
    5925          dbprint(printlevel - voice, "zero_rad passed");
    5926          // SL 2
    5927       }
    5928       else
    5929       {
    5930          ideal zero_rad = ideal(1);
    5931       }
    5932 
    5933       // We need the intersection of the ideals in the list quprimary with the
    5934       // polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
    5935       // but fi polynomials, then the intersection of q with the polynomialring
    5936       // is the saturation of the ideal generated by f1,...,fr with respect to
    5937       // h which is the lcm of the leading coefficients of the fi considered in
     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
    59386098      // the quotientring: this is coded in saturn
    59396099
    … …  
    60266186   } else {
    60276187      done = 0;
    6028    };
     6188   }
    60296189   // SL 2006.04.21 2
    60306190
    … …  
    60336193   return(result);
    60346194   // SL 2006.04.21 2
    6035 };
     6195}
    60366196
    60376197
    … …  
    60556215                good = 1 - rad_con(P[k], I);
    60566216                k++;
    6057         };
     6217        }
    60586218        k--;
    60596219        if (good == 0) {
    … …  
    60616221                list emptyList = list();
    60626222                return (emptyList);
    6063         };
     6223        }
    60646224        dbprint(printlevel - voice, "// That one was good!");
    60656225        dbprint(printlevel - voice, "// We saturate I with respect to this element.");
    … …  
    60696229                dbprint(printlevel - voice, "// The polynomial is 1, the saturation in not actually computed.");
    60706230                ideal J = I;
    6071         };
     6231        }
    60726232
    60736233        // We now call proc radicalNew;
    … …  
    60826242
    60836243        return(result[1]);
    6084 };
     6244}
    60856245
    60866246
    … …  
    61186278   } else {
    61196279           allMaximal = 1;
    6120    };
     6280   }
    61216281
    61226282   int nMax;
    … …  
    61256285   } else {
    61266286          nMax = 1;
    6127    };
     6287   }
    61286288
    61296289   for(n = 1; n <= nMax; n++)
    … …  
    64556615   pr;
    64566616}
     6617///////////////////////////////////////////////////////////////////////////////
     6618static proc newDecompStep(ideal i, list #)
     6619"USAGE:  newDecompStep(i); i ideal  (for primary decomposition)
     6620         newDecompStep(i,1);        (for the associated primes of dimension of i)
     6621         newDecompStep(i,2);        (for the minimal associated primes)
     6622         newDecompStep(i,3);        (for the absolute primary decomposition)
     6623         "oneIndep";        (for using only one max indep set)
     6624         "intersect";        (returns alse the intersection of the components founded)
     6625
     6626RETURN:  list = list of primary ideals and their associated primes
     6627         (at even positions in the list)
     6628         (resp. a list of the minimal associated primes)
     6629NOTE:    Algorithm of Gianni/Trager/Zacharias
     6630EXAMPLE: example newDecompStep; shows an example
     6631"
     6632{
     6633  intvec op,@vv;
     6634  def  @P = basering;
     6635  list primary,indep,ltras;
     6636  intvec @vh,isat,@w;
     6637  int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi,abspri,ab,nn;
     6638  ideal peek=i;
     6639  ideal ser,tras;
     6640  list data;
     6641  list result;
     6642  intvec @hilb;
     6643  int isS=(attrib(i,"isSB")==1);
     6644
     6645  // Debug
     6646  dbprint(printlevel - voice, "newDecompStep, v2.0");
     6647
     6648  string indepOption = "allIndep";
     6649  string intersectOption = "noIntersect";
     6650
     6651  if(size(#)>0)
     6652  {
     6653   int count = 1;
     6654   if(typeof(#[count]) == "string") {
     6655     if ((#[count] == "oneIndep") or (#[count] == "allIndep")){
     6656       indepOption = #[count];
     6657        count++;
     6658     }
     6659   }
     6660   if(typeof(#[count]) == "string") {
     6661     if ((#[count] == "intersect") or (#[count] == "noIntersect")){
     6662       intersectOption = #[count];
     6663        count++;
     6664     }
     6665   }
     6666     if((typeof(#[count]) == "int") or (typeof(#[count]) == "number"))
     6667     {
     6668       if ((#[count]==1)||(#[count]==2)||(#[count]==3))
     6669       {
     6670          @wr=#[count];
     6671          if(@wr==3){abspri = 1; @wr = 0;}
     6672          count++;
     6673       }
     6674     }
     6675     if(size(#)>count)
     6676     {
     6677          seri=1;
     6678          peek=#[count + 1];
     6679          ser=#[count + 2];
     6680     }
     6681  }
     6682  if(abspri)
     6683  {
     6684     list absprimary,abskeep,absprimarytmp,abskeeptmp;
     6685  }
     6686  homo=homog(i);
     6687  if(homo==1)
     6688  {
     6689    if(attrib(i,"isSB")!=1)
     6690    {
     6691      //ltras=mstd(i);
     6692      tras=groebner(i);
     6693      ltras=tras,tras;
     6694      attrib(ltras[1],"isSB",1);
     6695    }
     6696    else
     6697    {
     6698      ltras=i,i;
     6699      attrib(ltras[1],"isSB",1);
     6700    }
     6701    tras = ltras[1];
     6702    attrib(tras,"isSB",1);
     6703    if(dim(tras)==0)
     6704    {
     6705        primary[1]=ltras[2];
     6706        primary[2]=maxideal(1);
     6707        if(@wr>0)
     6708        {
     6709           list l;
     6710           l[1]=maxideal(1);
     6711           l[2]=maxideal(1);
     6712       if (intersectOption == "intersect") {
     6713         return(list(l, maxideal(1)));
     6714       } else {
     6715         return(l);
     6716       }
     6717        }
     6718        if (intersectOption == "intersect") {
     6719         return(list(primary, primary[1]));
     6720       } else {
     6721         return(primary);
     6722       }
     6723
     6724     }
     6725     for(@n=1;@n<=nvars(basering);@n++)
     6726     {
     6727        @w[@n]=ord(var(@n));
     6728     }
     6729     @hilb=hilb(tras,1,@w);
     6730     intvec keephilb=@hilb;
     6731  }
     6732
     6733  //----------------------------------------------------------------
     6734  //i is the zero-ideal
     6735  //----------------------------------------------------------------
     6736
     6737  if(size(i)==0)
     6738  {
     6739     primary=i,i;
     6740    if (intersectOption == "intersect") {
     6741       return(list(primary, i));
     6742    } else {
     6743       return(primary);
     6744   }
     6745  }
     6746
     6747  //----------------------------------------------------------------
     6748  //pass to the lexicographical ordering and compute a standardbasis
     6749  //----------------------------------------------------------------
     6750
     6751  int lp=islp();
     6752
     6753  execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);");
     6754  op=option(get);
     6755  option(redSB);
     6756
     6757  ideal ser=fetch(@P,ser);
     6758  if(homo==1)
     6759  {
     6760     if(!lp)
     6761     {
     6762       ideal @j=std(fetch(@P,i),@hilb,@w);
     6763     }
     6764     else
     6765     {
     6766        ideal @j=fetch(@P,tras);
     6767        attrib(@j,"isSB",1);
     6768     }
     6769  }
     6770  else
     6771  {
     6772      if(lp&&isS)
     6773      {
     6774        ideal @j=fetch(@P,i);
     6775        attrib(@j,"isSB",1);
     6776      }
     6777      else
     6778      {
     6779        ideal @j=groebner(fetch(@P,i));
     6780      }
     6781  }
     6782  option(set,op);
     6783  if(seri==1)
     6784  {
     6785    ideal peek=fetch(@P,peek);
     6786    attrib(peek,"isSB",1);
     6787  }
     6788  else
     6789  {
     6790    ideal peek=@j;
     6791  }
     6792  if((size(ser)==0)&&(!abspri))
     6793  {
     6794    ideal fried;
     6795    @n=size(@j);
     6796    for(@k=1;@k<=@n;@k++)
     6797    {
     6798      if(deg(lead(@j[@k]))==1)
     6799      {
     6800        fried[size(fried)+1]=@j[@k];
     6801        @j[@k]=0;
     6802      }
     6803    }
     6804    if(size(fried)==nvars(basering))
     6805    {
     6806       setring @P;
     6807       primary[1]=i;
     6808       primary[2]=i;
     6809     if (intersectOption == "intersect") {
     6810       return(list(primary, i));
     6811     } else {
     6812       return(primary);
     6813     }
     6814    }
     6815    if(size(fried)>0)
     6816    {
     6817       string newva;
     6818       string newma;
     6819       for(@k=1;@k<=nvars(basering);@k++)
     6820       {
     6821          @n1=0;
     6822          for(@n=1;@n<=size(fried);@n++)
     6823          {
     6824             if(leadmonom(fried[@n])==var(@k))
     6825             {
     6826                @n1=1;
     6827                break;
     6828             }
     6829          }
     6830          if(@n1==0)
     6831          {
     6832            newva=newva+string(var(@k))+",";
     6833            newma=newma+string(var(@k))+",";
     6834          }
     6835          else
     6836          {
     6837            newma=newma+string(0)+",";
     6838          }
     6839       }
     6840       newva[size(newva)]=")";
     6841       newma[size(newma)]=";";
     6842       execute("ring @deirf=("+charstr(gnir)+"),("+newva+",lp;");
     6843       execute("map @kappa=gnir,"+newma);
     6844       ideal @j= @kappa(@j);
     6845       @j=simplify(@j, 2);
     6846       attrib(@j,"isSB",1);
     6847       result = newDecompStep(@j, indepOption, intersectOption, @wr);
     6848       if (intersectOption == "intersect") {
     6849         list pr = result[1];
     6850         ideal intersection = result[2];
     6851       } else {
     6852         list pr = result;
     6853       }
     6854
     6855       setring gnir;
     6856       list pr=imap(@deirf,pr);
     6857       for(@k=1;@k<=size(pr);@k++)
     6858       {
     6859         @j=pr[@k]+fried;
     6860         pr[@k]=@j;
     6861       }
     6862       if (intersectOption == "intersect") {
     6863      ideal intersection = imap(@deirf, intersection);
     6864      @j = intersection + fried;
     6865      intersection = @j;
     6866     }
     6867     setring @P;
     6868     if (intersectOption == "intersect") {
     6869       return(list(imap(gnir,pr), imap(gnir,intersection)));
     6870     } else {
     6871       return(imap(gnir,pr));
     6872     }
     6873    }
     6874  }
     6875  //----------------------------------------------------------------
     6876  //j is the ring
     6877  //----------------------------------------------------------------
     6878
     6879  if (dim(@j)==-1)
     6880  {
     6881    setring @P;
     6882    primary=ideal(1),ideal(1);
     6883   if (intersectOption == "intersect") {
     6884     return(list(primary, ideal(1)));
     6885   } else {
     6886     return(primary);
     6887   }
     6888  }
     6889
     6890  //----------------------------------------------------------------
     6891  //  the case of one variable
     6892  //----------------------------------------------------------------
     6893
     6894  if(nvars(basering)==1)
     6895  {
     6896
     6897     list fac=factor(@j[1]);
     6898     list gprimary;
     6899     poly generator;
     6900     ideal gIntersection;
     6901     for(@k=1;@k<=size(fac[1]);@k++)
     6902     {
     6903        if(@wr==0)
     6904        {
     6905           gprimary[2*@k-1]=ideal(fac[1][@k]^fac[2][@k]);
     6906           gprimary[2*@k]=ideal(fac[1][@k]);
     6907        }
     6908        else
     6909        {
     6910           gprimary[2*@k-1]=ideal(fac[1][@k]);
     6911           gprimary[2*@k]=ideal(fac[1][@k]);
     6912        }
     6913        if (intersectOption == "intersect") {
     6914          generator = generator * fac[1][@k];
     6915        }
     6916     }
     6917   if (intersectOption == "intersect") {
     6918       gIntersection = generator;
     6919   }
     6920
     6921     setring @P;
     6922     primary=fetch(gnir,gprimary);
     6923   if (intersectOption == "intersect") {
     6924       ideal intersection = fetch(gnir,gIntersection);
     6925     }
     6926
     6927//HIER
     6928    if(abspri)
     6929    {
     6930       list resu,tempo;
     6931       string absotto;
     6932       for(ab=1;ab<=size(primary)/2;ab++)
     6933       {
     6934          absotto= absFactorize(primary[2*ab][1],77);
     6935          tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
     6936          resu[ab]=tempo;
     6937       }
     6938       primary=resu;
     6939       intersection = 1;
     6940       for(ab=1;ab<=size(primary);ab++) {
     6941          intersection = intersect(intersection, primary[ab][2]);
     6942        }
     6943     }
     6944     if (intersectOption == "intersect") {
     6945     return(list(primary, intersection));
     6946     } else {
     6947     return(primary);
     6948     }
     6949  }
     6950
     6951 //------------------------------------------------------------------
     6952 //the zero-dimensional case
     6953 //------------------------------------------------------------------
     6954  if (dim(@j)==0)
     6955  {
     6956    op=option(get);
     6957    option(redSB);
     6958    list gprimary= newZero_decomp(@j,ser,@wr);
     6959
     6960    setring @P;
     6961    primary=fetch(gnir,gprimary);
     6962
     6963    if(size(ser)>0)
     6964    {
     6965      primary=cleanPrimary(primary);
     6966    }
     6967//HIER
     6968    if(abspri)
     6969    {
     6970       list resu,tempo;
     6971       string absotto;
     6972       for(ab=1;ab<=size(primary)/2;ab++)
     6973       {
     6974          absotto= absFactorize(primary[2*ab][1],77);
     6975          tempo=primary[2*ab-1],primary[2*ab],absotto,string(var(nvars(basering)));
     6976          resu[ab]=tempo;
     6977       }
     6978       primary=resu;
     6979    }
     6980   if (intersectOption == "intersect") {
     6981     return(list(primary, fetch(gnir,@j)));
     6982   } else {
     6983     return(primary);
     6984   }
     6985
     6986  }
     6987
     6988  poly @gs,@gh,@p;
     6989  string @va,quotring;
     6990  list quprimary,htprimary,collectprimary,lsau,lnew,allindep,restindep;
     6991  ideal @h;
     6992  int jdim=dim(@j);
     6993  list fett;
     6994  int lauf,di,newtest;
     6995  //------------------------------------------------------------------
     6996  //search for a maximal independent set indep,i.e.
     6997  //look for subring such that the intersection with the ideal is zero
     6998  //j intersected with K[var(indep[3]+1),...,var(nvar] is zero,
     6999  //indep[1] is the new varstring and indep[2] the string for block-ordering
     7000  //------------------------------------------------------------------
     7001  if(@wr!=1)
     7002  {
     7003   allindep = newMaxIndependSetLp(@j, indepOption);
     7004     for(@m=1;@m<=size(allindep);@m++)
     7005     {
     7006        if(allindep[@m][3]==jdim)
     7007        {
     7008           di++;
     7009           indep[di]=allindep[@m];
     7010        }
     7011        else
     7012        {
     7013           lauf++;
     7014           restindep[lauf]=allindep[@m];
     7015        }
     7016     }
     7017   }
     7018   else
     7019   {
     7020      indep = newMaxIndependSetLp(@j, indepOption);
     7021   }
     7022
     7023  ideal jkeep=@j;
     7024  if(ordstr(@P)[1]=="w")
     7025  {
     7026     execute("ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),("+ordstr(@P)+");");
     7027  }
     7028  else
     7029  {
     7030     execute( "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);");
     7031  }
     7032
     7033  if(homo==1)
     7034  {
     7035    if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w")
     7036       ||(ordstr(@P)[3]=="w"))
     7037    {
     7038      ideal jwork=imap(@P,tras);
     7039      attrib(jwork,"isSB",1);
     7040    }
     7041    else
     7042    {
     7043      ideal jwork=std(imap(gnir,@j),@hilb,@w);
     7044    }
     7045
     7046  }
     7047  else
     7048  {
     7049    ideal jwork=groebner(imap(gnir,@j));
     7050  }
     7051  list hquprimary;
     7052  poly @p,@q;
     7053  ideal @h,fac,ser;
     7054//Aenderung================
     7055  ideal @Ptest=1;
     7056//=========================
     7057  di=dim(jwork);
     7058  keepdi=di;
     7059
     7060  ser = 1;
     7061
     7062  setring gnir;
     7063  for(@m=1; @m<=size(indep); @m++)
     7064  {
     7065    data[1] = indep[@m];
     7066    result = newReduction(@j, ser, @hilb, @w, jdim, abspri, @wr, data);
     7067    quprimary = quprimary + result[1];
     7068    if(abspri) {
     7069      absprimary = absprimary + result[2];
     7070      abskeep = abskeep + result[3];
     7071      }
     7072    @h = result[5];
     7073    ser = result[4];
     7074     if(size(@h)>0)
     7075     {
     7076           //---------------------------------------------------------------
     7077           //we change to @Phelp to have the ordering dp for saturation
     7078           //---------------------------------------------------------------
     7079
     7080        setring @Phelp;
     7081        @h=imap(gnir,@h);
     7082//Aenderung==================================
     7083        if(defined(@LL)){kill @LL;}
     7084        list @LL=minSat(jwork,@h);
     7085        @Ptest=intersect(@Ptest,@LL[1]);
     7086        ser = intersect(ser, @LL[1]);
     7087//===========================================
     7088
     7089        if(@wr!=1)
     7090        {
     7091//Aenderung==================================
     7092          @q=@LL[2];
     7093//===========================================
     7094         //@q=minSat(jwork,@h)[2];
     7095        }
     7096        else
     7097        {
     7098            fac=ideal(0);
     7099            for(lauf=1;lauf<=ncols(@h);lauf++)
     7100           {
     7101              if(deg(@h[lauf])>0)
     7102              {
     7103                 fac=fac+factorize(@h[lauf],1);
     7104              }
     7105           }
     7106           fac=simplify(fac,4);
     7107           @q=1;
     7108           for(lauf=1;lauf<=size(fac);lauf++)
     7109           {
     7110             @q=@q*fac[lauf];
     7111           }
     7112        }
     7113        jwork = std(jwork,@q);
     7114        keepdi = dim(jwork);
     7115        if(keepdi < di)
     7116        {
     7117           setring gnir;
     7118           @j = imap(@Phelp, jwork);
     7119           ser = imap(@Phelp, ser);
     7120           break;
     7121        }
     7122        if(homo == 1)
     7123        {
     7124           @hilb = hilb(jwork, 1, @w);
     7125        }
     7126
     7127        setring gnir;
     7128        ser = imap(@Phelp, ser);
     7129        @j = imap(@Phelp, jwork);
     7130     }
     7131  }
     7132
     7133  if((size(quprimary)==0)&&(@wr==1))
     7134  {
     7135     @j=ideal(1);
     7136     quprimary[1]=ideal(1);
     7137     quprimary[2]=ideal(1);
     7138  }
     7139  if((size(quprimary)==0))
     7140  {
     7141    keepdi = di - 1;
     7142    quprimary[1]=ideal(1);
     7143    quprimary[2]=ideal(1);
     7144  }
     7145  //---------------------------------------------------------------
     7146  //notice that j=sat(j,gh) intersected with (j,gh^n)
     7147  //we finished with sat(j,gh) and have to start with (j,gh^n)
     7148  //---------------------------------------------------------------
     7149  if((deg(@j[1])!=0)&&(@wr!=1))
     7150  {
     7151     if(size(quprimary)>0)
     7152     {
     7153        setring @Phelp;
     7154        ser=imap(gnir,ser);
     7155
     7156        hquprimary=imap(gnir,quprimary);
     7157        if(@wr==0)
     7158        {
     7159//Aenderung====================================================
     7160//HIER STATT DURCHSCHNITT SATURIEREN!
     7161           ideal htest=@Ptest;
     7162/*
     7163           ideal htest=hquprimary[1];
     7164           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
     7165           {
     7166              htest=intersect(htest,hquprimary[2*@n1-1]);
     7167           }
     7168*/
     7169//=============================================================
     7170        }
     7171        else
     7172        {
     7173           ideal htest=hquprimary[2];
     7174
     7175           for (@n1=2;@n1<=size(hquprimary)/2;@n1++)
     7176           {
     7177              htest=intersect(htest,hquprimary[2*@n1]);
     7178           }
     7179        }
     7180
     7181        if(size(ser)>0)
     7182        {
     7183           ser=intersect(htest,ser);
     7184        }
     7185        else
     7186        {
     7187          ser=htest;
     7188        }
     7189        setring gnir;
     7190        ser=imap(@Phelp,ser);
     7191     }
     7192     if(size(reduce(ser,peek,1))!=0)
     7193     {
     7194    for(@m=1;@m<=size(restindep);@m++)
     7195        {
     7196         // if(restindep[@m][3]>=keepdi)
     7197         // {
     7198           isat=0;
     7199           @n2=0;
     7200
     7201           if(restindep[@m][1]==varstr(basering))
     7202           //the good case, nothing to do, just to have the same notations
     7203           //change the ring
     7204           {
     7205              execute("ring gnir1 = ("+charstr(basering)+"),("+
     7206                varstr(basering)+"),("+ordstr(basering)+");");
     7207              ideal @j=fetch(gnir,jkeep);
     7208              attrib(@j,"isSB",1);
     7209           }
     7210           else
     7211           {
     7212              @va=string(maxideal(1));
     7213              execute("ring gnir1 = ("+charstr(basering)+"),("+
     7214                      restindep[@m][1]+"),(" +restindep[@m][2]+");");
     7215              execute("map phi=gnir,"+@va+";");
     7216              op=option(get);
     7217              option(redSB);
     7218              if(homo==1)
     7219              {
     7220                 ideal @j=std(phi(jkeep),keephilb,@w);
     7221              }
     7222              else
     7223              {
     7224                ideal @j=groebner(phi(jkeep));
     7225              }
     7226              ideal ser=phi(ser);
     7227              option(set,op);
     7228           }
     7229
     7230           for (lauf=1;lauf<=size(@j);lauf++)
     7231           {
     7232              fett[lauf]=size(@j[lauf]);
     7233           }
     7234           //------------------------------------------------------------------
     7235           //we have now the following situation:
     7236           //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may
     7237           //pass to this quotientring, j is their still a standardbasis, the
     7238           //leading coefficients of the polynomials  there (polynomials in
     7239           //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     7240           //we need their ggt, gh, because of the following:
     7241           //let (j:gh^n)=(j:gh^infinity) then
     7242           //j*K(var(nnp+1),..,var(nva))[..the rest..]
     7243           //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     7244           //on the other hand j=(j,gh^n) intersected with (j:gh^n)
     7245
     7246           //------------------------------------------------------------------
     7247
     7248           //the arrangement for the quotientring
     7249           // K(var(nnp+1),..,var(nva))[..the rest..]
     7250           //and the map phi:K[var(1),...,var(nva)] ---->
     7251           //--->K(var(nnpr+1),..,var(nva))[..the rest..]
     7252           //------------------------------------------------------------------
     7253
     7254           quotring=prepareQuotientring(nvars(basering)-restindep[@m][3]);
     7255
     7256           //------------------------------------------------------------------
     7257           //we pass to the quotientring  K(var(nnp+1),..,var(nva))[..rest..]
     7258           //------------------------------------------------------------------
     7259
     7260           execute(quotring);
     7261
     7262           // @j considered in the quotientring
     7263           ideal @j=imap(gnir1,@j);
     7264           ideal ser=imap(gnir1,ser);
     7265
     7266           kill gnir1;
     7267
     7268           //j is a standardbasis in the quotientring but usually not minimal
     7269           //here it becomes minimal
     7270           @j=clearSB(@j,fett);
     7271           attrib(@j,"isSB",1);
     7272
     7273           //we need later ggt(h[1],...)=gh for saturation
     7274           ideal @h;
     7275
     7276           for(@n=1;@n<=size(@j);@n++)
     7277           {
     7278              @h[@n]=leadcoef(@j[@n]);
     7279           }
     7280           //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..rest..]
     7281
     7282           op=option(get);
     7283           option(redSB);
     7284           list uprimary= newZero_decomp(@j,ser,@wr);
     7285//HIER
     7286           if(abspri)
     7287           {
     7288              ideal II;
     7289              ideal jmap;
     7290              map sigma;
     7291              nn=nvars(basering);
     7292              map invsigma=basering,maxideal(1);
     7293              for(ab=1;ab<=size(uprimary)/2;ab++)
     7294              {
     7295                 II=uprimary[2*ab];
     7296                 attrib(II,"isSB",1);
     7297                 if(deg(II[1])!=vdim(II))
     7298                 {
     7299                    jmap=randomLast(50);
     7300                    sigma=basering,jmap;
     7301                    jmap[nn]=2*var(nn)-jmap[nn];
     7302                    invsigma=basering,jmap;
     7303                    II=groebner(sigma(II));
     7304                  }
     7305                  absprimarytmp[ab]= absFactorize(II[1],77);
     7306                  II=var(nn);
     7307                  abskeeptmp[ab]=string(invsigma(II));
     7308                  invsigma=basering,maxideal(1);
     7309              }
     7310           }
     7311           option(set,op);
     7312
     7313           //we need the intersection of the ideals in the list quprimary with
     7314           //the polynomialring, i.e. let q=(f1,...,fr) in the quotientring
     7315           //such an ideal but fi polynomials, then the intersection of q with
     7316           //the polynomialring is the saturation of the ideal generated by
     7317           //f1,...,fr with respect toh which is the lcm of the leading
     7318           //coefficients of the fi considered in the quotientring:
     7319           //this is coded in saturn
     7320
     7321           list saturn;
     7322           ideal hpl;
     7323
     7324           for(@n=1;@n<=size(uprimary);@n++)
     7325           {
     7326              hpl=0;
     7327              for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
     7328              {
     7329                 hpl=hpl,leadcoef(uprimary[@n][@n1]);
     7330              }
     7331               saturn[@n]=hpl;
     7332           }
     7333           //------------------------------------------------------------------
     7334           //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..rest..]
     7335           //back to the polynomialring
     7336           //------------------------------------------------------------------
     7337           setring gnir;
     7338           collectprimary=imap(quring,uprimary);
     7339           lsau=imap(quring,saturn);
     7340           @h=imap(quring,@h);
     7341
     7342           kill quring;
     7343
     7344
     7345           @n2=size(quprimary);
     7346//================NEU=========================================
     7347           if(deg(quprimary[1][1])<=0){ @n2=0; }
     7348//============================================================
     7349
     7350           @n3=@n2;
     7351
     7352           for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
     7353           {
     7354              if(deg(collectprimary[2*@n1][1])>0)
     7355              {
     7356                 @n2++;
     7357                 quprimary[@n2]=collectprimary[2*@n1-1];
     7358                 lnew[@n2]=lsau[2*@n1-1];
     7359                 @n2++;
     7360                 lnew[@n2]=lsau[2*@n1];
     7361                 quprimary[@n2]=collectprimary[2*@n1];
     7362                 if(abspri)
     7363                 {
     7364                   absprimary[@n2/2]=absprimarytmp[@n1];
     7365                   abskeep[@n2/2]=abskeeptmp[@n1];
     7366                 }
     7367              }
     7368           }
     7369
     7370
     7371           //here the intersection with the polynomialring
     7372           //mentioned above is really computed
     7373
     7374           for(@n=@n3/2+1;@n<=@n2/2;@n++)
     7375           {
     7376              if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
     7377              {
     7378                 quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7379                 quprimary[2*@n]=quprimary[2*@n-1];
     7380              }
     7381              else
     7382              {
     7383                 if(@wr==0)
     7384                 {
     7385                    quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7386                 }
     7387                 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
     7388              }
     7389           }
     7390           if(@n2>=@n3+2)
     7391           {
     7392              setring @Phelp;
     7393              ser=imap(gnir,ser);
     7394              hquprimary=imap(gnir,quprimary);
     7395              for(@n=@n3/2+1;@n<=@n2/2;@n++)
     7396              {
     7397                if(@wr==0)
     7398                {
     7399                   ser=intersect(ser,hquprimary[2*@n-1]);
     7400                }
     7401                else
     7402                {
     7403                   ser=intersect(ser,hquprimary[2*@n]);
     7404                }
     7405              }
     7406              setring gnir;
     7407              ser=imap(@Phelp,ser);
     7408           }
     7409
     7410         // }
     7411        }
     7412//HIER
     7413        if(abspri)
     7414        {
     7415          list resu,tempo;
     7416          for(ab=1;ab<=size(quprimary)/2;ab++)
     7417          {
     7418             if (deg(quprimary[2*ab][1])!=0)
     7419             {
     7420               tempo=quprimary[2*ab-1],quprimary[2*ab],
     7421                         absprimary[ab],abskeep[ab];
     7422               resu[ab]=tempo;
     7423             }
     7424          }
     7425          quprimary=resu;
     7426          @wr=3;
     7427        }
     7428        if(size(reduce(ser,peek,1))!=0)
     7429        {
     7430           if(@wr>0)
     7431           {
     7432              // The following line was dropped to avoid the recursion step:
     7433              //htprimary=newDecompStep(@j,@wr,peek,ser);
     7434              htprimary = list();
     7435           }
     7436           else
     7437           {
     7438              // The following line was dropped to avoid the recursion step:
     7439              //htprimary=newDecompStep(@j,peek,ser);
     7440              htprimary = list();
     7441           }
     7442           // here we collect now both results primary(sat(j,gh))
     7443           // and primary(j,gh^n)
     7444           @n=size(quprimary);
     7445           if (deg(quprimary[1][1])<=0) { @n=0; }
     7446           for (@k=1;@k<=size(htprimary);@k++)
     7447           {
     7448              quprimary[@n+@k]=htprimary[@k];
     7449           }
     7450        }
     7451     }
     7452   }
     7453   else
     7454   {
     7455      if(abspri)
     7456      {
     7457        list resu,tempo;
     7458        for(ab=1;ab<=size(quprimary)/2;ab++)
     7459        {
     7460           tempo=quprimary[2*ab-1],quprimary[2*ab],
     7461                   absprimary[ab],abskeep[ab];
     7462           resu[ab]=tempo;
     7463        }
     7464        quprimary=resu;
     7465      }
     7466   }
     7467  //---------------------------------------------------------------------------
     7468  //back to the ring we started with
     7469  //the final result: primary
     7470  //---------------------------------------------------------------------------
     7471
     7472  setring @P;
     7473  primary=imap(gnir,quprimary);
     7474
     7475   if (intersectOption == "intersect") {
     7476     return(list(primary, imap(gnir, ser)));
     7477   } else {
     7478     return(primary);
     7479   }
     7480}
     7481example
     7482{ "EXAMPLE:"; echo = 2;
     7483   ring  r = 32003,(x,y,z),lp;
     7484   poly  p = z2+1;
     7485   poly  q = z4+2;
     7486   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     7487   list pr= newDecompStep(i);
     7488   pr;
     7489   testPrimary( pr, i);
     7490}
     7491
     7492//
     7493proc newReduction(ideal @j, ideal ser, intvec @hilb, intvec @w, int jdim, int abspri, int @wr, list data)
     7494{
     7495   string @va;
     7496   string quotring;
     7497   intvec op;
     7498   intvec @vv;
     7499   def gnir = basering;
     7500     ideal isat=0;
     7501     int @n;
     7502     int @n1 = 0;
     7503     int @n2 = 0;
     7504     int @n3 = 0;
     7505     int homo = homog(@j);
     7506     int lauf;
     7507     int @k;
     7508     list fett;
     7509     int keepdi;
     7510     list collectprimary;
     7511     list lsau;
     7512     list lnew;
     7513     ideal @h;
     7514
     7515     list indepInfo = data[1];
     7516     list quprimary = list();
     7517
     7518
     7519  //if(abspri)
     7520  //{
     7521     int ab;
     7522     list absprimarytmp,abskeeptmp;
     7523     list absprimary, abskeep;
     7524  //}
     7525  // Debug
     7526  dbprint(printlevel - voice, "newReduction, v2.0");
     7527
     7528     if((indepInfo[1]==varstr(basering)))  // &&(@m==1)
     7529     //this is the good case, nothing to do, just to have the same notations
     7530     //change the ring
     7531     {
     7532        execute("ring gnir1 = ("+charstr(basering)+"),("+varstr(basering)+"),("
     7533                              +ordstr(basering)+");");
     7534        ideal @j = fetch(gnir, @j);
     7535        attrib(@j,"isSB",1);
     7536        ideal ser = fetch(gnir, ser);
     7537     }
     7538     else
     7539     {
     7540        @va=string(maxideal(1));
     7541//Aenderung==============
     7542        //if(@m==1)
     7543        //{
     7544           //@j=fetch(@P,i);
     7545        //}
     7546//=======================
     7547        execute("ring gnir1 = ("+charstr(basering)+"),("+indepInfo[1]+"),("
     7548                              +indepInfo[2]+");");
     7549        execute("map phi=gnir,"+@va+";");
     7550        op=option(get);
     7551        option(redSB);
     7552        if(homo==1)
     7553        {
     7554           ideal @j=std(phi(@j),@hilb,@w);
     7555        }
     7556        else
     7557        {
     7558          ideal @j=groebner(phi(@j));
     7559        }
     7560        ideal ser=phi(ser);
     7561
     7562        option(set,op);
     7563     }
     7564     if((deg(@j[1])==0)||(dim(@j)<jdim))
     7565     {
     7566       setring gnir;
     7567       break;
     7568     }
     7569     for (lauf=1;lauf<=size(@j);lauf++)
     7570     {
     7571         fett[lauf]=size(@j[lauf]);
     7572     }
     7573     //------------------------------------------------------------------------
     7574     //we have now the following situation:
     7575     //j intersected with K[var(nnp+1),..,var(nva)] is zero so we may pass
     7576     //to this quotientring, j is their still a standardbasis, the
     7577     //leading coefficients of the polynomials  there (polynomials in
     7578     //K[var(nnp+1),..,var(nva)]) are collected in the list h,
     7579     //we need their ggt, gh, because of the following: let
     7580     //(j:gh^n)=(j:gh^infinity) then j*K(var(nnp+1),..,var(nva))[..the rest..]
     7581     //intersected with K[var(1),...,var(nva)] is (j:gh^n)
     7582     //on the other hand j=(j,gh^n) intersected with (j:gh^n)
     7583
     7584     //------------------------------------------------------------------------
     7585
     7586     //arrangement for quotientring K(var(nnp+1),..,var(nva))[..the rest..] and
     7587     //map phi:K[var(1),...,var(nva)] --->K(var(nnpr+1),..,var(nva))[..rest..]
     7588     //------------------------------------------------------------------------
     7589
     7590     quotring=prepareQuotientring(nvars(basering)-indepInfo[3]);
     7591
     7592     //---------------------------------------------------------------------
     7593     //we pass to the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     7594     //---------------------------------------------------------------------
     7595
     7596     ideal @jj=lead(@j);               //!! vorn vereinbaren
     7597     execute(quotring);
     7598
     7599     ideal @jj=imap(gnir1,@jj);
     7600     @vv=clearSBNeu(@jj,fett);  //!! vorn vereinbaren
     7601     setring gnir1;
     7602     @k=size(@j);
     7603     for (lauf=1;lauf<=@k;lauf++)
     7604     {
     7605         if(@vv[lauf]==1)
     7606         {
     7607            @j[lauf]=0;
     7608         }
     7609     }
     7610     @j=simplify(@j,2);
     7611     setring quring;
     7612      // @j considered in the quotientring
     7613     ideal @j=imap(gnir1,@j);
     7614
     7615     ideal ser=imap(gnir1,ser);
     7616
     7617     kill gnir1;
     7618
     7619     //j is a standardbasis in the quotientring but usually not minimal
     7620     //here it becomes minimal
     7621
     7622     attrib(@j,"isSB",1);
     7623
     7624     //we need later ggt(h[1],...)=gh for saturation
     7625     ideal @h;
     7626     if(deg(@j[1])>0)
     7627     {
     7628        for(@n=1;@n<=size(@j);@n++)
     7629        {
     7630           @h[@n]=leadcoef(@j[@n]);
     7631        }
     7632        //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..]
     7633        op=option(get);
     7634        option(redSB);
     7635
     7636       int zeroMinAss = @wr;
     7637       if (@wr == 2) {zeroMinAss = 1;}
     7638        list uprimary= newZero_decomp(@j, ser, zeroMinAss);
     7639
     7640//HIER
     7641        if(abspri)
     7642        {
     7643           ideal II;
     7644           ideal jmap;
     7645           map sigma;
     7646           nn=nvars(basering);
     7647           map invsigma=basering,maxideal(1);
     7648           for(ab=1;ab<=size(uprimary)/2;ab++)
     7649           {
     7650              II=uprimary[2*ab];
     7651              attrib(II,"isSB",1);
     7652              if(deg(II[1])!=vdim(II))
     7653              {
     7654                 jmap=randomLast(50);
     7655                 sigma=basering,jmap;
     7656                 jmap[nn]=2*var(nn)-jmap[nn];
     7657                 invsigma=basering,jmap;
     7658                 II=groebner(sigma(II));
     7659               }
     7660               absprimarytmp[ab]= absFactorize(II[1],77);
     7661               II=var(nn);
     7662               abskeeptmp[ab]=string(invsigma(II));
     7663               invsigma=basering,maxideal(1);
     7664           }
     7665        }
     7666        option(set,op);
     7667     }
     7668     else
     7669     {
     7670       list uprimary;
     7671       uprimary[1]=ideal(1);
     7672       uprimary[2]=ideal(1);
     7673     }
     7674     //we need the intersection of the ideals in the list quprimary with the
     7675     //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal
     7676     //but fi polynomials, then the intersection of q with the polynomialring
     7677     //is the saturation of the ideal generated by f1,...,fr with respect to
     7678     //h which is the lcm of the leading coefficients of the fi considered in
     7679     //in the quotientring: this is coded in saturn
     7680
     7681     list saturn;
     7682     ideal hpl;
     7683
     7684     for(@n=1;@n<=size(uprimary);@n++)
     7685     {
     7686        uprimary[@n]=interred(uprimary[@n]); // temporary fix
     7687        hpl=0;
     7688        for(@n1=1;@n1<=size(uprimary[@n]);@n1++)
     7689        {
     7690           hpl=hpl,leadcoef(uprimary[@n][@n1]);
     7691        }
     7692        saturn[@n]=hpl;
     7693     }
     7694
     7695     //--------------------------------------------------------------------
     7696     //we leave  the quotientring   K(var(nnp+1),..,var(nva))[..the rest..]
     7697     //back to the polynomialring
     7698     //---------------------------------------------------------------------
     7699     setring gnir;
     7700
     7701     collectprimary=imap(quring,uprimary);
     7702     lsau=imap(quring,saturn);
     7703     @h=imap(quring,@h);
     7704
     7705     kill quring;
     7706
     7707     @n2=size(quprimary);
     7708     @n3=@n2;
     7709
     7710     for(@n1=1;@n1<=size(collectprimary)/2;@n1++)
     7711     {
     7712        if(deg(collectprimary[2*@n1][1])>0)
     7713        {
     7714           @n2++;
     7715           quprimary[@n2]=collectprimary[2*@n1-1];
     7716           lnew[@n2]=lsau[2*@n1-1];
     7717           @n2++;
     7718           lnew[@n2]=lsau[2*@n1];
     7719           quprimary[@n2]=collectprimary[2*@n1];
     7720           if(abspri)
     7721           {
     7722              absprimary[@n2/2]=absprimarytmp[@n1];
     7723              abskeep[@n2/2]=abskeeptmp[@n1];
     7724           }
     7725        }
     7726     }
     7727
     7728     //here the intersection with the polynomialring
     7729     //mentioned above is really computed
     7730     for(@n=@n3/2+1;@n<=@n2/2;@n++)
     7731     {
     7732        if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n]))
     7733        {
     7734           quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7735           quprimary[2*@n]=quprimary[2*@n-1];
     7736        }
     7737        else
     7738        {
     7739           if(@wr==0)
     7740           {
     7741              quprimary[2*@n-1]=sat2(quprimary[2*@n-1],lnew[2*@n-1])[1];
     7742           }
     7743           quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1];
     7744        }
     7745     }
     7746
     7747     return(quprimary, absprimary, abskeep, ser, @h);
     7748}
     7749
     7750
     7751////////////////////////////////////////////////////////////////////////////
     7752
     7753
     7754
     7755
     7756///////////////////////////////////////////////////////////////////////////////
     7757// Based on minAssGTZ
     7758
     7759proc newMinAss(ideal i,list #)
     7760"USAGE:   newMinAss(i); i ideal
     7761          newMinAss(i,1); i ideal  does not use the factorizing Groebner
     7762RETURN:  a list, the minimal associated prime ideals of i.
     7763NOTE:    Designed for characteristic 0, works also in char k > 0 based
     7764         on an algorithm of Yokoyama
     7765EXAMPLE: example newMinAss; shows an example
     7766"
     7767{
     7768   string algorithm = "SL";
     7769   if(ord_test(basering) != 1)
     7770   {
     7771      ERROR(
     7772      "// Not implemented for this ordering, please change to global ordering."
     7773      );
     7774   }
     7775   if(minpoly!=0)
     7776   {
     7777      return(algeDeco(i, 2));
     7778   }
     7779   if(size(#)==0)
     7780   {
     7781      return(minAssPrimes(i, "facstd", algorithm));
     7782   }
     7783   return(minAssPrimes(i, algorithm));
     7784}
     7785example
     7786{ "EXAMPLE:";  echo = 2;
     7787   ring  r = 0, (x, y, z), dp;
     7788   poly  p = z2 + 1;
     7789   poly  q = z3 + 2;
     7790   ideal i = p * q^2, y - z2;
     7791   list pr = newMinAss(i);
     7792   pr;
     7793}
     7794
     7795
     7796
     7797
     7798
     7799
     7800
     7801
     7802///////////////////////////////////////////////////////////////////////////////
     7803//
     7804// Computes the minimal associated primes of I via the new algorithm,
     7805// using primary decomposition in the zero dimensional case.
     7806// For reduction to the zerodimensional case, it uses the procedure
     7807// decomp of primdec.lib, with some modifications to avoid the recursion.
     7808//
     7809
     7810proc newMinAssSL(ideal I)
     7811// Input = I, ideal
     7812// Output = primaryDec where primaryDec is the list of the minimal
     7813// associated primes and the primary components corresponding to these primes.
     7814
     7815{
     7816  ideal P = 1;
     7817  list pd = list();
     7818  int k;
     7819  int stop = 0;
     7820  list primaryDec = list();
     7821
     7822  while (stop == 0) {
     7823    // Debug
     7824    dbprint(printlevel - voice, "// We call newMinAssSLIteration to find new prime ideals!");
     7825    pd = newMinAssSLIteration(I, P);
     7826    // Debug
     7827    dbprint(printlevel - voice, "// Output of newMinAssSLIteration:");
     7828    dbprint(printlevel - voice, pd);
     7829    if (size(pd[1]) > 0) {
     7830      primaryDec = primaryDec + pd[1];
     7831      // Debug
     7832      dbprint(printlevel - voice, "// We intersect the prime ideals obtained.");
     7833      /*
     7834      for (k = 1; k <= size(pd); k++) {
     7835        P = intersect(P, pd[k]);
     7836      }
     7837      */
     7838      P = intersect(P, pd[2]);
     7839      // Debug
     7840      dbprint(printlevel - voice, "// Intersection finished.");
     7841    } else {
     7842      stop = 1;
     7843    }
     7844  }
     7845  //return(insert(primaryDec, P));
     7846  // Returns only the primary components, not the radical.
     7847  return(primaryDec);
     7848}
     7849
     7850///////////////////////////////////////////////////////////////////////////////
     7851// Given an ideal I and an ideal P (intersection of some minimal prime ideals
     7852// associated to I), it calculates new minimal prime ideals associated to I
     7853// which were not used to calculate P.
     7854// This version uses Primary Decomposition in the zerodimensional case.
     7855proc newMinAssSLIteration(ideal I, ideal P);
     7856{
     7857  int k = 1;
     7858  int good  = 0;
     7859  list primaryDec = list();
     7860  // Debug
     7861  dbprint (printlevel-voice, "// We search for an element in P - sqrt(I).");
     7862  while ((k <= size(P)) and (good == 0)) {
     7863    good = 1 - rad_con(P[k], I);
     7864    k++;
     7865  }
     7866  k--;
     7867  if (good == 0) {
     7868    // Debug
     7869    dbprint (printlevel - voice, "// No element was found, P = sqrt(I).");
     7870    return (list(primaryDec, ideal(0)));
     7871  }
     7872  // Debug
     7873  dbprint (printlevel - voice, "// We found h = ", P[k]);
     7874  dbprint (printlevel - voice, "// We calculate the saturation of I with respect to the element just founded.");
     7875  ideal J = sat(I, P[k])[1];
     7876
     7877  // Uses decomp from primdec, modified to avoid the recursion.
     7878  // Debug
     7879  dbprint(printlevel - voice, "// We do the reduction to the zerodimensional case, via decomp.");
     7880
     7881  primaryDec = newDecompStep(J, "oneIndep", "intersect", 2);
     7882  // Debug
     7883  dbprint(printlevel - voice, "// Proc decomp has found", size(primaryDec) / 2, "new primary components.");
     7884
     7885  return(primaryDec);
     7886}
     7887
     7888
     7889
     7890///////////////////////////////////////////////////////////////////////////////////
     7891// Based on maxIndependSet
     7892// Added list # as parameter
     7893// If the first element of # is 0, the output is only 1 max indep set.
     7894// If no list is specified or #[1] = 1, the output is all the max indep set of the
     7895// leading terms ideal. This is the original output of maxIndependSet
     7896
     7897// The ordering given in the output has been changed to block dp instead of lp.
     7898
     7899proc newMaxIndependSetLp(ideal j, list #)
     7900// #[1]=0 -->     computes only one max indep set
     7901"USAGE:   maxIndependentSet(i); i ideal
     7902RETURN:  list = #1. new varstring with the maximal independent set at the end,
     7903                #2. ordstring with the corresponding block ordering,
     7904                #3. the number of independent variables
     7905                // SL #3 explanation modified,
     7906                // it was: integer where the independent set starts in the varstring
     7907NOTE:
     7908EXAMPLE: example maxIndependentSetLp; shows an example
     7909"
     7910{
     7911   int n, k, di;
     7912   list resu, hilf;
     7913   string var1, var2;
     7914   list v = indepSet(j, 0);
     7915
     7916   // SL 2006.04.21 1 Lines modified to use only one independent Set
     7917   string indepOption;
     7918   if (size(#) > 0) {
     7919     indepOption = #[1];
     7920   } else {
     7921     indepOption = "allIndep";
     7922   }
     7923
     7924   int nMax;
     7925   if (indepOption == "allIndep") {
     7926      nMax = size(v);
     7927   } else {
     7928    nMax = 1;
     7929   }
     7930
     7931   for(n = 1; n <= nMax; n++)
     7932   // SL 2006.04.21 2
     7933   {
     7934      di = 0;
     7935      var1 = "";
     7936      var2 = "";
     7937      for(k = 1; k <= size(v[n]); k++)
     7938      {
     7939         if(v[n][k] != 0)
     7940         {
     7941            di++;
     7942            var2 = var2 + "var(" + string(k) + "), ";
     7943         }
     7944         else
     7945         {
     7946            var1 = var1 + "var(" + string(k) + "), ";
     7947         }
     7948      }
     7949      if(di > 0)
     7950      {
     7951         var1 = var1 + var2;
     7952         var1 = var1[1..size(var1) - 2];       // The "- 2" removes the trailer comma
     7953         hilf[1] = var1;
     7954         // SL 2006.21.04 1 The order is now block dp instead of lp
     7955         //hilf[2] = "dp(" + string(nvars(basering) - di) + "), dp(" + string(di) + ")";
     7956         // SL 2006.21.04 2
     7957         // For decomp, lp ordering is needed. Nothing is changed.
     7958         hilf[2] = "lp";
     7959         hilf[3] = di;
     7960         resu[n] = hilf;
     7961      }
     7962      else
     7963      {
     7964         resu[n] = varstr(basering), ordstr(basering), 0;
     7965      }
     7966   }
     7967   return(resu);
     7968}
     7969example
     7970{ "EXAMPLE:"; echo = 2;
     7971   ring s1 = (0, x, y), (a, b, c, d, e, f, g), lp;
     7972   ideal i = ea - fbg, fa + be, ec - fdg, fc + de;
     7973   i = std(i);
     7974   list l = newMaxIndependSetLp(i);
     7975   l;
     7976   i = i, g;
     7977   l = newMaxIndependSetLp(i);
     7978   l;
     7979
     7980   ring s = 0, (x, y, z), lp;
     7981   ideal i = z, yx;
     7982   list l = newMaxIndependSetLp(i);
     7983   l;
     7984}
     7985
     7986
     7987///////////////////////////////////////////////////////////////////////////////
     7988
     7989proc newZero_decomp (ideal j, ideal ser, int @wr, list #)
     7990"USAGE:   newZero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
     7991         (@wr=0 for primary decomposition, @wr=1 for computaion of associated
     7992         primes)
     7993         if #[1] = "nest", then #[2] indicates the nest level (number of recursive calls)
     7994         When the nest level is high it indicates that the computation is difficult,
     7995         and different methods are applied.
     7996RETURN:  list = list of primary ideals and their radicals (at even positions
     7997         in the list) if the input is zero-dimensional and a standardbases
     7998         with respect to lex-ordering
     7999         If ser!=(0) and ser is contained in j or if j is not zero-dimen-
     8000         sional then ideal(1),ideal(1) is returned
     8001NOTE:    Algorithm of Gianni/Trager/Zacharias
     8002EXAMPLE: example newZero_decomp; shows an example
     8003"
     8004{
     8005  def   @P = basering;
     8006  int uytrewq;
     8007  int nva = nvars(basering);
     8008  int @k,@s,@n,@k1,zz;
     8009  list primary,lres0,lres1,act,@lh,@wh;
     8010  map phi,psi,phi1,psi1;
     8011  ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1;
     8012  intvec @vh,@hilb;
     8013  string @ri;
     8014  poly @f;
     8015
     8016  // Debug
     8017  dbprint(printlevel - voice, "proc newZero_decomp");
     8018
     8019
     8020  if (dim(j)>0)
     8021  {
     8022     primary[1]=ideal(1);
     8023     primary[2]=ideal(1);
     8024     return(primary);
     8025  }
     8026  j=interred(j);
     8027
     8028  attrib(j,"isSB",1);
     8029
     8030  int nestLevel = 0;
     8031  if (size(#) > 0){
     8032    if (typeof(#[1]) == "string") {
     8033      if (#[1] == "nest") {
     8034        nestLevel = #[2];
     8035      }
     8036      # = list();
     8037    }
     8038  }
     8039
     8040  if(vdim(j)==deg(j[1]))
     8041  {
     8042     act=factor(j[1]);
     8043     for(@k=1;@k<=size(act[1]);@k++)
     8044     {
     8045       @qh=j;
     8046       if(@wr==0)
     8047       {
     8048          @qh[1]=act[1][@k]^act[2][@k];
     8049       }
     8050       else
     8051       {
     8052          @qh[1]=act[1][@k];
     8053       }
     8054       primary[2*@k-1]=interred(@qh);
     8055       @qh=j;
     8056       @qh[1]=act[1][@k];
     8057       primary[2*@k]=interred(@qh);
     8058       attrib( primary[2*@k-1],"isSB",1);
     8059
     8060       if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0))
     8061       {
     8062          primary[2*@k-1]=ideal(1);
     8063          primary[2*@k]=ideal(1);
     8064       }
     8065     }
     8066     return(primary);
     8067  }
     8068
     8069  if(homog(j)==1)
     8070  {
     8071     primary[1]=j;
     8072     if((size(ser)>0)&&(size(reduce(ser,j,1))==0))
     8073     {
     8074          primary[1]=ideal(1);
     8075          primary[2]=ideal(1);
     8076          return(primary);
     8077     }
     8078     if(dim(j)==-1)
     8079     {
     8080        primary[1]=ideal(1);
     8081        primary[2]=ideal(1);
     8082     }
     8083     else
     8084     {
     8085        primary[2]=maxideal(1);
     8086     }
     8087     return(primary);
     8088  }
     8089
     8090//the first element in the standardbase is factorized
     8091  if(deg(j[1])>0)
     8092  {
     8093    act=factor(j[1]);
     8094    testFactor(act,j[1]);
     8095  }
     8096  else
     8097  {
     8098     primary[1]=ideal(1);
     8099     primary[2]=ideal(1);
     8100     return(primary);
     8101  }
     8102
     8103//with the factors new ideals (hopefully the primary decomposition)
     8104//are created
     8105  if(size(act[1])>1)
     8106  {
     8107     if(size(#)>1)
     8108     {
     8109        primary[1]=ideal(1);
     8110        primary[2]=ideal(1);
     8111        primary[3]=ideal(1);
     8112        primary[4]=ideal(1);
     8113        return(primary);
     8114     }
     8115     for(@k=1;@k<=size(act[1]);@k++)
     8116     {
     8117        if(@wr==0)
     8118        {
     8119           primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
     8120
     8121        }
     8122        else
     8123        {
     8124           primary[2*@k-1]=std(j,act[1][@k]);
     8125        }
     8126        if((act[2][@k]==1)&&(vdim(primary[2*@k-1])==deg(act[1][@k])))
     8127        {
     8128           primary[2*@k]   = primary[2*@k-1];
     8129        }
     8130        else
     8131        {
     8132
     8133           primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
     8134
     8135        }
     8136     }
     8137  }
     8138  else
     8139  {
     8140     primary[1]=j;
     8141     if((size(#)>0)&&(act[2][1]>1))
     8142     {
     8143        act[2]=1;
     8144        primary[1]=std(primary[1],act[1][1]);
     8145     }
     8146     if(@wr!=0)
     8147     {
     8148        primary[1]=std(j,act[1][1]);
     8149     }
     8150     if((act[2][1]==1)&&(vdim(primary[1])==deg(act[1][1])))
     8151     {
     8152        primary[2]=primary[1];
     8153     }
     8154     else
     8155     {
     8156        primary[2]=primaryTest(primary[1],act[1][1]);
     8157     }
     8158  }
     8159
     8160  if(size(#)==0)
     8161  {
     8162     primary=splitPrimary(primary,ser,@wr,act);
     8163  }
     8164
     8165  if((voice>=6)&&(char(basering)<=181))
     8166  {
     8167     primary=splitCharp(primary);
     8168  }
     8169
     8170  if((@wr==2)&&(npars(basering)>0)&&(voice>=6)&&(char(basering)>0))
     8171  {
     8172  //the prime decomposition of Yokoyama in characteristic p
     8173     list ke,ek;
     8174     @k=0;
     8175     while(@k<size(primary)/2)
     8176     {
     8177        @k++;
     8178        if(size(primary[2*@k])==0)
     8179        {
     8180           ek=insepDecomp(primary[2*@k-1]);
     8181           primary=delete(primary,2*@k);
     8182           primary=delete(primary,2*@k-1);
     8183           @k--;
     8184        }
     8185        ke=ke+ek;
     8186     }
     8187     for(@k=1;@k<=size(ke);@k++)
     8188     {
     8189        primary[size(primary)+1]=ke[@k];
     8190        primary[size(primary)+1]=ke[@k];
     8191     }
     8192  }
     8193
     8194  if(nestLevel > 1){primary=extF(primary)}
     8195
     8196//test whether all ideals in the decomposition are primary and
     8197//in general position
     8198//if not after a random coordinate transformation of the last
     8199//variable the corresponding ideal is decomposed again.
     8200  if((npars(basering)>0)&&(nestLevel > 1))
     8201  {
     8202     poly randp;
     8203     for(zz=1;zz<nvars(basering);zz++)
     8204     {
     8205        randp=randp
     8206              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(zz);
     8207     }
     8208     randp=randp+var(nvars(basering));
     8209  }
     8210  @k=0;
     8211  while(@k<(size(primary)/2))
     8212  {
     8213    @k++;
     8214    if (size(primary[2*@k])==0)
     8215    {
     8216       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
     8217       {
     8218          attrib(primary[2*@k-1],"isSB",1);
     8219          if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
     8220          {
     8221             primary[2*@k]=primary[2*@k-1];
     8222          }
     8223       }
     8224    }
     8225  }
     8226
     8227  @k=0;
     8228  ideal keep;
     8229  while(@k<(size(primary)/2))
     8230  {
     8231    @k++;
     8232    if (size(primary[2*@k])==0)
     8233    {
     8234
     8235       jmap=randomLast(100);
     8236       jmap1=maxideal(1);
     8237       jmap2=maxideal(1);
     8238       @qht=primary[2*@k-1];
     8239       if((npars(basering)>0)&&(nestLevel > 1))
     8240       {
     8241          jmap[size(jmap)]=randp;
     8242       }
     8243
     8244       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
     8245       {
     8246          if(deg(lead(primary[2*@k-1][@n]))==1)
     8247          {
     8248             for(zz=1;zz<=nva;zz++)
     8249             {
     8250                if(lead(primary[2*@k-1][@n])/var(zz)!=0)
     8251                {
     8252                 jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n]
     8253                   +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]);
     8254                    jmap2[zz]=primary[2*@k-1][@n];
     8255                    @qht[@n]=var(zz);
     8256
     8257                }
     8258             }
     8259             jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0);
     8260          }
     8261       }
     8262       if(size(subst(jmap[nva],var(1),0)-var(nva))!=0)
     8263       {
     8264         // jmap[nva]=subst(jmap[nva],var(1),0);
     8265         //hier geaendert +untersuchen!!!!!!!!!!!!!!
     8266       }
     8267       phi1=@P,jmap1;
     8268       phi=@P,jmap;
     8269       for(@n=1;@n<=nva;@n++)
     8270       {
     8271          jmap[@n]=-(jmap[@n]-2*var(@n));
     8272       }
     8273       psi=@P,jmap;
     8274       psi1=@P,jmap2;
     8275       @qh=phi(@qht);
     8276
     8277//=================== the new part ============================
     8278
     8279       @qh=groebner(@qh);
     8280
     8281//=============================================================
     8282//       if(npars(@P)>0)
     8283//       {
     8284//          @ri= "ring @Phelp ="
     8285//                  +string(char(@P))+",
     8286//                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
     8287//       }
     8288//       else
     8289//       {
     8290//          @ri= "ring @Phelp ="
     8291//                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
     8292//       }
     8293//       execute(@ri);
     8294//       ideal @qh=homog(imap(@P,@qht),@t);
     8295//
     8296//       ideal @qh1=std(@qh);
     8297//       @hilb=hilb(@qh1,1);
     8298//       @ri= "ring @Phelp1 ="
     8299//                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
     8300//       execute(@ri);
     8301//       ideal @qh=homog(imap(@P,@qh),@t);
     8302//       kill @Phelp;
     8303//       @qh=std(@qh,@hilb);
     8304//       @qh=subst(@qh,@t,1);
     8305//       setring @P;
     8306//       @qh=imap(@Phelp1,@qh);
     8307//       kill @Phelp1;
     8308//       @qh=clearSB(@qh);
     8309//       attrib(@qh,"isSB",1);
     8310//=============================================================
     8311
     8312       ser1=phi1(ser);
     8313       @lh=newZero_decomp (@qh,phi(ser1),@wr, list("nest", nestLevel + 1));
     8314
     8315       kill lres0;
     8316       list lres0;
     8317       if(size(@lh)==2)
     8318       {
     8319          helpprim=@lh[2];
     8320          lres0[1]=primary[2*@k-1];
     8321          ser1=psi(helpprim);
     8322          lres0[2]=psi1(ser1);
     8323          if(size(reduce(lres0[2],lres0[1],1))==0)
     8324          {
     8325             primary[2*@k]=primary[2*@k-1];
     8326             continue;
     8327          }
     8328       }
     8329       else
     8330       {
     8331
     8332          lres1=psi(@lh);
     8333          lres0=psi1(lres1);
     8334       }
     8335
     8336//=================== the new part ============================
     8337
     8338       primary=delete(primary,2*@k-1);
     8339       primary=delete(primary,2*@k-1);
     8340       @k--;
     8341       if(size(lres0)==2)
     8342       {
     8343          lres0[2]=groebner(lres0[2]);
     8344       }
     8345       else
     8346       {
     8347          for(@n=1;@n<=size(lres0)/2;@n++)
     8348          {
     8349             if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1)
     8350             {
     8351                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
     8352                lres0[2*@n]=lres0[2*@n-1];
     8353                attrib(lres0[2*@n],"isSB",1);
     8354             }
     8355             else
     8356             {
     8357                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
     8358                lres0[2*@n]=groebner(lres0[2*@n]);
     8359             }
     8360          }
     8361       }
     8362       primary=primary+lres0;
     8363
     8364//=============================================================
     8365//       if(npars(@P)>0)
     8366//       {
     8367//          @ri= "ring @Phelp ="
     8368//                  +string(char(@P))+",
     8369//                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
     8370//       }
     8371//       else
     8372//       {
     8373//          @ri= "ring @Phelp ="
     8374//                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
     8375//       }
     8376//       execute(@ri);
     8377//       list @lvec;
     8378//       list @lr=imap(@P,lres0);
     8379//       ideal @lr1;
     8380//
     8381//       if(size(@lr)==2)
     8382//       {
     8383//          @lr[2]=homog(@lr[2],@t);
     8384//          @lr1=std(@lr[2]);
     8385//          @lvec[2]=hilb(@lr1,1);
     8386//       }
     8387//       else
     8388//       {
     8389//          for(@n=1;@n<=size(@lr)/2;@n++)
     8390//          {
     8391//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
     8392//             {
     8393//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
     8394//                @lr1=std(@lr[2*@n-1]);
     8395//                @lvec[2*@n-1]=hilb(@lr1,1);
     8396//                @lvec[2*@n]=@lvec[2*@n-1];
     8397//             }
     8398//             else
     8399//             {
     8400//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
     8401//                @lr1=std(@lr[2*@n-1]);
     8402//                @lvec[2*@n-1]=hilb(@lr1,1);
     8403//                @lr[2*@n]=homog(@lr[2*@n],@t);
     8404//                @lr1=std(@lr[2*@n]);
     8405//                @lvec[2*@n]=hilb(@lr1,1);
     8406//
     8407//             }
     8408//         }
     8409//       }
     8410//       @ri= "ring @Phelp1 ="
     8411//                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
     8412//       execute(@ri);
     8413//       list @lr=imap(@Phelp,@lr);
     8414//
     8415//       kill @Phelp;
     8416//       if(size(@lr)==2)
     8417//      {
     8418//          @lr[2]=std(@lr[2],@lvec[2]);
     8419//          @lr[2]=subst(@lr[2],@t,1);
     8420//
     8421//       }
     8422//       else
     8423//       {
     8424//          for(@n=1;@n<=size(@lr)/2;@n++)
     8425//          {
     8426//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
     8427//             {
     8428//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
     8429//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
     8430//                @lr[2*@n]=@lr[2*@n-1];
     8431//                attrib(@lr[2*@n],"isSB",1);
     8432//             }
     8433//             else
     8434//             {
     8435//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
     8436//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
     8437//                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
     8438//                @lr[2*@n]=subst(@lr[2*@n],@t,1);
     8439//             }
     8440//          }
     8441//       }
     8442//       kill @lvec;
     8443//       setring @P;
     8444//       lres0=imap(@Phelp1,@lr);
     8445//       kill @Phelp1;
     8446//       for(@n=1;@n<=size(lres0);@n++)
     8447//       {
     8448//          lres0[@n]=clearSB(lres0[@n]);
     8449//          attrib(lres0[@n],"isSB",1);
     8450//       }
     8451//
     8452//       primary[2*@k-1]=lres0[1];
     8453//       primary[2*@k]=lres0[2];
     8454//       @s=size(primary)/2;
     8455//       for(@n=1;@n<=size(lres0)/2-1;@n++)
     8456//       {
     8457//         primary[2*@s+2*@n-1]=lres0[2*@n+1];
     8458//         primary[2*@s+2*@n]=lres0[2*@n+2];
     8459//       }
     8460//       @k--;
     8461//=============================================================
     8462     }
     8463  }
     8464  return(primary);
     8465}
     8466example
     8467{ "EXAMPLE:"; echo = 2;
     8468   ring  r = 0,(x,y,z),lp;
     8469   poly  p = z2+1;
     8470   poly  q = z4+2;
     8471   ideal i = p^2*q^3,(y-z3)^3,(x-yz+z4)^4;
     8472   i=std(i);
     8473   list  pr= newZero_decomp(i,ideal(0),0);
     8474   pr;
     8475}
     8476///////////////////////////////////////////////////////////////////////////////
     8477
    64578478////////////////////////////////////////////////////////////////////////////
    64588479/*
Note: See TracChangeset for help on using the changeset viewer.