Changeset d88470 in git for Singular/LIB/primdec.lib


Ignore:
Timestamp:
Nov 13, 2006, 1:57:33 PM (17 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
1aecaec05d17a34b8ff6a44d71e1684cc74281e4
Parents:
11a4417493a03a79effa983f174d085bb0c77f5e
Message:
*hannes: format, removed ord_test


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/primdec.lib

    r11a441 rd88470  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: primdec.lib,v 1.127 2006-11-08 13:45:02 Singular Exp $";
     2version="$Id: primdec.lib,v 1.128 2006-11-13 12:57:33 Singular Exp $";
    33category="Commutative Algebra";
    44info="
    … …  
    23242324"
    23252325{
    2326   if(ord_test(basering)!=1)
     2326  if(attrib(basering,"global")!=1)
    23272327  {
    23282328      ERROR(
    … …  
    24302430"
    24312431{
    2432   if(ord_test(basering)!=1)
     2432  if(attrib(basering,"global")!=1)
    24332433  {
    24342434      ERROR(
    … …  
    37993799"
    38003800{
    3801    if(ord_test(basering)!=1)
     3801   if(attrib(basering,"global")!=1)
    38023802   {
    38033803      ERROR(
    … …  
    50805080"
    50815081{
    5082    if(ord_test(basering)!=1)
     5082   if(attrib(basering,"global")!=1)
    50835083   {
    50845084      ERROR(
    … …  
    51325132  }
    51335133
    5134    if(ord_test(basering)!=1)
     5134   if(attrib(basering,"global")!=1)
    51355135   {
    51365136      ERROR(
    … …  
    52315231"
    52325232{
    5233    if(ord_test(basering)!=1)
     5233   if(attrib(basering,"global")!=1)
    52345234   {
    52355235      ERROR(
    … …  
    53085308   }
    53095309
    5310    if(ord_test(basering)!=1)
     5310   if(attrib(basering,"global")!=1)
    53115311   {
    53125312      ERROR(
    … …  
    53455345"
    53465346{
    5347    if(ord_test(basering)!=1)
     5347   if(attrib(basering,"global")!=1)
    53485348   {
    53495349      ERROR(
    … …  
    53745374"
    53755375{
    5376    if(ord_test(basering)!=1)
    5377    {
    5378       ERROR(
    5379       "// Not implemented for this ordering, please change to global ordering."
    5380       );
    5381    }
    5382    return(radical(i, 1));
     5376  if(attrib(basering,"global")!=1)
     5377  {
     5378     ERROR(
     5379     "// Not implemented for this ordering, please change to global ordering."
     5380     );
     5381  }
     5382  return(radical(i, 1));
    53835383}
    53845384example
    … …  
    54085408"
    54095409{
    5410    dbprint(printlevel - voice, "Radical, version 2006.05.08");
    5411    if(ord_test(basering)!=1)
    5412    {
    5413       ERROR(
    5414       "// Not implemented for this ordering, please change to global ordering."
    5415       );
    5416    }
    5417    if(size(i) == 0){return(ideal(0));}
    5418    int j;
    5419    def P0 = basering;
    5420    list Pl=ringlist(P0);
    5421    intvec dp_w;
    5422    for(j=nvars(P0);j>0;j--) {dp_w[j]=1;}
    5423    Pl[3]=list(list("dp",dp_w),list("C",0));
    5424    def @P=ring(Pl);
    5425    setring @P;
    5426    ideal i=imap(P0,i);
    5427 
    5428    int il;
    5429    string algorithm;
    5430    int useFac;
    5431 
    5432    // Set input parameters
    5433    algorithm = "SL";                                 // Default: SL algorithm
    5434    il = 0;                                           // Default: Full radical (not only equiRadical)
    5435    if (homog(i) == 1) {   // Default: facStd is used, except if the ideal is homogeneous.
    5436       useFac = 0;
    5437    } else {
    5438       useFac = 1;
    5439    }
    5440    if(size(#) > 0){
    5441       int valid;
    5442       for(j = 1; j <= size(#); j++){
    5443          valid = 0;
    5444          if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
    5445          {
    5446             il = #[j];          // If il == 1, equiRadical is computed
    5447             valid = 1;
    5448          }
    5449          if(typeof(#[j]) == "string"){
    5450             if(#[j] == "KL") {
    5451                algorithm = "KL";
    5452                valid = 1;}
    5453             if(#[j] == "SL") {
    5454                algorithm = "SL";
    5455                valid = 1;}
    5456             if(#[j] == "noFacstd") {
    5457                useFac = 0;
    5458                valid = 1;}
    5459             if(#[j] == "facstd") {
    5460                useFac = 1;
    5461                valid = 1;}
    5462             if(#[j] == "equiRad") {
    5463                il = 1;
    5464                valid = 1;}
    5465             if(#[j] == "fullRad") {
    5466                il = 0;
    5467                valid = 1;}
    5468          }
    5469          if(valid == 0)
    5470          {
    5471             dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
    5472          }
    5473       }
    5474    }
    5475 
    5476    ideal rad = 1;
    5477    intvec op = option(get);
    5478    list qr = simplifyIdeal(i);
    5479    map phi = @P, qr[2];
    5480 
    5481    option(redSB);
    5482    i = groebner(qr[1]);
    5483    option(set, op);
    5484    int di = dim(i);
    5485 
    5486    if(di == 0)
    5487    {
    5488       i = zeroRad(i, qr[1]);
    5489       i=interred(phi(i));
    5490       setring(P0);
    5491       i=imap(@P,i);
    5492       return(i);
    5493    }
    5494 
    5495    option(redSB);
    5496    list pr;
    5497    if(useFac == 1)
    5498    {
    5499       pr = facstd(i);
    5500    } else {
    5501       pr = i;
    5502    }
    5503    option(set, op);
    5504    int s = size(pr);
    5505    if(useFac == 1) {
    5506       dbprint(printlevel - voice, "Number of components returned by facstd: ", s);
    5507    }
    5508    for(j = 1; j <= s; j++)
    5509    {
    5510       attrib(pr[s + 1 - j], "isSB", 1);
    5511       if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il))
    5512       {
    5513          // SL Debug messages
    5514          dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
    5515          dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));
    5516 
    5517          if(algorithm == "KL") {
    5518             rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
    5519          }
    5520          if(algorithm == "SL") {
    5521             rad = intersect(rad, radicalSL(pr[s + 1 - j], il));
    5522          }
    5523       } else
    5524       {
    5525          // SL Debug
    5526          dbprint(printlevel-voice, "The radical of this component is not needed.");
    5527          dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))",
    5528                  size(reduce(rad, pr[s + 1 - j], 1)));
    5529          dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
    5530          dbprint(printlevel-voice, "il", il);
    5531       }
    5532    }
    5533    rad=interred(phi(rad));
    5534    setring(P0);
    5535    i=imap(@P,rad);
    5536    return(i);
     5410  dbprint(printlevel - voice, "Radical, version 2006.05.08");
     5411  if(attrib(basering,"global")!=1)
     5412  {
     5413     ERROR(
     5414     "// Not implemented for this ordering, please change to global ordering."
     5415     );
     5416  }
     5417  if(size(i) == 0){return(ideal(0));}
     5418  int j;
     5419  def P0 = basering;
     5420  list Pl=ringlist(P0);
     5421  intvec dp_w;
     5422  for(j=nvars(P0);j>0;j--) {dp_w[j]=1;}
     5423  Pl[3]=list(list("dp",dp_w),list("C",0));
     5424  def @P=ring(Pl);
     5425  setring @P;
     5426  ideal i=imap(P0,i);
     5427
     5428  int il;
     5429  string algorithm;
     5430  int useFac;
     5431
     5432  // Set input parameters
     5433  algorithm = "SL";                                 // Default: SL algorithm
     5434  il = 0;                                           // Default: Full radical (not only equiRadical)
     5435  if (homog(i) == 1)
     5436  {   // Default: facStd is used, except if the ideal is homogeneous.
     5437     useFac = 0;
     5438  } else {
     5439     useFac = 1;
     5440  }
     5441  if(size(#) > 0){
     5442    int valid;
     5443    for(j = 1; j <= size(#); j++)
     5444    {
     5445      valid = 0;
     5446      if((typeof(#[j]) == "int") or (typeof(#[j]) == "number"))
     5447      {
     5448        il = #[j];          // If il == 1, equiRadical is computed
     5449        valid = 1;
     5450      }
     5451      if(typeof(#[j]) == "string"){
     5452        if(#[j] == "KL") {
     5453          algorithm = "KL";
     5454          valid = 1;
     5455        }
     5456        if(#[j] == "SL") {
     5457          algorithm = "SL";
     5458          valid = 1;
     5459        }
     5460        if(#[j] == "noFacstd") {
     5461          useFac = 0;
     5462          valid = 1;}
     5463        if(#[j] == "facstd") {
     5464          useFac = 1;
     5465          valid = 1;}
     5466        if(#[j] == "equiRad") {
     5467          il = 1;
     5468          valid = 1;}
     5469        if(#[j] == "fullRad") {
     5470          il = 0;
     5471          valid = 1;}
     5472      }
     5473      if(valid == 0)
     5474      {
     5475        dbprint(1, "Warning! The following input parameter was not recognized:", #[j]);
     5476      }
     5477    }
     5478  }
     5479
     5480  ideal rad = 1;
     5481  intvec op = option(get);
     5482  list qr = simplifyIdeal(i);
     5483  map phi = @P, qr[2];
     5484
     5485  option(redSB);
     5486  i = groebner(qr[1]);
     5487  option(set, op);
     5488  int di = dim(i);
     5489
     5490  if(di == 0)
     5491  {
     5492    i = zeroRad(i, qr[1]);
     5493    i=interred(phi(i));
     5494    setring(P0);
     5495    i=imap(@P,i);
     5496    return(i);
     5497  }
     5498
     5499  option(redSB);
     5500  list pr;
     5501  if(useFac == 1)
     5502  {
     5503    pr = facstd(i);
     5504  } else {
     5505    pr = i;
     5506  }
     5507  option(set, op);
     5508  int s = size(pr);
     5509  if(useFac == 1) {
     5510    dbprint(printlevel - voice, "Number of components returned by facstd: ", s);
     5511  }
     5512  for(j = 1; j <= s; j++)
     5513  {
     5514    attrib(pr[s + 1 - j], "isSB", 1);
     5515    if((size(reduce(rad, pr[s + 1 - j], 1)) != 0) && ((dim(pr[s + 1 - j]) == di) || !il))
     5516    {
     5517      // SL Debug messages
     5518      dbprint(printlevel-voice, "We shall compute the radical of ", pr[s + 1 - j]);
     5519      dbprint(printlevel-voice, "The dimension is: ", dim(pr[s+1-j]));
     5520
     5521      if(algorithm == "KL")
     5522      {
     5523        rad = intersect(rad, radicalKL(pr[s + 1 - j], rad, il));
     5524      }
     5525      if(algorithm == "SL") {
     5526        rad = intersect(rad, radicalSL(pr[s + 1 - j], il));
     5527      }
     5528    }
     5529    else
     5530    {
     5531      // SL Debug
     5532      dbprint(printlevel-voice, "The radical of this component is not needed.");
     5533      dbprint(printlevel-voice, "size(reduce(rad, pr[s + 1 - j], 1))",
     5534              size(reduce(rad, pr[s + 1 - j], 1)));
     5535      dbprint(printlevel-voice, "dim(pr[s + 1 - j])", dim(pr[s + 1 - j]));
     5536      dbprint(printlevel-voice, "il", il);
     5537    }
     5538  }
     5539  rad=interred(phi(rad));
     5540  setring(P0);
     5541  i=imap(@P,rad);
     5542  return(i);
    55375543}
    55385544example
    … …  
    61216127"
    61226128{
    6123   if(ord_test(basering)!=1)
     6129  if(attrib(basering,"global")!=1)
    61246130  {
    61256131      ERROR(
    … …  
    61696175"
    61706176{
    6171   if(ord_test(basering)!=1)
     6177  if(attrib(basering,"global")!=1)
    61726178  {
    61736179      ERROR(
    … …  
    62376243"
    62386244{
    6239    if(ord_test(basering)!=1)
    6240    {
    6241       ERROR(
    6242       "// Not implemented for this ordering, please change to global ordering."
    6243       );
    6244    }
    6245    def R=basering;
    6246    poly q;
    6247    int j,time;
    6248    matrix m;
    6249    list re;
    6250    poly va=var(1);
    6251    ideal J=groebner(I);
    6252    ideal ba=kbase(J);
    6253    int d=vdim(J);
    6254    dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d));
     6245  if(attrib(basering,"global")!=1)
     6246  {
     6247    ERROR(
     6248    "// Not implemented for this ordering, please change to global ordering."
     6249    );
     6250  }
     6251  def R=basering;
     6252  poly q;
     6253  int j,time;
     6254  matrix m;
     6255  list re;
     6256  poly va=var(1);
     6257  ideal J=groebner(I);
     6258  ideal ba=kbase(J);
     6259  int d=vdim(J);
     6260  dbprint(printlevel-voice+2,"// multiplicity of ideal : "+ string(d));
    62556261//------ compute matrix of multiplication on R/I with generic element p -----
    6256    int e=nvars(basering);
    6257    poly p=randomLast(100)[e]+random(-50,50);     //the generic element
    6258    matrix n[d][d];
    6259    time = timer;
    6260    for(j=2;j<=e;j++)
    6261    {
    6262       va=va*var(j);
    6263    }
    6264    for(j=1;j<=d;j++)
    6265    {
    6266       q=reduce(p*ba[j],J);
    6267       m=coeffs(q,ba,va);
    6268       n[j,1..d]=m[1..d,1];
    6269    }
    6270    dbprint(printlevel-voice+2,
    6271       "// time for computing multiplication matrix (with generic element) : "+
    6272       string(timer-time));
     6262  int e=nvars(basering);
     6263  poly p=randomLast(100)[e]+random(-50,50);     //the generic element
     6264  matrix n[d][d];
     6265  time = timer;
     6266  for(j=2;j<=e;j++)
     6267  {
     6268    va=va*var(j);
     6269  }
     6270  for(j=1;j<=d;j++)
     6271  {
     6272    q=reduce(p*ba[j],J);
     6273    m=coeffs(q,ba,va);
     6274    n[j,1..d]=m[1..d,1];
     6275  }
     6276  dbprint(printlevel-voice+2,
     6277    "// time for computing multiplication matrix (with generic element) : "+
     6278    string(timer-time));
    62736279//---------------- compute characteristic polynomial of matrix --------------
    6274    execute("ring P1=("+charstr(R)+"),T,dp;");
    6275    matrix n=imap(R,n);
    6276    time = timer;
    6277    poly charpol=det(n-T*freemodule(d));
    6278    dbprint(printlevel-voice+2,"// time for computing char poly: "+
    6279            string(timer-time));
     6280  execute("ring P1=("+charstr(R)+"),T,dp;");
     6281  matrix n=imap(R,n);
     6282  time = timer;
     6283  poly charpol=det(n-T*freemodule(d));
     6284  dbprint(printlevel-voice+2,"// time for computing char poly: "+
     6285         string(timer-time));
    62806286//------------------- factorize characteristic polynomial -------------------
    62816287//check first if constant term of charpoly is != 0 (which is true for
    62826288//sufficiently generic element)
    6283    if(charpol[size(charpol)]!=0)
    6284    {
    6285      time = timer;
    6286      list fac=factor(charpol);
    6287      testFactor(fac,charpol);
    6288      dbprint(printlevel-voice+2,"// time for factorizing char poly: "+
    6289              string(timer-time));
    6290      int f=size(fac[1]);
     6289  if(charpol[size(charpol)]!=0)
     6290  {
     6291    time = timer;
     6292    list fac=factor(charpol);
     6293    testFactor(fac,charpol);
     6294    dbprint(printlevel-voice+2,"// time for factorizing char poly: "+
     6295            string(timer-time));
     6296    int f=size(fac[1]);
    62916297//--------------------------- the irreducible case --------------------------
    6292      if(f==1)
    6293      {
    6294        setring R;
    6295        re=I;
    6296        return(re);
    6297      }
     6298    if(f==1)
     6299    {
     6300      setring R;
     6301      re=I;
     6302      return(re);
     6303    }
    62986304//---------------------------- the reducible case ---------------------------
    62996305//if f_i are the irreducible factors of charpoly, mult=ri, then <I,g_i^ri>
    … …  
    63026308//Hence it is better to simultaneously reduce with I. For this we need a new
    63036309//ring.
    6304      execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);");
    6305      list rfac=imap(P1,fac);
    6306      intvec ov=option(get);;
    6307      option(redSB);
    6308      list re1;
    6309      ideal new = T-imap(R,p),imap(R,J);
    6310      attrib(new, "isSB",1);    //we know that new is a standard basis
    6311      for(j=1;j<=f;j++)
    6312      {
    6313         re1[j]=reduce(rfac[1][j]^rfac[2][j],new);
    6314      }
    6315      setring R;
    6316      re = imap(P,re1);
    6317      for(j=1;j<=f;j++)
    6318      {
    6319         J=I,re[j];
    6320         re[j]=interred(J);
    6321      }
    6322      option(set,ov);
    6323      return(re);
     6310    execute("ring P=("+charstr(R)+"),(T,"+varstr(R)+"),(dp(1),dp);");
     6311    list rfac=imap(P1,fac);
     6312    intvec ov=option(get);;
     6313    option(redSB);
     6314    list re1;
     6315    ideal new = T-imap(R,p),imap(R,J);
     6316    attrib(new, "isSB",1);    //we know that new is a standard basis
     6317    for(j=1;j<=f;j++)
     6318    {
     6319      re1[j]=reduce(rfac[1][j]^rfac[2][j],new);
     6320    }
     6321    setring R;
     6322    re = imap(P,re1);
     6323    for(j=1;j<=f;j++)
     6324    {
     6325      J=I,re[j];
     6326      re[j]=interred(J);
     6327    }
     6328    option(set,ov);
     6329    return(re);
    63246330  }
    63256331  else
    63266332//------------------- choice of generic element failed -------------------
    63276333  {
    6328      dbprint(printlevel-voice+2,"// try new generic element!");
    6329      setring R;
    6330      return(zerodec(I));
     6334    dbprint(printlevel-voice+2,"// try new generic element!");
     6335    setring R;
     6336    return(zerodec(I));
    63316337  }
    63326338}
Note: See TracChangeset for help on using the changeset viewer.