Changeset 6ea057 in git


Ignore:
Timestamp:
Sep 20, 2010, 12:38:34 PM (14 years ago)
Author:
Thomas Markwig <keilen@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '3720ae8bfcff4a4649ee98a15552089151d2d59b')
Children:
5d651dc1723c0ceef3e8c472352840a2c713d8a4
Parents:
cfd2e6ba5e3f99f46e106336cafaef280361dd61
Message:
puiseuxExpansion und displayPuiseuxExpansion hinzugefuegt.


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/tropical.lib

    rcfd2e6 r6ea057  
    957957proc puiseuxExpansion (poly f,int n,list #)
    958958"USAGE:  puiseuxExpansion(f,n,#); f poly, n int, # list
    959 ASSUME:  f is a polynomial in two variables and the base field is either
    960          the field of rational numbers or a finite extension thereof;
    961          monomial ordering is assumed to be local;
     959ASSUME:  f is a non-constant polynomial in two variables which is not
     960         divisible by the first variable and which is squarefree
     961         as a power series over the complex numbers;
     962         the base field is either the field of rational numbers or a finite extension thereof;
     963         monomial ordering is assumed to be local;
    962964         the optional parameter # can be the string 'subst'
    963965RETURN:  list, where each entry of the list l describes the Newton-Puiseux
     
    967969@*             l[i][1] = is a ring
    968970@*             l[i][2] = int
    969 @*             l[i][3] = list
     971@*             l[i][3] = string
    970972@*
    971973@*       WE NOW DESCRIBE THE LIST ENTRIES l[i] IN MORE DETAIL:
     
    977979           and LIFT[2] describes how the old parameter can be computed from the new one
    978980@*       - if the option subst was set l[i][3] contains the polynomial where
    979            y has been substituted by y(t^{1/N})
     981           y has been substituted by y(t^{1/N}) as a string
    980982REMARK:  - it is best to use the procedure displayPuiseuxExpansion to
    981983           display the result
    982984@*       - the procedure requires the @sc{Singular} procedure absPrimdecGTZ to be
    983985           present in the package primdec.lib
     986@*       - if f is not squarefree it will be replaced by its squarefree part
    984987EXAMPLE:   example puiseuxExpansion;   shows an example"
    985988{
     989  if (deg(f)<=0)
     990  {
     991    ERROR("The input polynomial is not allowed to be constant!");
     992  }
    986993  if (char(basering)!=0)
    987994  {
     
    9951002  {
    9961003    ERROR("The base ring should depend on exactly two variables.");
     1004  }
     1005  if (f==(f / var(1))*var(1))
     1006  {
     1007    ERROR("The input polynomial is not allowed to be divisible by the first variable.");
    9971008  }
    9981009  // if the ordering is not local, change to a local ordering
     
    10001011  {
    10011012    def GLOBALRING=basering;
     1013    number mp=minpoly;
    10021014    execute("ring LOCALRING=("+charstr(basering)+"),("+varstr(basering)+"),ds;");
    10031015    poly f=imap(GLOBALRING,f);
     1016    minpoly=imap(GLOBALRING,mp);
    10041017  }
    10051018  // check if a substitution is necessary
     
    10121025    }
    10131026  }
    1014   // compute the Newton polygon
    1015   int jj,zw;
    1016   intvec w;
    1017   int ggteiler;
    1018   list NewtP=newtonpoly(f);
    1019   list tls;
    1020   list pexp;
    1021   // if the base field has a minimal polynomial change the base ring and the ideal
    1022   if (minpoly!=0)
    1023   {
    1024     poly mp=minpoly;
    1025     def OLDRING=basering;
    1026     execute("ring NEWRING=0,("+varstr(basering)+","+parstr(basering)+"),ds;");
    1027     ideal I=imap(OLDRING,mp),imap(OLDRING,f);
    1028   }
    1029   else
    1030   {
    1031     ideal I=f;
    1032   }
    1033   // for each facet of the Newton polygon compute the corresponding parametrisations
    1034   // using tropicalLifting with the option "puiseux" which avoids gfan
    1035   for (jj=1;jj<=size(NewtP)-1;jj++)
    1036   {
    1037     w=NewtP[jj]-NewtP[jj+1];
    1038     ggteiler=gcd(w[1],w[2]);
    1039     zw=w[1]/ggteiler;
    1040     w[1]=w[2]/ggteiler;
    1041     w[2]=zw;
    1042     // if we have introduced a third variable for the parameter, then w needs a third component
    1043     if (nvars(basering)==3)
    1044     {
    1045       w[3]=0;
    1046     }
    1047     if (noResubst==0)
    1048     {
    1049       tls=tropicalLifting(I,w,n,"findAll","puiseux");
     1027  // replace f by its squarefree part
     1028  f=squarefree(f);
     1029  // check if var(1) or var(2) divide f
     1030  int coordinatebranchtwo;
     1031  if (f==(f / var(2))*var(2))
     1032  {
     1033    coordinatebranchtwo=1;
     1034    f=f/var(2);
     1035  }
     1036  int jj;
     1037  // if f is now constant then we should skip the next part
     1038  if (deg(f)!=0)
     1039  {
     1040    // compute the Newton polygon
     1041    int zw;
     1042    intvec w;
     1043    int ggteiler;
     1044    list NewtP=newtonpoly(f);
     1045    list tls;
     1046    list pexp;
     1047    // if the base field has a minimal polynomial change the base ring and the ideal
     1048    if (minpoly!=0)
     1049    {
     1050      poly mp=minpoly;
     1051      def OLDRING=basering;
     1052      execute("ring NEWRING=0,("+varstr(basering)+","+parstr(basering)+"),ds;");
     1053      ideal I=imap(OLDRING,mp),imap(OLDRING,f);
    10501054    }
    10511055    else
    10521056    {
    1053       tls=tropicalLifting(I,w,n,"findAll","puiseux","noResubst");
    1054     }
    1055     pexp=pexp+tls;
    1056   }
    1057   // kill rings that are no longer needed
    1058   if (defined(NEWRING))
    1059   {
    1060     setring OLDRING;
    1061     kill NEWRING;
     1057      ideal I=f;
     1058    }
     1059    // for each facet of the Newton polygon compute the corresponding parametrisations
     1060    // using tropicalLifting with the option "puiseux" which avoids gfan
     1061    for (jj=1;jj<=size(NewtP)-1;jj++)
     1062    {
     1063      w=NewtP[jj]-NewtP[jj+1];
     1064      ggteiler=gcd(w[1],w[2]);
     1065      zw=w[1]/ggteiler;
     1066      w[1]=w[2]/ggteiler;
     1067      w[2]=zw;
     1068      // if we have introduced a third variable for the parameter, then w needs a third component
     1069      if (nvars(basering)==3)
     1070      {
     1071        w[3]=0;
     1072      }
     1073      if (noResubst==0)
     1074      {
     1075        tls=tropicalLifting(I,w,n,"findAll","puiseux");
     1076      }
     1077      else
     1078      {
     1079        tls=tropicalLifting(I,w,n,"findAll","puiseux","noResubst");
     1080      }
     1081      pexp=pexp+tls;
     1082    }
     1083    // kill rings that are no longer needed
     1084    if (defined(NEWRING))
     1085    {
     1086      setring OLDRING;
     1087      kill NEWRING;
     1088    } 
     1089    if (defined(GLOBALRING))
     1090    {
     1091      setring GLOBALRING;
     1092      kill LOCALRING;
     1093    } 
     1094    // remove the third entry in the list of parametrisations since we know
     1095    // that none of the exponents in the parametrisation will be negative
     1096    for (jj=1;jj<=size(pexp);jj++)
     1097    {
     1098      pexp[jj]=delete(pexp[jj],3);
     1099      pexp[jj][3]=pexp[jj][3][size(pexp[jj][3])]; // replace the list in pexp[jj][3] by its first entry
     1100    }
     1101  }
     1102  else // if f was reduced to a constant we still have to introduce the list pexp
     1103  {
     1104    list pexp;
     1105  }
     1106  // we have to add the parametrisations for the branches var(1) resp. var(2) if
     1107  // we removed them
     1108  def BASERING=basering;
     1109  if (coordinatebranchtwo==1)
     1110  {
     1111    ring brring=0,t,ds;
     1112    ideal LIFT=0;
     1113    export(LIFT);
     1114    setring BASERING;
     1115    list pexpentry;
     1116    pexpentry[1]=brring;
     1117    pexpentry[2]=1;
     1118    pexpentry[3]="0";
     1119    pexp[size(pexp)+1]=pexpentry;
     1120    kill brring;
     1121  }
     1122  // check if all branches have been found
     1123  int numberofbranchesfound;
     1124  for (jj=1;jj<=size(pexp);jj++)
     1125  {
     1126    def countring=pexp[jj][1];
     1127    setring countring;
     1128    if (minpoly==0)
     1129    {
     1130      numberofbranchesfound=numberofbranchesfound+1;
     1131    }
     1132    else
     1133    {
     1134      poly mp=minpoly;
     1135      ring degreering=0,a,dp;
     1136      poly mp=imap(countring,mp);
     1137      numberofbranchesfound=numberofbranchesfound+deg(mp);
     1138      setring countring;
     1139      kill degreering;
     1140    }
     1141    setring BASERING;
     1142    kill countring;
    10621143  } 
    1063   if (defined(GLOBALRING))
    1064   {
    1065     setring GLOBALRING;
    1066     kill LOCALRING;
    1067   } 
    1068   // remove the third entry in the list of parametrisations since we know
    1069   // that none of the exponents in the parametrisation will be negative
    1070   for (jj=1;jj<=size(pexp);jj++)
    1071   {
    1072     pexp[jj]=delete(pexp[jj],3);
     1144  // give a warning if not all branches have been found
     1145  if (numberofbranchesfound!=ord(subst(f,x,0))+coordinatebranchtwo)
     1146  {
     1147    "!!!! WARNING: The number of terms computed in the Puiseux expansion were";
     1148    "!!!!          not enough to find all branches of the curve singularity!";
    10731149  }
    10741150  return(pexp);
     
    11181194      variablen[j-1]=string(var(j));
    11191195    }
     1196    def BASERING=basering;
    11201197    def LIFTRing=puiseux[1];
    11211198    int N=puiseux[2];
     
    11421219    "";
    11431220    "The expansion has the form:";
    1144     for (j=1;j<=size(LIFT);j++)
    1145     {
    1146       if (ncols(LIFT)==size(variablen))
    1147       {
    1148         variablen[j]+"="+displaypoly(LIFT[j],N,0,1);
    1149       }
    1150       else
    1151       {
    1152         "var("+string(j)+")="+displaypoly(LIFT[j],N,0,1);
     1221    // treat the case LIFT==0 separately
     1222    if (size(LIFT)==0)
     1223    {
     1224      variablen[1]+"=0";
     1225    }
     1226    else
     1227    {
     1228      for (j=1;j<=size(LIFT);j++)
     1229      {
     1230        if (ncols(LIFT)==size(variablen))
     1231        {
     1232          variablen[j]+"="+displaypoly(LIFT[j],N,0,1);
     1233        }
     1234        else
     1235        {
     1236          "var("+string(j)+")="+displaypoly(LIFT[j],N,0,1);
     1237        }
    11531238      }
    11541239    }
     
    11571242      if (#[1]=="subst")
    11581243      {
     1244        setring BASERING;
    11591245        "";
    11601246        "Substituting the expansion into the polynomial gives:";
    1161         "f="+puiseux[3][size(puiseux[3])];
     1247        "f="+puiseux[3];
    11621248      }
    11631249    }
     
    56285714            { // we have to pass to a ring with two variables to compute the newtonpoly
    56295715              def PRENPRing=basering;
    5630               poly NPpoly=i[size(i)];
     5716              poly NPpoly=i[1];
    56315717              ring NPRING=(0,@a),(x(1),t),ds;
    56325718              poly NPpoly=imap(PRENPRing,NPpoly);
     
    58605946    // ideal <t,x_1,...,[no x_lastvar],...,x_n> so that
    58615947    // there is a solution which has strictly positive valuation
     5948    // ADDENDUM:
     5949    // if i (without m) has only one polynomial, then we can divide
     5950    // i by t as long as possible to compute the saturation with respect to t
    58625951/*
    58635952    // DER NACHFOLGENDE TEIL IST MUELL UND WIRD NICHT MEHR GAMACHT
     
    58965985    if (size(LI)==0) // if no power of t is in lead(I) (where @a is considered as a field element)
    58975986*/
    5898     I=reduce(sat(ideal(var(lastvar)+i),t)[1],std(m)); // get rid of the minimal
    5899                                                       // polynomial for the test
     5987    if (size(i)-size(m)!=1)
     5988    {
     5989      I=reduce(sat(std(ideal(var(lastvar)+i)),t)[1],std(m)); // get rid of the minimal
     5990                                                             // polynomial for the test
     5991    }
     5992    else
     5993    {
     5994      I=subst(i,var(lastvar),0);
     5995      while (subst(I[1],t,0)==0)
     5996      {
     5997        I[1]=I[1]/t;
     5998      }
     5999      I=reduce(I,std(m));
     6000    }
    59006001    LI=subst(I,var(nvars(basering)),0);
    59016002    //size(deletedvariables)=anzahlvariablen(before elimination)
     
    59496050    // all solutions with that component zero;
    59506051    // the checking is then done as above
    5951     i=std(sat(i,var(lastvar))[1]);
     6052    // ADDENDUM
     6053    // if the ideal i (without m) is generated by a single polynomial
     6054    // then we saturate by successively dividing by the variable
     6055    if (size(i)-size(m)==1)
     6056    {
     6057      while (subst(i[1],var(lastvar),0)==0)
     6058      {
     6059        i[1]=i[1]/var(lastvar);
     6060      }
     6061    }
     6062    else
     6063    {
     6064      i=std(sat(i,var(lastvar))[1]);
     6065    }
    59526066    I=reduce(i,std(m)); // "remove" m from i
    59536067    // test first, if i still is in the ideal <t,x_1,...,x_n> -- if not, then
     
    59596073      LI=subst(LI,var(j),0);
    59606074    }
    5961     if (size(LI)==0) // the entries of i have not constant part
     6075    if (size(LI)==0) // the entries of i have no constant part
    59626076    {     
    59636077      // check now, if the leading ideal of i contains an element
     
    60046118    for (j=anzahlvariablen-1;j>=1;j--)  // the variable t is the last one !!!
    60056119    {
    6006       I=sat(ideal(var(j)+i),t)[1];
     6120      I=sat(ideal(std(var(j)+i)),t)[1];
    60076121      LI=subst(I,var(nvars(basering)),0);
    60086122      //size(deletedvariables)=anzahlvariablen-1(before elimination)
     
    62256339      {
    62266340        if (l!=1) // corresponds to  x_(l-1) --  note t is the last variable
    6227         {       
    6228           i[k]=substitute(i[k],var(l-1),(zeroneu[l]+var(l-1))*t^(-w[l]));
     6341        {
     6342          i[k]=subst(i[k],var(l-1),(zeroneu[l]+var(l-1))*t^(-w[l]));
    62296343        }
    62306344        else // corresponds to t
    62316345        {
    6232           i[k]=substitute(i[k],var(nvars(basering)),var(nvars(basering))^(-w[l]));
     6346          i[k]=subst(i[k],var(nvars(basering)),var(nvars(basering))^(-w[l]));
    62336347        }
    62346348      }
Note: See TracChangeset for help on using the changeset viewer.