Changeset 24f458 in git


Ignore:
Timestamp:
Feb 12, 2004, 10:54:06 AM (20 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
56ae4f7d3475c8a359c53899f648c760bbf91fc7
Parents:
e5b6d23c31318491b9f16c63a82fe33cafc3e9bf
Message:
*hannes: from 2-0


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/primdec.lib

    re5b6d2 r24f458  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: primdec.lib,v 1.101 2002-01-21 09:30:10 Singular Exp $";
     2version="$Id: primdec.lib,v 1.102 2004-02-12 09:54:06 Singular Exp $";
    33category="Commutative Algebra";
    44info="
    … …  
    1515 @* The procedures are implemented to be used in characteristic 0.
    1616 @* They also work in positive characteristic >> 0.
    17  @* In small characteristic and for algebraic extensions, primdecGTZ and
    18     minAssGTZ may not terminate, while primdecSY and minAssChar may not give
    19     a complete decomposition.
     17 @* In small characteristic and for algebraic extensions, primdecGTZ
     18    may not terminate.
    2019    Algorithms for the computation of the radical based on the ideas of
    2120    Krick, Logar and Kemper (implementation by Gerhard Pfister).
    2221
    2322PROCEDURES:
     23 Ann(M);           annihilator of R^n/M, R=basering, M in R^n
    2424 primdecGTZ(I);    complete primary decomposition via Gianni,Trager,Zacharias
    2525 primdecSY(I...);  complete primary decomposition via Shimoyama-Yokoyama
    … …  
    4343LIB "inout.lib";
    4444LIB "matrix.lib";
     45LIB "triang.lib";
    4546///////////////////////////////////////////////////////////////////////////////
    4647//
    … …  
    144145///////////////////////////////////////////////////////////////////////////////
    145146
    146 static proc minSat(ideal inew, ideal h)
     147
     148proc minSat(ideal inew, ideal h)
    147149{
    148150   int i,k;
    … …  
    372374   list  l = factor(p);
    373375   l;
    374 
    375 }
    376 
    377 
     376}
    378377
    379378///////////////////////////////////////////////////////////////////////////////
    … …  
    424423}
    425424///////////////////////////////////////////////////////////////////////////////
    426 
    427425
    428426proc primaryTest (ideal i, poly p)
    … …  
    736734  }
    737735  def P=basering;
    738   int i,j,k,m,q,d;
     736  int i,j,k,m,q,d,o;
    739737  int n=nvars(basering);
    740738  ideal s,t,u,sact;
    741739  poly ni;
    742740  string minp,gnir,va;
    743   list sa,keep;
     741  list sa,keep,rp,keep1;
    744742  for(i=1;i<=size(l)/2;i++)
    745743  {
    … …  
    756754    if(size(l[2*i])==0)
    757755    {
    758       s=factorize(l[2*i-1][1],1);
     756      s=factorize(l[2*i-1][1],1);   //vermeiden!!!
    759757      t=l[2*i-1];
    760758      m=size(t);
    … …  
    766764        while(va[j]!=","){j--;}
    767765        va=va[1..j-1];
    768         gnir="ring RL=("+string(char(P))+","+string(var(n))+"),("+va+"),dp;";
     766        gnir="ring RL=("+string(char(P))+","+string(var(n))+"),("+va+"),lp;";
    769767        execute(gnir);
    770768        minpoly=leadcoef(imap(P,ni));
    771769        ideal act;
    772770        ideal t=imap(P,t);
     771
    773772        for(k=2;k<=m;k++)
    774773        {
    775774           act=factorize(t[k],1);
    776            if(size(act)>1)
    777            {
    778 
    779              break;
    780            }
     775           if(size(act)>1){break;}
    781776        }
    782777        setring P;
    783778        sact=imap(RL,act);
    784         kill RL;
     779
    785780        if(size(sact)>1)
    786781        {
    … …  
    790785           l[2*i]=primaryTest(sa[1],sa[1][1]);
    791786        }
     787        if((size(sact)==1)&&(m==2))
     788        {
     789           l[2*i]=l[2*i-1];
     790           attrib(l[2*i],"isSB",1);
     791        }
     792        if((size(sact)==1)&&(m>2))
     793        {
     794           setring RL;
     795           option(redSB);
     796           t=std(t);
     797
     798           list sp=zero_decomp(t,0,0);
     799
     800           setring P;
     801           rp=imap(RL,sp);
     802           for(o=1;o<=size(rp);o++)
     803           {
     804              rp[o]=interred(simplify(rp[o],1)+ideal(ni));
     805           }
     806           l[2*i-1]=rp[1];
     807           l[2*i]=rp[2];
     808           rp=delete(rp,1);
     809           rp=delete(rp,1);
     810           keep1=keep1+rp;
     811           option(noredSB);
     812        }
     813        kill RL;
    792814      }
    793815    }
    … …  
    804826    }
    805827  }
     828  l=l+keep1;
    806829  return(l);
    807830}
    … …  
    809832///////////////////////////////////////////////////////////////////////////////
    810833
    811 static proc zero_decomp (ideal j,ideal ser,int @wr,list #)
     834proc zero_decomp (ideal j,ideal ser,int @wr,list #)
    812835"USAGE:   zero_decomp(j,ser,@wr); j,ser ideals, @wr=0 or 1
    813836         (@wr=0 for primary decomposition, @wr=1 for computaion of associated
    … …  
    841864
    842865  attrib(j,"isSB",1);
     866
    843867  if(vdim(j)==deg(j[1]))
    844868  {
    … …  
    906930//with the factors new ideals (hopefully the primary decomposition)
    907931//are created
    908 
    909932  if(size(act[1])>1)
    910933  {
    … …  
    922945        {
    923946           primary[2*@k-1]=std(j,act[1][@k]^act[2][@k]);
     947
    924948        }
    925949        else
    … …  
    933957        else
    934958        {
     959
    935960           primary[2*@k]   = primaryTest(primary[2*@k-1],act[1][@k]);
     961
    936962        }
    937963     }
    … …  
    960986     primary=splitPrimary(primary,ser,@wr,act);
    961987  }
    962   primary=splitCharp(primary);
     988
     989  if((voice>=6)&&(char(basering)<=181))
     990  {
     991     primary=splitCharp(primary);
     992  }
     993
     994  if((@wr==2)&&(npars(basering)>0)&&(voice>=6)&&(char(basering)>0))
     995  {
     996  //the prime decomposition of Yokoyama in characteristic p
     997     list ke,ek;
     998     @k=0;
     999     while(@k<size(primary)/2)
     1000     {
     1001        @k++;
     1002        if(size(primary[2*@k])==0)
     1003        {
     1004           ek=insepDecomp(primary[2*@k-1]);
     1005           primary=delete(primary,2*@k);
     1006           primary=delete(primary,2*@k-1);
     1007           @k--;
     1008        }
     1009        ke=ke+ek;
     1010     }
     1011     for(@k=1;@k<=size(ke);@k++)
     1012     {
     1013        primary[size(primary)+1]=ke[@k];
     1014        primary[size(primary)+1]=ke[@k];
     1015     }
     1016  }
     1017
     1018  if(voice>=8){primary=extF(primary)};
     1019
    9631020//test whether all ideals in the decomposition are primary and
    9641021//in general position
    9651022//if not after a random coordinate transformation of the last
    9661023//variable the corresponding ideal is decomposed again.
    967 
     1024  if((npars(basering)>0)&&(voice>=8))
     1025  {
     1026     poly randp;
     1027     for(zz=1;zz<nvars(basering);zz++)
     1028     {
     1029        randp=randp
     1030              +(random(0,5)*par(1)^2+random(0,5)*par(1)+random(0,5))*var(zz);
     1031     }
     1032     randp=randp+var(nvars(basering));
     1033  }
    9681034  @k=0;
    9691035  while(@k<(size(primary)/2))
    … …  
    9741040       for(zz=1;zz<size(primary[2*@k-1])-1;zz++)
    9751041       {
     1042          attrib(primary[2*@k-1],"isSB",1);
    9761043          if(vdim(primary[2*@k-1])==deg(primary[2*@k-1][zz]))
    9771044          {
    … …  
    9941061       jmap2=maxideal(1);
    9951062       @qht=primary[2*@k-1];
     1063       if((npars(basering)>0)&&(voice>=10))
     1064       {
     1065          jmap[size(jmap)]=randp;
     1066       }
    9961067
    9971068       for(@n=2;@n<=size(primary[2*@k-1]);@n++)
    … …  
    10271098       psi1=@P,jmap2;
    10281099       @qh=phi(@qht);
    1029        if(npars(@P)>0)
    1030        {
    1031           @ri= "ring @Phelp ="
    1032                   +string(char(@P))+",
    1033                    ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
    1034        }
    1035        else
    1036        {
    1037           @ri= "ring @Phelp ="
    1038                   +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
    1039        }
    1040        execute(@ri);
    1041        ideal @qh=homog(imap(@P,@qht),@t);
    1042 
    1043        ideal @qh1=std(@qh);
    1044        @hilb=hilb(@qh1,1);
    1045        @ri= "ring @Phelp1 ="
    1046                   +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
    1047        execute(@ri);
    1048        ideal @qh=homog(imap(@P,@qh),@t);
    1049        kill @Phelp;
    1050        @qh=std(@qh,@hilb);
    1051        @qh=subst(@qh,@t,1);
    1052        setring @P;
    1053        @qh=imap(@Phelp1,@qh);
    1054        kill @Phelp1;
    1055        @qh=clearSB(@qh);
    1056        attrib(@qh,"isSB",1);
     1100
     1101//=================== the new part ============================
     1102
     1103       @qh=groebner(@qh);
     1104
     1105//=============================================================
     1106//       if(npars(@P)>0)
     1107//       {
     1108//          @ri= "ring @Phelp ="
     1109//                  +string(char(@P))+",
     1110//                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
     1111//       }
     1112//       else
     1113//       {
     1114//          @ri= "ring @Phelp ="
     1115//                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
     1116//       }
     1117//       execute(@ri);
     1118//       ideal @qh=homog(imap(@P,@qht),@t);
     1119//
     1120//       ideal @qh1=std(@qh);
     1121//       @hilb=hilb(@qh1,1);
     1122//       @ri= "ring @Phelp1 ="
     1123//                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
     1124//       execute(@ri);
     1125//       ideal @qh=homog(imap(@P,@qh),@t);
     1126//       kill @Phelp;
     1127//       @qh=std(@qh,@hilb);
     1128//       @qh=subst(@qh,@t,1);
     1129//       setring @P;
     1130//       @qh=imap(@Phelp1,@qh);
     1131//       kill @Phelp1;
     1132//       @qh=clearSB(@qh);
     1133//       attrib(@qh,"isSB",1);
     1134//=============================================================
     1135
    10571136       ser1=phi1(ser);
     1137
    10581138       @lh=zero_decomp (@qh,phi(ser1),@wr);
    1059 //       @lh=zero_decomp (@qh,psi(ser),@wr);
     1139
    10601140
    10611141       kill lres0;
    … …  
    10751155       else
    10761156       {
    1077           //act=factor(@qh[1]);
    1078           //if(2*size(act[1])==size(@lh))
    1079           //{
    1080 
    1081          //   for(@n=1;@n<=size(act[1]);@n++)
    1082          //    {
    1083          //       @f=act[1][@n]^act[2][@n];
    1084          //       ser1=psi(@f);
    1085          //       lres0[2*@n-1]=interred(primary[2*@k-1]+psi1(ser1));
    1086          //       helpprim=@lh[2*@n];
    1087          //      ser1=psi(helpprim);
    1088          //       lres0[2*@n]=psi1(ser1);
    1089          //    }
    1090          // }
    1091          // else
    1092          // {
    1093              lres1=psi(@lh);
    1094              lres0=psi1(lres1);
    1095          //}
     1157
     1158          lres1=psi(@lh);
     1159          lres0=psi1(lres1);
    10961160       }
    1097        if(npars(@P)>0)
     1161
     1162//=================== the new part ============================
     1163
     1164       primary=delete(primary,2*@k-1);
     1165       primary=delete(primary,2*@k-1);
     1166       @k--;
     1167       if(size(lres0)==2)
    10981168       {
    1099           @ri= "ring @Phelp ="
    1100                   +string(char(@P))+",
    1101                    ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
     1169          lres0[2]=groebner(lres0[2]);
    11021170       }
    11031171       else
    11041172       {
    1105           @ri= "ring @Phelp ="
    1106                   +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
    1107        }
    1108        execute(@ri);
    1109        list @lvec;
    1110        list @lr=imap(@P,lres0);
    1111        ideal @lr1;
    1112 
    1113        if(size(@lr)==2)
    1114        {
    1115           @lr[2]=homog(@lr[2],@t);
    1116           @lr1=std(@lr[2]);
    1117           @lvec[2]=hilb(@lr1,1);
    1118        }
    1119        else
    1120        {
    1121           for(@n=1;@n<=size(@lr)/2;@n++)
     1173          for(@n=1;@n<=size(lres0)/2;@n++)
    11221174          {
    1123              if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
     1175             if(specialIdealsEqual(lres0[2*@n-1],lres0[2*@n])==1)
    11241176             {
    1125                 @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
    1126                 @lr1=std(@lr[2*@n-1]);
    1127                 @lvec[2*@n-1]=hilb(@lr1,1);
    1128                 @lvec[2*@n]=@lvec[2*@n-1];
     1177                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
     1178                lres0[2*@n]=lres0[2*@n-1];
     1179                attrib(lres0[2*@n],"isSB",1);
    11291180             }
    11301181             else
    11311182             {
    1132                 @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
    1133                 @lr1=std(@lr[2*@n-1]);
    1134                 @lvec[2*@n-1]=hilb(@lr1,1);
    1135                 @lr[2*@n]=homog(@lr[2*@n],@t);
    1136                 @lr1=std(@lr[2*@n]);
    1137                 @lvec[2*@n]=hilb(@lr1,1);
    1138 
    1139              }
    1140          }
    1141        }
    1142        @ri= "ring @Phelp1 ="
    1143                   +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
    1144        execute(@ri);
    1145        list @lr=imap(@Phelp,@lr);
    1146 
    1147        kill @Phelp;
    1148        if(size(@lr)==2)
    1149        {
    1150           @lr[2]=std(@lr[2],@lvec[2]);
    1151           @lr[2]=subst(@lr[2],@t,1);
    1152 
    1153        }
    1154        else
    1155        {
    1156           for(@n=1;@n<=size(@lr)/2;@n++)
    1157           {
    1158              if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
    1159              {
    1160                 @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
    1161                 @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
    1162                 @lr[2*@n]=@lr[2*@n-1];
    1163                 attrib(@lr[2*@n],"isSB",1);
    1164              }
    1165              else
    1166              {
    1167                 @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
    1168                 @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
    1169                 @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
    1170                 @lr[2*@n]=subst(@lr[2*@n],@t,1);
     1183                lres0[2*@n-1]=groebner(lres0[2*@n-1]);
     1184                lres0[2*@n]=groebner(lres0[2*@n]);
    11711185             }
    11721186          }
    11731187       }
    1174        kill @lvec;
    1175        setring @P;
    1176        lres0=imap(@Phelp1,@lr);
    1177        kill @Phelp1;
    1178        for(@n=1;@n<=size(lres0);@n++)
    1179        {
    1180           lres0[@n]=clearSB(lres0[@n]);
    1181           attrib(lres0[@n],"isSB",1);
    1182        }
    1183 
    1184        primary[2*@k-1]=lres0[1];
    1185        primary[2*@k]=lres0[2];
    1186        @s=size(primary)/2;
    1187        for(@n=1;@n<=size(lres0)/2-1;@n++)
    1188        {
    1189          primary[2*@s+2*@n-1]=lres0[2*@n+1];
    1190          primary[2*@s+2*@n]=lres0[2*@n+2];
    1191        }
    1192        @k--;
     1188       primary=primary+lres0;
     1189
     1190//=============================================================
     1191//       if(npars(@P)>0)
     1192//       {
     1193//          @ri= "ring @Phelp ="
     1194//                  +string(char(@P))+",
     1195//                   ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);";
     1196//       }
     1197//       else
     1198//       {
     1199//          @ri= "ring @Phelp ="
     1200//                  +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);";
     1201//       }
     1202//       execute(@ri);
     1203//       list @lvec;
     1204//       list @lr=imap(@P,lres0);
     1205//       ideal @lr1;
     1206//
     1207//       if(size(@lr)==2)
     1208//       {
     1209//          @lr[2]=homog(@lr[2],@t);
     1210//          @lr1=std(@lr[2]);
     1211//          @lvec[2]=hilb(@lr1,1);
     1212//       }
     1213//       else
     1214//       {
     1215//          for(@n=1;@n<=size(@lr)/2;@n++)
     1216//          {
     1217//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
     1218//             {
     1219//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
     1220//                @lr1=std(@lr[2*@n-1]);
     1221//                @lvec[2*@n-1]=hilb(@lr1,1);
     1222//                @lvec[2*@n]=@lvec[2*@n-1];
     1223//             }
     1224//             else
     1225//             {
     1226//                @lr[2*@n-1]=homog(@lr[2*@n-1],@t);
     1227//                @lr1=std(@lr[2*@n-1]);
     1228//                @lvec[2*@n-1]=hilb(@lr1,1);
     1229//                @lr[2*@n]=homog(@lr[2*@n],@t);
     1230//                @lr1=std(@lr[2*@n]);
     1231//                @lvec[2*@n]=hilb(@lr1,1);
     1232//
     1233//             }
     1234//         }
     1235//       }
     1236//       @ri= "ring @Phelp1 ="
     1237//                  +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);";
     1238//       execute(@ri);
     1239//       list @lr=imap(@Phelp,@lr);
     1240//
     1241//       kill @Phelp;
     1242//       if(size(@lr)==2)
     1243//      {
     1244//          @lr[2]=std(@lr[2],@lvec[2]);
     1245//          @lr[2]=subst(@lr[2],@t,1);
     1246//
     1247//       }
     1248//       else
     1249//       {
     1250//          for(@n=1;@n<=size(@lr)/2;@n++)
     1251//          {
     1252//             if(specialIdealsEqual(@lr[2*@n-1],@lr[2*@n])==1)
     1253//             {
     1254//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
     1255//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
     1256//                @lr[2*@n]=@lr[2*@n-1];
     1257//                attrib(@lr[2*@n],"isSB",1);
     1258//             }
     1259//             else
     1260//             {
     1261//                @lr[2*@n-1]=std(@lr[2*@n-1],@lvec[2*@n-1]);
     1262//                @lr[2*@n-1]=subst(@lr[2*@n-1],@t,1);
     1263//                @lr[2*@n]=std(@lr[2*@n],@lvec[2*@n]);
     1264//                @lr[2*@n]=subst(@lr[2*@n],@t,1);
     1265//             }
     1266//          }
     1267//       }
     1268//       kill @lvec;
     1269//       setring @P;
     1270//       lres0=imap(@Phelp1,@lr);
     1271//       kill @Phelp1;
     1272//       for(@n=1;@n<=size(lres0);@n++)
     1273//       {
     1274//          lres0[@n]=clearSB(lres0[@n]);
     1275//          attrib(lres0[@n],"isSB",1);
     1276//       }
     1277//
     1278//       primary[2*@k-1]=lres0[1];
     1279//       primary[2*@k]=lres0[2];
     1280//       @s=size(primary)/2;
     1281//       for(@n=1;@n<=size(lres0)/2-1;@n++)
     1282//       {
     1283//         primary[2*@s+2*@n-1]=lres0[2*@n+1];
     1284//         primary[2*@s+2*@n]=lres0[2*@n+2];
     1285//       }
     1286//       @k--;
     1287//=============================================================
    11931288     }
    11941289  }
    … …  
    12051300   pr;
    12061301}
     1302///////////////////////////////////////////////////////////////////////////////
     1303proc extF(list l,list #)
     1304{
     1305//zero_dimensional primary decomposition after finite field extension
     1306   def R=basering;
     1307   int p=char(R);
     1308
     1309   if((p==0)||(p>13)||(npars(R)>0)){return(l);}
     1310
     1311   int ex=3;
     1312   if(size(#)>0){ex=#[1];}
     1313
     1314   list peek,peek1;
     1315   while(size(l)>0)
     1316   {
     1317      if(size(l[2])==0)
     1318      {
     1319         peek[size(peek)+1]=l[1];
     1320      }
     1321      else
     1322      {
     1323         peek1[size(peek1)+1]=l[1];
     1324         peek1[size(peek1)+1]=l[2];
     1325      }
     1326      l=delete(l,1);
     1327      l=delete(l,1);
     1328   }
     1329   if(size(peek)==0){return(peek1);}
     1330
     1331   string gnir="ring RH=("+string(p)+"^"+string(ex)+",a),("+varstr(R)+"),lp;";
     1332   execute(gnir);
     1333   string mp="minpoly="+string(minpoly)+";";
     1334   gnir="ring RL=("+string(p)+",a),("+varstr(R)+"),lp;";
     1335   execute(gnir);
     1336   execute(mp);
     1337   list L=imap(R,peek);
     1338   list pr, keep;
     1339   int i;
     1340   for(i=1;i<=size(L);i++)
     1341   {
     1342      attrib(L[i],"isSB",1);
     1343      pr=zero_decomp(L[i],0,0);
     1344      keep=keep+pr;
     1345   }
     1346   for(i=1;i<=size(keep);i++)
     1347   {
     1348      keep[i]=simplify(keep[i],1);
     1349   }
     1350   mp="poly pp="+string(minpoly)+";";
     1351
     1352   string gnir1="ring RS="+string(p)+",("+varstr(R)+",a),lp;";
     1353   execute(gnir1);
     1354   execute(mp);
     1355   list L=imap(RL,keep);
     1356
     1357   for(i=1;i<=size(L);i++)
     1358   {
     1359      L[i]=eliminate(L[i]+ideal(pp),a);
     1360   }
     1361   i=0;
     1362   int j;
     1363   while(i<size(L)/2-1)
     1364   {
     1365      i++;
     1366      j=i;
     1367      while(j<size(L)/2)
     1368      {
     1369         j++;
     1370         if(idealsEqual(L[2*i-1],L[2*j-1]))
     1371         {
     1372            L=delete(L,2*j-1);
     1373            L=delete(L,2*j-1);
     1374            j--;
     1375         }
     1376      }
     1377   }
     1378   setring R;
     1379   list re=imap(RS,L);
     1380   re=re+peek1;
     1381
     1382   return(extF(re,ex+1));
     1383}
     1384
     1385///////////////////////////////////////////////////////////////////////////////
     1386proc zeroSp(ideal i)
     1387{
     1388//preparation for the separable closure
     1389//decomposition into ideals of special type
     1390//i.e. the minimal polynomials of every variable mod i are irreducible
     1391//returns a list of 2 lists: rr=pe,qe
     1392//the ideals in pe[l] are special, their special elements are in qe[l]
     1393//pe[l] is a dp-Groebnerbasis
     1394//the radical of the intersection of the pe[l] is equal to the radical of i
     1395
     1396   def R=basering;
     1397
     1398   //i has to be a reduced groebner basis
     1399   ideal F=finduni(i);
     1400
     1401   int j,k,l,ready;
     1402   list fa;
     1403   fa[1]=factorize(F[1],1);
     1404   poly te,ti;
     1405   ideal tj;
     1406   //avoid factorization of the same polynomial
     1407   for(j=2;j<=size(F);j++)
     1408   {
     1409      for(k=1;k<=j-1;k++)
     1410      {
     1411         ti=F[k];
     1412         te=subst(ti,var(k),var(j));
     1413         if(te==F[j])
     1414         {
     1415            tj=fa[k];
     1416            fa[j]=subst(tj,var(k),var(j));
     1417            ready=1;
     1418            break;
     1419         }
     1420      }
     1421      if(!ready)
     1422      {
     1423         fa[j]=factorize(F[j],1);
     1424      }
     1425      ready=0;
     1426   }
     1427   execute( "ring P=("+charstr(R)+"),("+varstr(R)+"),(C,dp);");
     1428   ideal i=imap(R,i);
     1429   if(npars(basering)==0)
     1430   {
     1431      ideal J=fglm(R,i);
     1432   }
     1433   else
     1434   {
     1435      ideal J=groebner(i);
     1436   }
     1437   list fa=imap(R,fa);
     1438   list qe=J;          //collects a dp-Groebnerbasis of the special ideals
     1439   list keep=ideal(0); //collects the special elements
     1440
     1441   list re,em,ke;
     1442   ideal K,L;
     1443
     1444   for(j=1;j<=nvars(basering);j++)
     1445   {
     1446      for(l=1;l<=size(qe);l++)
     1447      {
     1448         for(k=1;k<=size(fa[j]);k++)
     1449         {
     1450            L=std(qe[l],fa[j][k]);
     1451            K=keep[l],fa[j][k];
     1452            if(deg(L[1])>0)
     1453            {
     1454               re[size(re)+1]=L;
     1455               ke[size(ke)+1]=K;
     1456            }
     1457         }
     1458      }
     1459      qe=re;
     1460      re=em;
     1461      keep=ke;
     1462      ke=em;
     1463   }
     1464
     1465   setring R;
     1466   list qe=imap(P,keep);
     1467   list pe=imap(P,qe);
     1468   for(l=1;l<=size(qe);l++)
     1469   {
     1470      qe[l]=simplify(qe[l],2);
     1471    }
     1472    list rr=pe,qe;
     1473    return(rr);
     1474}
     1475///////////////////////////////////////////////////////////////////////////////
     1476
     1477proc zeroSepClos(ideal I,ideal F)
     1478{
     1479//computes the separable closure of the special ideal I
     1480//F is the set of special elements of I
     1481//returns the separable closure sc(I) of I and an intvec v
     1482//such that sc(I)=preimage(frobenius definde by v)
     1483//i.e. var(i)----->var(i)^(p^v[i])
     1484
     1485   if(homog(I)==1){return(maxideal(1));}
     1486
     1487   //assume F[i] irreducible in I and depending only on var(i)
     1488
     1489   def R=basering;
     1490   int n=nvars(R);
     1491   int p=char(R);
     1492   intvec v;
     1493   v[n]=0;
     1494   int i,k;
     1495   list l;
     1496
     1497   for(i=1;i<=n;i++)
     1498   {
     1499      l[i]=sep(F[i],i);
     1500      F[i]=l[i][1];
     1501      if(l[i][2]>k){k=l[i][2];}
     1502   }
     1503
     1504   if(k==0){return(list(I,v));}        //the separable case
     1505   ideal m;
     1506
     1507   for(i=1;i<=n;i++)
     1508   {
     1509      m[i]=var(i)^(p^l[i][2]);
     1510      v[i]=l[i][2];
     1511   }
     1512   map phi=R,m;
     1513   ideal J=preimage(R,phi,I);
     1514   return(list(J,v));
     1515}
     1516///////////////////////////////////////////////////////////////////////////////
     1517
     1518proc insepDecomp(ideal i)
     1519{
     1520//decomposes i into special ideals
     1521//computes the prime decomposition of the special ideals
     1522//and transforms it back to a decomposition of i
     1523
     1524   def R=basering;
     1525   list pr=zeroSp(i);
     1526   int l,k;
     1527   list re,wo,qr;
     1528   ideal m=maxideal(1);
     1529   ideal K;
     1530   map phi=R,m;
     1531   int p=char(R);
     1532   intvec op=option(get);
     1533
     1534   for(l=1;l<=size(pr[1]);l++)
     1535   {
     1536      wo=zeroSepClos(pr[1][l],pr[2][l]);
     1537      for(k=1;k<=nvars(basering);k++)
     1538      {
     1539         m[k]=var(k)^(p^wo[2][k]);
     1540      }
     1541      phi=R,m;
     1542      qr=decomp(wo[1],2);
     1543
     1544      option(redSB);
     1545      for(k=1;k<=size(qr)/2;k++)
     1546      {
     1547         K=qr[2*k];
     1548         K=phi(K);
     1549         K=groebner(K);
     1550         re[size(re)+1]=zeroRad(K);
     1551      }
     1552      option(noredSB);
     1553   }
     1554   option(set,op);
     1555   return(re);
     1556}
     1557
     1558
    12071559///////////////////////////////////////////////////////////////////////////////
    12081560
    … …  
    16131965}
    16141966
    1615 
    1616 
    1617 
     1967static proc primT(ideal i)
     1968{
     1969   //assumes that all generators of i are irreducible
     1970   //i is standard basis
     1971
     1972   attrib(i,"isSB",1);
     1973   int j=size(i);
     1974   int k;
     1975   while(j>0)
     1976   {
     1977     if(deg(i[j])>1){break;}
     1978     j--;
     1979   }
     1980   if(j==0){return(1);}
     1981   if(deg(i[j])==vdim(i)){return(1);}
     1982   return(0);
     1983}
    16181984
    16191985
    … …  
    16281994   list q=simplifyIdeal(i);
    16291995   list re=maxideal(1);
    1630    int j,a;
     1996   int j,a,k;
    16311997   intvec op=option(get);
    16321998   map phi=P,q[2];
    16331999
    1634    i=groebner(q[1]);
     2000   if(npars(P)==0){option(redSB);}
     2001
     2002   i=std(q[1]);
     2003   if(dim(i)==-1){re=ideal(1);return(re);}
     2004   if((dim(i)==0)&&(npars(P)==0))
     2005   {
     2006      int di=vdim(i);
     2007      execute ("ring gnir=("+charstr(P)+"),("+varstr(P)+"),lp;");
     2008      ideal J=interred(imap(P,i));
     2009      attrib(J,"isSB",1);
     2010      if(vdim(J)!=di)
     2011      {
     2012         J=fglm(P,i);
     2013      }
     2014      list pr=triangMH(J,2);
     2015      list qr,re;
     2016
     2017      for(k=1;k<=size(pr);k++)
     2018      {
     2019         if(primT(pr[k]))
     2020         {
     2021            re[size(re)+1]=pr[k];
     2022         }
     2023         else
     2024         {
     2025            attrib(pr[k],"isSB",1);
     2026            qr=decomp(pr[k],2);
     2027            for(j=1;j<=size(qr)/2;j++)
     2028            {
     2029               re[size(re)+1]=qr[2*j];
     2030            }
     2031         }
     2032      }
     2033      setring P;
     2034      re=imap(gnir,re);
     2035      option(set,op);
     2036      return(phi(re));
     2037   }
     2038
    16352039   if((size(#)==0)||(dim(i)==0))
    16362040   {
    16372041      re[1]=decomp(i,2);
    16382042      re=union(re);
     2043      option(set,op);
    16392044      return(phi(re));
    16402045   }
    1641 
    1642    option(redSB);
    16432046
    16442047   q=facstd(i);
    … …  
    19832386   ideal i = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5));
    19842387   equidimMax(i);
     2388}
     2389///////////////////////////////////////////////////////////////////////////////
     2390static proc islp()
     2391{
     2392   string s=ordstr(basering);
     2393   int n=find(s,"lp");
     2394   if(!n){return(0);}
     2395   int k=find(s,",");
     2396   string t=s[k+1..size(s)];
     2397   int l=find(t,",");
     2398   t=s[1..k-1];
     2399   int m=find(t,",");
     2400   if(l+m){return(0);}
     2401   return(1);
     2402}
     2403///////////////////////////////////////////////////////////////////////////////
     2404
     2405proc algeDeco(ideal i, int w)
     2406{
     2407//reduces primery decomposition over algebraic extensions to
     2408//the other cases
     2409   def R=basering;
     2410   int n=nvars(R);
     2411   string mp="poly p="+string(minpoly)+";";
     2412   string gnir="ring RH="+string(char(R))+",("+varstr(R)+","+string(par(1))
     2413                +"),dp;";
     2414   execute(gnir);
     2415   execute(mp);
     2416   ideal i=imap(R,i);
     2417         i=i,p;
     2418   list pr;
     2419   int j;
     2420
     2421   if(w==0)
     2422   {
     2423      pr=decomp(i);
     2424   }
     2425   if(w==1)
     2426   {
     2427      pr=prim_dec(i,1);
     2428      pr=reconvList(pr);
     2429   }
     2430   if(w==2)
     2431   {
     2432      pr=minAssPrimes(i);
     2433   }
     2434
     2435   gnir="ring RS="+string(char(R))+",("+varstr(RH)
     2436                +"),(dp("+string(n)+"),lp);";
     2437   execute(gnir);
     2438   list pr=imap(RH,pr);
     2439   ideal K;
     2440   for(j=1;j<=size(pr);j++)
     2441   {
     2442      K=groebner(pr[j]);
     2443      K=K[2..size(K)];
     2444      pr[j]=K;
     2445   }
     2446   setring R;
     2447   list pr=imap(RS,pr);
     2448      list re;
     2449   if(w==2)
     2450   {
     2451      re=pr;
     2452   }
     2453   else
     2454   {
     2455      re=convList(pr);
     2456   }
     2457   return(re);
    19852458}
    19862459
    … …  
    20042477  ideal peek=i;
    20052478  ideal ser,tras;
     2479  int isS=(attrib(i,"isSB")==1);
    20062480
    20072481  if(size(#)>0)
    … …  
    20362510    {
    20372511      ltras=i,i;
     2512      attrib(ltras[1],"isSB",1);
    20382513    }
    20392514    tras=ltras[1];
     2515    attrib(tras,"isSB",1);
    20402516    if(dim(tras)==0)
    20412517    {
    … …  
    20732549  //----------------------------------------------------------------
    20742550
     2551  int lp=islp();
     2552
    20752553  execute("ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);");
    20762554  op=option(get);
    … …  
    20812559  if(homo==1)
    20822560  {
    2083      if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)"))
     2561     if(!lp)
    20842562     {
    20852563       ideal @j=std(fetch(@P,i),@hilb,@w);
    … …  
    20932571  else
    20942572  {
    2095      ideal @j=groebner(fetch(@P,i));
     2573      if(lp&&isS)
     2574      {
     2575        ideal @j=fetch(@P,i);
     2576        attrib(@j,"isSB",1);
     2577      }
     2578      else
     2579      {
     2580        ideal @j=groebner(fetch(@P,i));
     2581      }
    20962582  }
    20972583  option(set,op);
    … …  
    22142700    op=option(get);
    22152701    option(redSB);
     2702
    22162703    list gprimary= zero_decomp(@j,ser,@wr);
    22172704    option(set,op);
    … …  
    28543341
    28553342///////////////////////////////////////////////////////////////////////////////
    2856 static proc zeroRad(ideal I,list #)
     3343 proc zeroRad(ideal I,list #)
    28573344"USAGE:  zeroRad(I) , I a zero-dimensional ideal
    28583345 RETURN: the radical of I
    … …  
    28893376      F[i]=powerCoeffs(F[i],k-l[i][2]);
    28903377   }
     3378
    28913379   string cR="ring @R="+string(p)+",("+parstr(R)+","+varstr(R)+"),dp;";
    28923380   execute(cR);
    28933381   ideal F=imap(R,F);
     3382
    28943383   string nR="ring @S="+string(p)+",(y(1..m),"+varstr(R)+","+parstr(R)+"),dp;";
    28953384   execute(nR);
    28963385
    28973386   ideal G=fetch(@R,F);    //G[i](t(1)^(p^-k),...,t(m)^(p^-k),x(i))=sep(F[i])
     3387
    28983388   ideal I=imap(R,I);
    28993389   ideal J=I+G;
     3390   poly el=1;
    29003391   k=p^k;
    29013392   for(i=1;i<=m;i++)
    29023393   {
    29033394      J=J,var(i)^k-var(m+n+i);
    2904    }
    2905    J=eliminate(J,y(1..m));
    2906 
     3395      el=el*y(i);
     3396   }
     3397
     3398   J=eliminate(J,el);
    29073399   setring R;
    29083400   ideal J=imap(@S,J);
    … …  
    29133405   ring R=(5,t),(x,y),dp;
    29143406   ideal I=x^5-t,y^5-t;
    2915    zeroradp(I);
     3407   zeroRad(I);
    29163408}
    29173409
    … …  
    32583750///////////////////////////////////////////////////////////////////////////////
    32593751
    3260 static proc Ann(module M)
     3752proc Ann(module M)
    32613753{
    32623754  M=prune(M);  //to obtain a small embedding
    … …  
    41374629
    41384630      U[size(U)]=")";           // we compute the extractor of I (w.r.t. U)
    4139       execute("ring RAU="+charstr(basering)+",("+A+U+",(dp("+string(a)+"),dp);");
     4631      execute("ring RAU=("+charstr(basering)+"),("+A+U+",(dp("+string(a)+"),dp);");
    41404632      ideal I=imap(R,SI);
    41414633      //I=std(I,hv);            // the standard basis in (R[U])[A]
    … …  
    44364928      );
    44374929   }
    4438    return(convList(decomp(i)));
     4930   if(minpoly!=0)
     4931   {
     4932      return(algeDeco(i,0));
     4933      ERROR(
     4934      "// Not implemented yet for algebraic extensions.Simulate the ring extension by adding the minpoly to the ideal"
     4935      );
     4936   }
     4937  return(convList(decomp(i)));
    44394938}
    44404939example
    … …  
    44634962   if c=3,  minAssGTZ and facstd are used.
    44644963@end format
    4465          Due to a bug in the factorization, the result may be not completely
    4466          decomposed in small characteristic.
    44674964EXAMPLE: example primdecSY; shows an example
    44684965"
    … …  
    44754972   }
    44764973   i=simplify(i,2);
    4477    if (i[1]==0)
    4478    {
    4479      list L=list(ideal(0),ideal(0));
     4974   if ((i[1]==0)||(i[1]==1))
     4975   {
     4976     list L=list(ideal(i[1]),ideal(i[1]));
    44804977     return(list(L));
    4481    } 
     4978   }
     4979   if(minpoly!=0)
     4980   {
     4981      return(algeDeco(i,1));
     4982   }
    44824983   if (size(#)==1)
    44834984   { return(prim_dec(i,#[1])); }
    … …  
    44995000          minAssGTZ(i,1); i ideal  does not use the factorizing Groebner
    45005001RETURN:  a list, the minimal associated prime ideals of i.
    4501 NOTE:    Designed for characteristic 0, works also in char k > 0 if it
    4502          terminates, may result in an infinite loop in small characteristic.
     5002NOTE:    Designed for characteristic 0, works also in char k > 0 based
     5003         on an algorithm of Yokoyama
    45035004EXAMPLE: example minAssGTZ; shows an example
    45045005"
    … …  
    45105011      );
    45115012   }
    4512    if(size(#)==0)
     5013   if(minpoly!=0)
     5014   {
     5015      return(algeDeco(i,2));
     5016      ERROR(
     5017      "// Not implemented for algebraic extensions. Simulate the ring extension by adding the minpoly to the ideal"
     5018      );
     5019   }
     5020  if(size(#)==0)
    45135021   {
    45145022      return(minAssPrimes(i,1));
    … …  
    46205128
    46215129   option(redSB);
    4622    list pr=facstd(i);
     5130   list pr=i;
     5131   if (!homog(i))
     5132   {
     5133     pr=facstd(i);
     5134   }
    46235135   option(set,op);
    46245136   int s=size(pr);
    … …  
    50115523time=timer; list pr =primdecSY(gls); timer-time;
    50125524time=timer; ideal ra =radical(gls); timer-time;size(pr);
     5525LIB"all.lib";
     5526
     5527ring R=0,(a,b,c,d,e,f),dp;
     5528ideal I=cyclic(6);
     5529minAssGTZ(I);
     5530
     5531
     5532ring S=(2,a,b),(x,y),lp;
     5533ideal I=x8-b,y4+a;
     5534minAssGTZ(I);
     5535
     5536ring S1=2,(x,y,a,b),lp;
     5537ideal I=x8-b,y4+a;
     5538minAssGTZ(I);
     5539
     5540
     5541ring S2=(2,z),(x,y),dp;
     5542minpoly=z2+z+1;
     5543ideal I=y3+y+1,x4+x+1;
     5544primdecGTZ(I);
     5545minAssGTZ(I);
     5546
     5547ring S3=2,(x,y,z),dp;
     5548ideal I=y3+y+1,x4+x+1,z2+z+1;
     5549primdecGTZ(I);
     5550minAssGTZ(I);
     5551
     5552
     5553ring R1=2,(x,y,z),lp;
     5554ideal I=y6+y5+y3+y2+1,x4+x+1,z2+z+1;
     5555primdecGTZ(I);
     5556minAssGTZ(I);
     5557
     5558
     5559ring R2=(2,z),(x,y),lp;
     5560minpoly=z3+z+1;
     5561ideal I=y2+y+(z2+z+1),x4+x+1;
     5562primdecGTZ(I);
     5563minAssGTZ(I);
     5564
    50135565*/
Note: See TracChangeset for help on using the changeset viewer.