Changeset c05763 in git


Ignore:
Timestamp:
Oct 4, 2019, 6:07:35 PM (5 years ago)
Author:
Sachin <sachinkm308@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'c5facdfddea2addfd91babd8b9019161dea4b695')
Children:
6546d7d800293348c495c339b2e608da49fd41ab
Parents:
f082e7584b5654aa2820b8af68ba222b71f4d903
git-author:
Sachin <sachinkm308@gmail.com>2019-10-04 18:07:35+02:00
git-committer:
Sachin <sachinkm308@gmail.com>2019-10-04 18:09:28+02:00
Message:
replacing execute with create_ring(8)
Location:
Singular/LIB
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/autgradalg.lib

    rf082e7 rc05763  
    12431243  }
    12441244
    1245   string Sstr = "ring S = (" + charstr(basering) + "),(T(1..n0), Y(1..n0*n)),dp;";
    1246   execute(Sstr); // creates the ring 0,(T(1..n0), Y(1..n0*n)),dp;
     1245  int ii;
     1246  list l3;
     1247  for (ii = 1; ii <= n0; ii++)
     1248  {
     1249   l3[ii] = "T("+string(ii)+")";
     1250  }
     1251  for (ii = 1; ii <= n0*n; ii++)
     1252  {
     1253   l3[n0+ii] = "Y("+string(ii)+")";
     1254  }
     1255  ring S = create_ring(ringlist(basering)[1], l3, "dp("+string(n0+n0*n)+")", "no_minpoly");
    12471256  setring S;
    12481257  setBaseMultigrading(QS);
    … …  
    17881797  }
    17891798
    1790   string Sstr = "ring S = (" + charstr(basering) + "),(T(1..n0), Y(1..n1),Z),(dp(n0+n1),dp(1));";
    1791   execute(Sstr);
     1799  int ii;
     1800  list l4;
     1801  for (ii = 1; ii <= n0; ii++)
     1802  {
     1803   l4[ii] = "T("+string(ii)+")";
     1804  }
     1805  for (ii = 1; ii <= n1; ii++)
     1806  {
     1807   l4[n0+ii] = "Y("+string(ii)+")";
     1808  }
     1809  l4[size(l4)+1] = "Z";
     1810  ring S = create_ring(ringlist(basering)[1], l4, "(dp("+string(n0+n1)+"),dp(1))", "no_minpoly");
    17921811  setring S;
    17931812  setBaseMultigrading(QS);
    … …  
    30353054  }
    30363055
    3037   string Sstr = "ring S = (" + charstr(RR) + "),(T(1..size(V))),dp;";
    3038   execute(Sstr); // creates the ring 0,(T(1..size(V))),dp;
     3056  list l5;
     3057  for (int ii = 1; ii <= size(V); ii++)
     3058  {
     3059   l5[ii] = "T("+string(ii)+")";
     3060  }
     3061  ring S = create_ring(ringlist(RR)[1], l5, "dp", "no_minpoly");
    30393062  setring S;
    30403063  setring RR;

  • Singular/LIB/classify.lib

    rf082e7 rc05763  
    160160  @ringdisplay = "setring RingB;";
    161161  if(defined(RingB)!=0) { kill RingB; }
    162   execute ("ring RingB="+charstr(basering)+",("+A_Z("x", n)+"),(c,ds);");
     162  ring RingB = create_ring(ringlist(basering)[1], "("+A_Z("x", n)+")", "(c,ds)", "no_minpoly");
    163163  export RingB;
    164164  setring ring_top;

  • Singular/LIB/ehv.lib

    rf082e7 rc05763  
    821821  def base = basering;
    822822  j = j[1,size(j)-1];
    823   j = "ring @r = (" + charstr(basering) + "),(" + j + "),(" + ordstr(basering) + ");";
    824   execute(j);
     823  j;
     824  ring @r = create_ring(ringlist(basering)[1], (" + j + "), "(" + ordstr(basering) + ")", "no_minpoly");
    825825  ideal J = imap(base,J);
    826826
    … …  
    10441044                l2[i] = "x("+string(ii)+")";
    10451045              }
    1046               l2[n+1] = "t";
     1046              l2[size(l2)+1] = "t";
    10471047              ring newR = create_ring(ringlist(base)[1], l2, "dp", "no_minpoly");
    10481048              ideal homoK = fetch(homoR,homoJ);

  • Singular/LIB/equising.lib

    rf082e7 rc05763  
    1 ///////////////////////////////////////////////////////////////////////
    2 version="version equising.lib 4.1.2.0 Feb_2019 "; // $Id$
    3 category="Singularities";
    4 info="
    5 LIBRARY:  equising.lib  Equisingularity Stratum of a Family of Plane Curves
    6 AUTHOR:   Christoph Lossen, lossen@mathematik.uni-kl.de
    7           Andrea Mindnich, mindnich@mathematik.uni-kl.de
    8 
    9 PROCEDURES:
    10  tau_es(f);             codim of mu-const stratum in semi-universal def. base
    11  esIdeal(f);            (Wahl's) equisingularity ideal of f
    12  esStratum(F[,m,L]);    equisingularity stratum of a family F
    13  isEquising(F[,m,L]);   tests if a given deformation is equisingular
    14 
    15  control_Matrix(M);     computes list of blowing-up data
    16 ";
    17 
    18 LIB "hnoether.lib";
    19 LIB "poly.lib";
    20 LIB "elim.lib";
    21 LIB "deform.lib";
    22 LIB "sing.lib";
    23 
    24 ////////////////////////////////////////////////////////////////////////////////
    25 //
    26 //  The following (static) procedures are used by esComputation
    27 //
    28 ////////////////////////////////////////////////////////////////////////////////
    29 // COMPUTES  a weight vector. x and y get weight 1 and all other
    30 //           variables get weight 0.
    31 static proc xyVector()
    32 {
    33   intvec iv ;
    34   iv[nvars(basering)]=0 ;
    35   iv[rvar(x)] =1;
    36   iv[rvar(y)] =1;
    37   return (iv);
    38 }
    39 
    40 ////////////////////////////////////////////////////////////////////////////////
    41 // exchanges the variables x and y in the polynomial f
    42 static proc swapXY(poly f)
    43 {
    44   def r_base = basering;
    45   ideal MI = maxideal(1);
    46   MI[rvar(x)]=y;
    47   MI[rvar(y)]=x;
    48   map phi = r_base, MI;
    49   f=phi(f);
    50   return (f);
    51 }
    52 
    53 ////////////////////////////////////////////////////////////////////////////////
    54 // computes m-jet w.r.t. the variables x,y (other variables weighted 0
    55 static proc m_Jet(poly F,int m);
    56 {
    57   intvec w=xyVector();
    58   poly Fd=jet(F,m,w);
    59   return(Fd);
    60 }
    61 
    62 
    63 ////////////////////////////////////////////////////////////////////////////////
    64 // computes the 4 control matrices (input is multsequence(L))
    65 proc control_Matrix(list M);
    66 "USAGE:   control_Matrix(L); L list
    67 ASSUME:  L is the output of multsequence(hnexpansion(f)).
    68 RETURN:  list M of 4 intmat's:
    69 @format
    70   M[1] contains the multiplicities at the respective infinitely near points
    71        p[i,j] (i=step of blowup+1, j=branch) -- if branches j=k,...,k+m pass
    72        through the same p[i,j] then the multiplicity is stored in M[1][k,j],
    73        while M[1][k+1]=...=M[1][k+m]=0.
    74   M[2] contains the number of branches meeting at p[i,j] (again, the information
    75        is stored according to the above rule)
    76   M[3] contains the information about the splitting of M[1][i,j] with respect to
    77        different tangents of branches at p[i,j] (information is stored only for
    78        minimal j>=k corresponding to a new tangent direction).
    79        The entries are the sum of multiplicities of all branches with the
    80        respective tangent.
    81   M[4] contains the maximal sum of higher multiplicities for a branch passing
    82        through p[i,j] ( = degree Bound for blowing up)
    83 @end format
    84 NOTE:    the branches are ordered in such a way that only consecutive branches
    85          can meet at an infinitely near point. @*
    86          the final rows of the matrices M[1],...,M[3] is (1,1,1,...,1), and
    87          correspond to infinitely near points such that the strict transforms
    88          of the branches are smooth and intersect the exceptional divisor
    89          transversally.
    90 SEE ALSO: multsequence
    91 EXAMPLE: example control_Matrix; shows an example
    92 "
    93 {
    94   int i,j,k,dummy;
    95 
    96   dummy=0;
    97   for (j=1;j<=ncols(M[2]);j++)
    98   {
    99     dummy=dummy+M[1][nrows(M[1])-1,j]-M[1][nrows(M[1]),j];
    100   }
    101   intmat S[nrows(M[1])+dummy][ncols(M[1])];
    102   intmat T[nrows(M[1])+dummy][ncols(M[1])];
    103   intmat U[nrows(M[1])+dummy][ncols(M[1])];
    104   intmat maxDeg[nrows(M[1])+dummy][ncols(M[1])];
    105 
    106   for (i=1;i<=nrows(M[2]);i++)
    107   {
    108     dummy=1;
    109     for (j=1;j<=ncols(M[2]);j++)
    110     {
    111       for (k=dummy;k<dummy+M[2][i,j];k++)
    112       {
    113         T[i,dummy]=T[i,dummy]+1;
    114         S[i,dummy]=S[i,dummy]+M[1][i,k];
    115         if (i>1)
    116         {
    117           U[i-1,dummy]=U[i-1,dummy]+M[1][i-1,k];
    118         }
    119       }
    120       dummy=k;
    121     }
    122   }
    123 
    124   // adding an extra row (in some cases needed to control ES-Stratum
    125   // computation)
    126   for (i=nrows(M[1]);i<=nrows(S);i++)
    127   {
    128     for (j=1;j<=ncols(M[2]);j++)
    129     {
    130       S[i,j]=1;
    131       T[i,j]=1;
    132       U[i,j]=1;
    133     }
    134   }
    135 
    136   // Computing the degree Bounds to be stored in M[4]:
    137   for (i=1;i<=nrows(S);i++)
    138   {
    139     dummy=1;
    140     for (j=1;j<=ncols(S);j++)
    141     {
    142       for (k=dummy;k<dummy+T[i,j];k++)
    143       {
    144         maxDeg[i,k]=S[i,dummy];  // multiplicity at i-th blowup
    145       }
    146       dummy=k;
    147     }
    148   }
    149   // adding up multiplicities
    150   for (i=nrows(S);i>=2;i--)
    151   {
    152     for (j=1;j<=ncols(S);j++)
    153     {
    154       maxDeg[i-1,j]=maxDeg[i-1,j]+maxDeg[i,j];
    155     }
    156   }
    157 
    158   list L=S,T,U,maxDeg;
    159   return(L);
    160 }
    161 
    162 
    163 ////////////////////////////////////////////////////////////////////////////////
    164 //  matrix of higher tangent directions:
    165 //  returns list: 1) tangent directions
    166 //                2) swapping information (x <--> y)
    167 static proc inf_Tangents(list L,int s); // L aus hnexpansion,
    168 {
    169   int nv=nvars(basering);
    170   matrix M;
    171   matrix B[s][size(L)];
    172   intvec V;
    173   intmat Mult=multsequence(L)[1];
    174 
    175   int i,j,k,counter,e;
    176   for (k=1;k<=size(L);k++)
    177   {
    178     V[k]=L[k][3];  // switch: 0 --> tangent 2nd parameter
    179                    //         1 --> tangent 1st parameter
    180     e=0;
    181     M=L[k][1];
    182     counter=1;
    183     B[counter,k]=M[1,1];
    184 
    185     for (i=1;i<=nrows(M);i++)
    186     {
    187       for (j=2;j<=ncols(M);j++)
    188       {
    189         counter=counter+1;
    190         if (M[i,j]==var(nv-1))
    191         {
    192           if (i<>nrows(M))
    193           {
    194             B[counter,k]=M[i,j];
    195             j=ncols(M)+1; // goto new row of HNmatrix...
    196             if (counter<>s)
    197             {
    198               if (counter+1<=nrows(Mult))
    199               {
    200                 e=Mult[counter-1,k]-Mult[counter,k]-Mult[counter+1,k];
    201               }
    202               else
    203               {
    204                 e=Mult[counter-1,k]-Mult[counter,k]-1;
    205               }
    206             }
    207           }
    208           else
    209           {
    210             B[counter,k]=0;
    211             j=ncols(M)+1; // goto new row of HNmatrix...
    212           }
    213         }
    214         else
    215         {
    216           if (e<=0)
    217           {
    218             B[counter,k]=M[i,j];
    219           }
    220           else  // point is still proximate to an earlier point
    221           {
    222             B[counter,k]=y; // marking proximity (without swap....)
    223             if (counter+1<=nrows(Mult))
    224             {
    225               e=e-Mult[counter+1,k];
    226             }
    227             else
    228             {
    229               e=e-1;
    230             }
    231           }
    232         }
    233 
    234         if (counter==s) // given number of points determined
    235         {
    236             j=ncols(M)+1;
    237             i=nrows(M)+1;
    238             // leave procedure
    239         }
    240       }
    241     }
    242   }
    243   L=B,V;
    244   return(L);
    245 }
    246 
    247 ////////////////////////////////////////////////////////////////////////////////
    248 // compute "good" upper bound for needed number of help variables
    249 //
    250 static proc Determine_no_b(intmat U,matrix B)
    251 // U is assumed to be 3rd output of control_Matrix
    252 // B is assumed to be 1st output of inf_Tangents
    253 {
    254   int nv=nvars(basering);
    255   int i,j,counter;
    256   for (j=1;j<=ncols(U);j++)
    257   {
    258     for (i=1;i<=nrows(U);i++)
    259     {
    260       if (U[i,j]>1)
    261       {
    262         if (B[i,j]<>var(nv-1) and B[i,j]<>var(nv))
    263         {
    264           counter=counter+1;
    265         }
    266       }
    267 
    268     }
    269   }
    270   counter=counter+ncols(U);
    271   return(counter);
    272 }
    273 
    274 ////////////////////////////////////////////////////////////////////////////////
    275 // compute number of infinitely near free points corresponding to non-zero
    276 // entries in control_Matrix[1] (except first row)
    277 //
    278 static proc no_freePoints(intmat Mult,matrix B)
    279 // Mult is assumed to be 1st output of control_Matrix
    280 // U is assumed to be 3rd output of control_Matrix
    281 // B is assumed to be 1st output of inf_Tangents
    282 {
    283   int i,j,k,counter;
    284   for (j=1;j<=ncols(Mult);j++)
    285   {
    286     for (i=2;i<=nrows(Mult);i++)
    287     {
    288       if (Mult[i,j]>=1)
    289       {
    290         if (B[i-1,j]<>x and B[i-1,j]<>y)
    291         {
    292           counter=counter+1;
    293         }
    294       }
    295     }
    296   }
    297   return(counter);
    298 }
    299 
    300 
    301 ///////////////////////////////////////////////////////////////////////////////
    302 // COMPUTES string(minpoly) and substitutes the parameter by newParName
    303 static proc makeMinPolyString (string newParName)
    304 {
    305   int i;
    306   string parName = parstr(basering);
    307   int parNameSize = size(parName);
    308 
    309   string oldMinPolyStr = string (minpoly);
    310   int minPolySize = size(oldMinPolyStr);
    311 
    312   string newMinPolyStr = "";
    313 
    314   for (i=1;i <= minPolySize; i++)
    315   {
    316     if (oldMinPolyStr[i,parNameSize] == parName)
    317     {
    318       newMinPolyStr = newMinPolyStr + newParName;
    319       i = i + parNameSize-1;
    320     }
    321     else
    322     {
    323       newMinPolyStr = newMinPolyStr + oldMinPolyStr[i];
    324     }
    325   }
    326 
    327   return(newMinPolyStr);
    328 }
    329 
    330 
    331 ///////////////////////////////////////////////////////////////////////////////
    332 //
    333 // DEFINES: A new basering, "myRing",
    334 //          with new names for the parameters and variables.
    335 //          The new names for the parameters are a(1..k),
    336 //          and t(1..s),x,y for the variables
    337 //          The ring ordering is ordStr.
    338 // NOTE:    This proc uses 'execute'.
    339 static proc createMyRing_new(poly p_F, string ordStr,
    340                                 string minPolyStr, int no_b)
    341 {
    342   def r_old = basering;
    343 
    344   int chara = char(basering);
    345   string charaStr;
    346   int i;
    347   string helpStr;
    348   int nDefParams = nvars(r_old)-2;
    349 
    350   ideal qIdeal = ideal(basering);
    351 
    352   if ((npars(basering)==0) and (minPolyStr==""))
    353   {
    354     helpStr = "ring myRing1 ="
    355               + string(chara)+ ", (t(1..nDefParams), x, y),("+ ordStr +");";
    356     execute(helpStr);
    357   }
    358   else
    359   {
    360     charaStr = charstr(basering);
    361     if ((charaStr == string(chara) + "," + parstr(basering)) or (minPolyStr<>""))
    362     {
    363       if (minPolyStr<>"")
    364       {
    365         helpStr = "ring myRing1 =
    366                  (" + string(chara) + ",a),
    367                  (t(1..nDefParams), x, y),(" + ordStr + ");";
    368         execute(helpStr);
    369 
    370         execute (minPolyStr);
    371       }
    372       else // no minpoly given
    373       {
    374         helpStr = "ring myRing1 =
    375                   (" + string(chara) + ",a(1..npars(basering)) ),
    376                   (t(1..nDefParams), x, y),(" + ordStr + ");";
    377         execute(helpStr);
    378       }
    379     }
    380     else
    381     {
    382       // ground field is of type (p^k,a)....
    383       i = find (charaStr,",");
    384       helpStr = "ring myRing1 = (" + charaStr[1,i] + "a),
    385               (t(1..nDefParams), x, y),(" + ordStr + ");";
    386       execute (helpStr);
    387     }
    388   }
    389 
    390   ideal mIdeal = maxideal(1);
    391   ideal qIdeal = fetch(r_old, qIdeal);
    392   poly p_F = fetch(r_old, p_F);
    393   export p_F,mIdeal;
    394 
    395   // Extension by no_b auxiliary variables
    396   if (no_b>0)
    397   {
    398     if (npars(basering) == 0)
    399     {
    400       ordStr = "(dp("+string(no_b)+"),"+ordStr+")";
    401       helpStr = "ring myRing ="
    402                 + string(chara)+ ", (b(1..no_b), t(1..nDefParams), x, y),"
    403                 + ordStr +";";
    404       execute(helpStr);
    405     }
    406     else
    407     {
    408       charaStr = charstr(basering);
    409       if (charaStr == string(chara) + "," + parstr(basering))
    410       {
    411         if (minpoly !=0)
    412         {
    413           ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")";
    414           minPolyStr = makeMinPolyString("a");
    415           helpStr = "ring myRing =
    416                    (" + string(chara) + ",a),
    417                    (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";";
    418           execute(helpStr);
    419 
    420           helpStr = "minpoly =" + minPolyStr + ";";
    421           execute (helpStr);
    422         }
    423         else // no minpoly given
    424         {
    425           ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")";
    426           helpStr = "ring myRing =
    427                     (" + string(chara) + ",a(1..npars(basering)) ),
    428                     (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";";
    429           execute(helpStr);
    430         }
    431       }
    432       else
    433       {
    434         i = find (charaStr,",");
    435         ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")";
    436         helpStr = "ring myRing =
    437                 (" + charaStr[1,i] + "a),
    438                 (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";";
    439         execute (helpStr);
    440       }
    441     }
    442     ideal qIdeal = imap(myRing1, qIdeal);
    443 
    444     if(size(qIdeal) != 0)
    445     {
    446       def r_base = basering;
    447       setring r_base;
    448       kill myRing;
    449       qring myRing = std(qIdeal);
    450     }
    451 
    452     poly p_F = imap(myRing1, p_F);
    453     ideal mIdeal = imap(myRing1, mIdeal);
    454     export p_F,mIdeal;
    455     kill myRing1;
    456   }
    457   else
    458   {
    459     if(size(qIdeal) != 0)
    460     {
    461       def r_base = basering;
    462       setring r_base;
    463       kill myRing1;
    464       qring myRing = std(qIdeal);
    465       poly p_F = imap(myRing1, p_F);
    466       ideal mIdeal = imap(myRing1, mIdeal);
    467       export p_F,mIdeal;
    468     }
    469     else
    470     {
    471       def myRing=myRing1;
    472     }
    473     kill myRing1;
    474   }
    475 
    476   setring r_old;
    477   return(myRing);
    478 }
    479 
    480 ////////////////////////////////////////////////////////////////////////////////
    481 // returns list of coef, leadmonomial
    482 //
    483 static proc determine_coef (poly Fm)
    484 {
    485   def r_base = basering; // is assumed to be the result of createMyRing
    486 
    487   int chara = char(basering);
    488   string charaStr;
    489   int i;
    490   string minPolyStr = "";
    491   string helpStr = "";
    492 
    493   if (npars(basering) == 0)
    494   {
    495     helpStr = "ring myRing1 ="
    496               + string(chara)+ ", (y,x),ds;";
    497     execute(helpStr);
    498   }
    499   else
    500   {
    501     charaStr = charstr(basering);
    502     if (charaStr == string(chara) + "," + parstr(basering))
    503     {
    504       if (minpoly !=0)
    505       {
    506         minPolyStr = makeMinPolyString("a");
    507         helpStr = "ring myRing1 = (" + string(chara) + ",a), (y,x),ds;";
    508         execute(helpStr);
    509 
    510         helpStr = "minpoly =" + minPolyStr + ";";
    511         execute (helpStr);
    512       }
    513       else // no minpoly given
    514       {
    515         helpStr = "ring myRing1 =
    516                   (" + string(chara) + ",a(1..npars(basering)) ), (y,x),ds;";
    517         execute(helpStr);
    518       }
    519     }
    520     else
    521     {
    522       i = find (charaStr,",");
    523 
    524       helpStr = " ring myRing1 = (" + charaStr[1,i] + "a), (y,x),ds;";
    525       execute (helpStr);
    526     }
    527   }
    528   poly f=imap(r_base,Fm);
    529   poly g=leadmonom(f);
    530   setring r_base;
    531   poly g=imap(myRing1,g);
    532   kill myRing1;
    533   def M=coef(Fm,xy);
    534 
    535   for (i=1; i<=ncols(M); i++)
    536   {
    537     if (M[1,i]==g)
    538     {
    539       poly h=M[2,i];  // determine coefficient of leading monomial (in K[t])
    540       i=ncols(M)+1;
    541     }
    542   }
    543   return(list(h,g));
    544 }
    545 
    546 ///////////////////////////////////////////////////////////////////////////////
    547 // RETURNS: 1, if p_f = 0 or char(basering) divides the order of p_f
    548 //             or p_f is not squarefree.
    549 //          0, otherwise
    550 static proc checkPoly (poly p_f)
    551 {
    552   int i_print = printlevel - voice + 3;
    553   int i_ord;
    554 
    555   if (p_f == 0)
    556   {
    557     print("Input is a 'deformation'  of the zero polynomial!");
    558     return(1);
    559   }
    560 
    561   i_ord = mindeg1(p_f);
    562 
    563   if (number(i_ord) == 0)
    564   {
    565     print("Characteristic of coefficient field "
    566             +"divides order of zero-fiber !");
    567     return(1);
    568   }
    569 
    570   if (squarefree(p_f) != p_f)
    571   {
    572     print("Original polynomial (= zero-fiber) is not reduced!");
    573     return(1);
    574   }
    575 
    576   return(0);
    577 }
    578 
    579 ////////////////////////////////////////////////////////////////////////////////
    580 static proc make_ring_small(ideal J)
    581 // returns varstr for new ring, the map and the number of vars
    582 {
    583   attrib(J,"isSB",1);
    584   int counter=0;
    585   ideal newmap;
    586   string newvar="";
    587   for (int i=1; i<=nvars(basering); i++)
    588   {
    589     if (reduce(var(i),J)<>0)
    590     {
    591       newmap[i]=var(i);
    592 
    593       if (newvar=="")
    594       {
    595         newvar=newvar+string(var(i));
    596         counter=counter+1;
    597       }
    598       else
    599       {
    600         newvar=newvar+","+string(var(i));
    601         counter=counter+1;
    602       }
    603     }
    604     else
    605     {
    606       newmap[i]=0;
    607     }
    608   }
    609   list L=newvar,newmap,counter;
    610   attrib(J,"isSB",0);
    611   return(L);
    612 }
    613 
    614 ///////////////////////////////////////////////////////////////////////////////
    615 //  The following procedure is called by esStratum (typ=0), resp. by
    616 //  isEquising (typ=1)
    617 ///////////////////////////////////////////////////////////////////////////////
    618 
    619 static proc esComputation (int typ, poly p_F, list #)
    620 {
    621   intvec ov=option(get);  // store options set at beginning
    622   option(redSB);
    623   // Initialize variables
    624   int branch=1;
    625   int blowup=1;
    626   int auxVar=1;
    627   int nVars;
    628 
    629   intvec upper_bound, upper_bound_old, fertig, soll;
    630   list blowup_string;
    631   int i_print= printlevel-voice+2;
    632 
    633   int no_b, number_of_branches, swapped;
    634   int i,j,k,m, counter, dummy;
    635   string helpStr = "";
    636   string ordStr = "";
    637   string MinPolyStr = "";
    638 
    639   if (nvars(basering)<=2)
    640   {
    641     print("family is trivial (no deformation parameters)!");
    642     if (typ==1) //isEquising
    643     {
    644       option(set,ov);
    645       return(1);
    646     }
    647     else
    648     {
    649       option(set,ov);
    650       return(list(ideal(0),0));
    651     }
    652   }
    653 
    654   if (size(#)>0)
    655   {
    656     if (typeof(#[1])=="int")
    657     {
    658       def artin_bd=#[1];  // compute modulo maxideal(artin_bd)
    659       if (artin_bd <= 1)
    660       {
    661         print("Do you really want to compute over Basering/maxideal("
    662               +string(artin_bd)+") ?");
    663         print("No computation performed !");
    664         if (typ==1) //isEquising
    665         {
    666           option(set,ov);
    667           return(1);
    668         }
    669         else
    670         {
    671           option(set,ov);
    672           return(list(ideal(0),int(1)));
    673         }
    674       }
    675       if (size(#)>1)
    676       {
    677         if (typeof(#[2])=="list")
    678         {
    679           def @L=#[2];  // is assumed to be the Hamburger-Noether matrix
    680         }
    681       }
    682     }
    683     else
    684     {
    685       if (typeof(#)=="list")
    686       {
    687         def @L=#;  // is assumed to be the Hamburger-Noether matrix
    688       }
    689     }
    690   }
    691   int ring_is_changed;
    692   def old_ring=basering;
    693   if(defined(@L)<=0)
    694   {
    695     // define a new ring without deformation-parameters and change to it:
    696     string str;
    697     string minpolyStr = string(minpoly);
    698     str = " ring HNERing = (" + charstr(basering) + "), (x,y), ls;";
    699     execute (str);
    700     if (minpolyStr!="0")
    701     {str = "minpoly ="+ minpolyStr+";";execute(str);}
    702     ring_is_changed=1;
    703     // Basering changed to HNERing (variables x,y, with ls ordering)
    704 
    705     k=nvars(old_ring);
    706     matrix Map_Phi[1][k];
    707     Map_Phi[1,k-1]=x;
    708     Map_Phi[1,k]=y;
    709     map phi=old_ring,Map_Phi;
    710     poly f=phi(p_F);
    711 
    712     // Heuristics: if x,y are transversal parameters then computation of HNE
    713     // can be much faster when exchanging variables...!
    714     if (2*size(coeffs(f,x))<size(coeffs(f,y)))
    715     {
    716       swapped=1;
    717       f=swapXY(f);
    718     }
    719 
    720     int error=checkPoly(f);
    721     if (error)
    722     {
    723       setring old_ring;
    724       if (typ==1) //isEquising
    725       {
    726         print("Return value (=0) has no meaning!");
    727         option(set,ov);
    728         return(0);
    729       }
    730       else
    731       {
    732         option(set,ov);
    733         return(list( ideal(0),error));
    734       }
    735     }
    736 
    737     dbprint(i_print,"// ");
    738     dbprint(i_print,"// Compute HN expansion");
    739     dbprint(i_print,"// ---------------------");
    740     i=printlevel;
    741     printlevel=printlevel-5;
    742     list LLL=hnexpansion(f);
    743 
    744     if (size(LLL)==0) { // empty list returned by hnexpansion
    745       setring old_ring;
    746       print(i_print,"Unable to compute HN expansion !");
    747       if (typ==1) //isEquising
    748       {
    749         print("Return value (=0) has no meaning!");
    750         option(set,ov);
    751         return(0);
    752       }
    753       else
    754       {
    755         option(set,ov);
    756         return(list(ideal(0),int(1)));
    757       }
    758       option(set,ov);
    759       return(0);
    760     }
    761     else
    762     {
    763       if (typeof(LLL[1])=="ring") {
    764         def HNering = LLL[1];
    765         setring HNering;
    766         def @L=stripHNE(hne);
    767       }
    768       else {
    769         def @L=stripHNE(LLL);
    770       }
    771     }
    772     printlevel=i;
    773     dbprint(i_print,"// finished");
    774     dbprint(i_print,"// ");
    775   }
    776   def HNEring=basering;
    777   list M=multsequence(@L);
    778   M=control_Matrix(M);     // this returns the 4 control matrices
    779   def maxDeg=M[4];
    780 
    781   list L1=inf_Tangents(@L,nrows(M[1]));
    782   matrix B=L1[1];
    783   intvec V=L1[2];
    784   kill L1;
    785 
    786   // if we have computed the HNE for f after swapping x and y, we have
    787   // to reinterprete the (swap) matrix V:
    788   if (swapped==1)
    789   {
    790     for (i=1;i<=size(V);i++) { V[i]=V[i]-1; } // turns 0 into -1, 1 into 0
    791   }
    792 
    793   // Determine maximal number of needed auxiliary parameters (free tangents):
    794   no_b=Determine_no_b(M[3],B);
    795 
    796   // test whether HNexpansion needed field extension....
    797   string minPolyStr = "";
    798   if (minpoly !=0)
    799   {
    800     minPolyStr = makeMinPolyString("a");
    801     minPolyStr = "minpoly =" + minPolyStr + ";";
    802   }
    803 
    804   setring old_ring;
    805 
    806   def myRing=createMyRing_new(p_F,"dp",minPolyStr,no_b);
    807   setring myRing;  // comes with mIdeal
    808   map hole=HNEring,mIdeal;
    809   // basering has changed to myRing, in particular, the "old"
    810   // variable names, e.g., A,B,C,z,y are replaced by t(1),t(2),t(3),x,y
    811 
    812   ideal bNodes;
    813 
    814   // Initialize some variables:
    815   map phi;
    816   poly G, F_save;
    817   poly b_dummy;
    818   ideal J,Jnew,final_Map;
    819   number_of_branches=ncols(M[1]);
    820   for (i=1;i<=number_of_branches;i++)
    821   {
    822     poly F(i);
    823     ideal bl_Map(i);
    824   }
    825   upper_bound[number_of_branches]=0;
    826   upper_bound[1]=number_of_branches;
    827   upper_bound_old=upper_bound;
    828   fertig[number_of_branches]=0;
    829   for (i=1;i<=number_of_branches;i++){ soll[i]=1; }
    830 
    831   // Hole:  B = matrix of blowup points
    832   if (ring_is_changed==0) { matrix B=hole(B); }
    833   else                    { matrix B=imap(HNEring,B); }
    834   m=M[1][blowup,branch];    // multiplicity at 0
    835 
    836   // now, we start by checking equimultiplicity along trivial section
    837   poly Fm=m_Jet(p_F,m-1);
    838 
    839   matrix coef_Mat = coef(Fm,xy);
    840   Jnew=coef_Mat[2,1..ncols(coef_Mat)];
    841   J=J,Jnew;
    842 
    843   if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
    844                          // deformation parameters can be cutted off
    845   {
    846     J=jet(J,artin_bd-1);
    847   }
    848 
    849   J=interred(J);
    850   if (defined(artin_bd)) { J=jet(J,artin_bd-1); }
    851 
    852 // J=std(J);
    853 
    854   if (typ==1) // isEquising
    855   {
    856     if(ideal(nselect(J,1..no_b))<>0)
    857     {
    858       setring old_ring;
    859       option(set,ov);
    860       return(0);
    861     }
    862   }
    863 
    864   F(1)=p_F;
    865 
    866   // and reduce the remaining terms in F(1):
    867   bl_Map(1)=maxideal(1);
    868 
    869   attrib(J,"isSB",1);
    870   bl_Map(1)=reduce(bl_Map(1),J);
    871   attrib(J,"isSB",0);
    872 
    873   phi=myRing,bl_Map(1);
    874   F(1)=phi(F(1));
    875 
    876   // simplify F(1)
    877   attrib(J,"isSB",1);
    878   F(1)=reduce(F(1),J);
    879   attrib(J,"isSB",0);
    880 
    881   // now we compute the m-jet:
    882   Fm=m_Jet(F(1),m);
    883 
    884   G=1;
    885   counter=branch;
    886   k=upper_bound[branch];
    887 
    888   F_save=F(1);  // is truncated differently in the following loop
    889 
    890   while(counter<=k)
    891   {
    892     F(counter)=m_Jet(F_save,maxDeg[blowup,counter]);
    893     if (V[counter]==0) // 2nd ring variable is tangent to this branch
    894     {
    895       G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]);
    896     }
    897     else // 1st ring variable is tangent to this branch
    898     {
    899       G=G*(x-(b(auxVar)+B[blowup,counter])*y)^(M[3][blowup,counter]);
    900       F(counter)=swapXY(F(counter));
    901     }
    902     bl_Map(counter)=maxideal(1);
    903     bl_Map(counter)[nvars(basering)]=xy+(b(auxVar)+B[blowup,counter])*x;
    904 
    905     bNodes[counter]=b(auxVar);
    906 
    907     auxVar=auxVar+1;
    908     upper_bound[counter]=counter+M[2][blowup+1,counter]-1;
    909     counter=counter+M[2][blowup+1,counter];
    910 
    911   }
    912 
    913   list LeadDataFm=determine_coef(Fm);
    914   def LeadDataG=coef(G,xy);
    915 
    916   for (i=1; i<=ncols(LeadDataG); i++)
    917   {
    918     if (LeadDataG[1,i]==LeadDataFm[2])
    919     {
    920       poly LeadG = LeadDataG[2,i];  // determine the coefficient of G
    921       i=ncols(LeadDataG)+1;
    922     }
    923   }
    924 
    925   G=LeadDataFm[1]*G-LeadG*Fm;  // leading terms in y should cancel...
    926 
    927   coef_Mat = coef(G,xy);
    928   Jnew=coef_Mat[2,1..ncols(coef_Mat)];
    929 
    930   // simplification of Jnew
    931 
    932   if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
    933                          // deformation parameters can be cutted off
    934   {
    935     Jnew=jet(Jnew,artin_bd-1);
    936   }
    937   Jnew=interred(Jnew);
    938   if (defined(artin_bd)) { Jnew=jet(Jnew,artin_bd-1); }
    939   J=J,Jnew;
    940 
    941   if (typ==1) // isEquising
    942   {
    943     if(ideal(nselect(J,1..no_b))<>0)
    944     {
    945       setring old_ring;
    946       option(set,ov);
    947       return(0);
    948     }
    949   }
    950 
    951   while (fertig<>soll and blowup<nrows(M[3]))
    952   {
    953     upper_bound_old=upper_bound;
    954     dbprint(i_print,"// Blowup Step "+string(blowup)+" completed");
    955     blowup=blowup+1;
    956 
    957     for (branch=1;branch<=number_of_branches;branch=branch+1)
    958     {
    959       Jnew=0;
    960 
    961       // First we check if the branch still has to be considered:
    962       if (branch==upper_bound_old[branch] and fertig[branch]<>1)
    963       {
    964         if (M[3][blowup-1,branch]==1 and
    965                ((B[blowup,branch]<>x and B[blowup,branch]<>y)
    966             or (blowup==nrows(M[3])) ))
    967         {
    968           fertig[branch]=1;
    969           dbprint(i_print,"// 1 branch finished");
    970         }
    971       }
    972 
    973       if (branch<=upper_bound_old[branch] and fertig[branch]<>1)
    974       {
    975         for (i=branch;i>=1;i--)
    976         {
    977           if (M[1][blowup-1,i]<>0)
    978           {
    979             m=M[1][blowup-1,i]; // multiplicity before blowup
    980             i=0;
    981           }
    982         }
    983 
    984         // we blow up the branch and take the strict transform:
    985         attrib(J,"isSB",1);
    986         bl_Map(branch)=reduce(bl_Map(branch),J);
    987         attrib(J,"isSB",0);
    988 
    989         phi=myRing,bl_Map(branch);
    990         F(branch)=phi(F(branch))/x^m;
    991 
    992         // simplify F
    993         attrib(Jnew,"isSB",1);
    994 
    995         F(branch)=reduce(F(branch),Jnew);
    996         attrib(Jnew,"isSB",0);
    997 
    998         m=M[1][blowup,branch]; // multiplicity after blowup
    999         Fm=m_Jet(F(branch),m); // homogeneous part of lowest degree
    1000 
    1001 
    1002         // we check for Fm=F[k]*...*F[k+s] where
    1003         //
    1004         //    F[j]=(y-b'(j)*x)^m(j), respectively F[j]=(-b'(j)*y+x)^m(j)
    1005         //
    1006         // according to the entries m(j)= M[3][blowup,j] and
    1007         //                          b'(j) mod m_A = B[blowup,j]
    1008         // computed from the HNE of the special fibre of the family:
    1009         G=1;
    1010         counter=branch;
    1011         k=upper_bound[branch];
    1012 
    1013         F_save=F(branch);
    1014 
    1015         while(counter<=k)
    1016         {
    1017           F(counter)=m_Jet(F_save,maxDeg[blowup,counter]);
    1018 
    1019           if (B[blowup,counter]<>x and B[blowup,counter]<>y)
    1020           {
    1021             G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]);
    1022             bl_Map(counter)=maxideal(1);
    1023             bl_Map(counter)[nvars(basering)]=
    1024                                         xy+(b(auxVar)+B[blowup,counter])*x;
    1025             bNodes[counter]=b(auxVar);
    1026             auxVar=auxVar+1;
    1027           }
    1028           else
    1029           {
    1030             if (B[blowup,counter]==x)
    1031             {
    1032               G=G*x^(M[3][blowup,counter]);  // branch has tangent x !!
    1033               F(counter)=swapXY(F(counter)); // will turn x to y for blow up
    1034               bl_Map(counter)=maxideal(1);
    1035               bl_Map(counter)[nvars(basering)]=xy;
    1036             }
    1037             else
    1038             {
    1039               G=G*y^(M[3][blowup,counter]); // tangent has to be y
    1040               bl_Map(counter)=maxideal(1);
    1041               bl_Map(counter)[nvars(basering)]=xy;
    1042             }
    1043             bNodes[counter]=0;
    1044           }
    1045           upper_bound[counter]=counter+M[2][blowup+1,counter]-1;
    1046           counter=counter+M[2][blowup+1,counter];
    1047         }
    1048         G=determine_coef(Fm)[1]*G-Fm;  // leading terms in y should cancel
    1049         coef_Mat = coef(G,xy);
    1050         Jnew=coef_Mat[2,1..ncols(coef_Mat)];
    1051         if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
    1052                               // deformation parameters can be cutted off
    1053         {
    1054           Jnew=jet(Jnew,artin_bd-1);
    1055         }
    1056 
    1057         // simplification of J
    1058         Jnew=interred(Jnew);
    1059 
    1060         J=J,Jnew;
    1061         if (typ==1) // isEquising
    1062         {
    1063           if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    1064           if(ideal(nselect(J,1..no_b))<>0)
    1065           {
    1066             setring old_ring;
    1067             option(set,ov);
    1068             return(0);
    1069           }
    1070         }
    1071       }
    1072     }
    1073     if (number_of_branches>=2)
    1074     {
    1075       J=interred(J);
    1076       if (typ==1) // isEquising
    1077       {
    1078         if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    1079         if(ideal(nselect(J,1..no_b))<>0)
    1080         {
    1081           setring old_ring;
    1082           option(set,ov);
    1083           return(0);
    1084         }
    1085       }
    1086     }
    1087   }
    1088 
    1089   // Computation for all equimultiple sections being trivial (I^s(f))
    1090   ideal Jtriv=J;
    1091   for (i=1;i<=no_b; i++)
    1092   {
    1093     if (reduce(b(i),std(bNodes))!=0){
    1094       Jtriv=subst(Jtriv,b(i),0);
    1095     }
    1096   }
    1097   Jtriv=std(Jtriv);
    1098 
    1099 
    1100 
    1101   dbprint(i_print,"// ");
    1102   dbprint(i_print,"// Elimination starts:");
    1103   dbprint(i_print,"// -------------------");
    1104 
    1105   poly gg;
    1106   int b_left=no_b;
    1107 
    1108   for (i=1;i<=no_b; i++)
    1109   {
    1110     attrib(J,"isSB",1);
    1111     gg=reduce(b(i),J);
    1112     if (gg==0)
    1113     {
    1114       b_left = b_left-1;  // another b(i) has to be 0
    1115     }
    1116     J = subst(J, b(i), gg);
    1117     attrib(J,"isSB",0);
    1118   }
    1119   J=simplify(J,10);
    1120   if (typ==1) // isEquising
    1121   {
    1122     if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    1123     if(ideal(nselect(J,1..no_b))<>0)
    1124     {
    1125       setring old_ring;
    1126       option(set,ov);
    1127       return(0);
    1128     }
    1129   }
    1130 
    1131   //new CL 11/06:  check in which equations b(k) appears and remove those b(k)
    1132   //               which appear in exactly one of the equations (by removing this
    1133   //               equation)
    1134   dbprint(i_print,"// ");
    1135   dbprint(i_print,"// Remove superfluous equations:");
    1136   dbprint(i_print,"// -----------------------------");
    1137   int Z,App_in;
    1138   ideal J_Tmp;
    1139   int ncJ=ncols(J);
    1140 
    1141   intmat Mdet[ncJ][1];
    1142   for (Z=1;Z<=ncJ;Z++){ Mdet[Z,1]=Z; }
    1143 
    1144   for (i=1;i<=no_b; i++)
    1145   {
    1146     ideal b_appears_in(i);            // Eintraege sind spaeter 1 oder 0
    1147     intmat b_app_in(i)[1][ncJ];       // Eintraege sind spaeter 1 oder 0
    1148     b_appears_in(i)[ncJ]=0;
    1149     J_Tmp = matrix(J)-subst(J,b(i),0);
    1150     for (Z=1; Z<=ncJ; Z++) {
    1151       if (J_Tmp[Z]<>0) {      // b(i) appear in J_Tmp[Z]
    1152         b_appears_in(i)[Z]=1;
    1153         b_app_in(i)[1,Z]=1;
    1154       }
    1155     }
    1156     if (size(b_appears_in(i))==1) { //b(i) appears only in one J_Tmp[Z]
    1157       App_in = (b_app_in(i)*Mdet)[1,1];  // determines Z
    1158       J[App_in]=0;
    1159       b_appears_in(i)[App_in]=0;
    1160       b_app_in(i)[1,App_in]=0;
    1161     }
    1162   }
    1163 
    1164   for (i=1;i<=no_b; i++)
    1165   {
    1166     if (size(b_appears_in(i))==1) { //b(i) appears only in one J_Tmp[Z]
    1167       App_in = (b_app_in(i)*Mdet)[1,1];  // determines Z
    1168       J[App_in]=0;
    1169       b_appears_in(i)[App_in]=0;
    1170       b_app_in(i)[1,Z]=1;
    1171       i=0;
    1172     }
    1173   }
    1174 
    1175   Jtriv = nselect(Jtriv,1..no_b);
    1176   ideal J_no_b = nselect(J,1..no_b);
    1177   if (size(J) > size(J_no_b))
    1178   {
    1179     dbprint(i_print,"// std computation started");
    1180     // some b(i) didn't appear in linear conditions and have to be eliminated
    1181     if (defined(artin_bd))
    1182     {
    1183       // first we make the ring smaller (removing variables, which are
    1184       // forced to 0 by J
    1185       list LL=make_ring_small(J);
    1186       ideal Shortmap=LL[2];
    1187       minPolyStr = "";
    1188       if (minpoly !=0)
    1189       {
    1190         minPolyStr = "minpoly = "+string(minpoly);
    1191       }
    1192       ordStr = "dp(" + string(b_left) + "),dp";
    1193       ideal qId = ideal(basering);
    1194 
    1195       helpStr = "ring Shortring = ("
    1196                 + charstr(basering) + "),("+ LL[1] +") , ("+ ordStr  +");";
    1197       execute(helpStr);
    1198       execute(minPolyStr);
    1199       // ring has changed to "Shortring"
    1200 
    1201       ideal MM=maxideal(artin_bd);
    1202       MM=subst(MM,x,0);
    1203       MM=subst(MM,y,0);
    1204       MM=simplify(MM,2);
    1205       dbprint(i_print-1,"// maxideal("+string(artin_bd)+") has "
    1206                          +string(size(MM))+" elements");
    1207       dbprint(i_print-1,"//");
    1208 
    1209       // we change to the qring mod m^artin_bd
    1210       // first, we have to check if we were in a qring when starting
    1211       ideal qId = imap(myRing, qId);
    1212       if (size(qId) == 0)
    1213       {
    1214          attrib(MM,"isSB",1);
    1215          qring QQ=MM;
    1216       }
    1217       else
    1218       {
    1219          qId=qId,MM;
    1220          qring QQ = std(qId);
    1221       }
    1222 
    1223       ideal Shortmap=imap(myRing,Shortmap);
    1224       map phiphi=myRing,Shortmap;
    1225 
    1226       ideal J=phiphi(J);
    1227       option(redSB);
    1228       J=std(J);
    1229       J=nselect(J,1..no_b);
    1230 
    1231       setring myRing;
    1232       // back to "myRing"
    1233 
    1234       J=nselect(J,1..no_b);
    1235       Jnew=imap(QQ,J);
    1236 
    1237       J=J,Jnew;
    1238       J=interred(J);
    1239       if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
    1240     }
    1241     else
    1242     {
    1243       J=std(J);
    1244       J=nselect(J,1..no_b);
    1245       if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
    1246     }
    1247   }
    1248 
    1249   dbprint(i_print,"// finished");
    1250   dbprint(i_print,"// ");
    1251 
    1252   minPolyStr = "";option(set,ov);
    1253   if (minpoly !=0)
    1254   {
    1255    minPolyStr = "minpoly = "+string(minpoly);
    1256   }
    1257 
    1258   kill HNEring;
    1259 
    1260   if (typ==1) // isEquising
    1261   {
    1262     if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    1263     if(J<>0)
    1264     {
    1265       setring old_ring;
    1266       option(set,ov);
    1267       return(0);
    1268     }
    1269     else
    1270     {
    1271       setring old_ring;
    1272       option(set,ov);
    1273       return(1);
    1274     }
    1275   }
    1276 
    1277   setring old_ring;
    1278   // we are back in the original ring
    1279 
    1280   if (npars(myRing)<>0)
    1281   {
    1282     ideal qIdeal = ideal(basering);
    1283     helpStr = "ring ESSring = ("
    1284                  + string(char(basering))+ "," + parstr(myRing) +
    1285                  ") , ("+ varstr(basering)+") , ("+ ordstr(basering) +");";
    1286     execute(helpStr);
     1ring ESSring = create_ring(s3, "("+ varstr(basering)+")", "("+ ordstr(basering) +")");
    12872    execute(minPolyStr);
    12883    // basering has changed to ESSring
Note: See TracChangeset for help on using the changeset viewer.