Changeset 6ea057 in git
 Timestamp:
 Sep 20, 2010, 12:38:34 PM (13 years ago)
 Branches:
 (u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
 Children:
 5d651dc1723c0ceef3e8c472352840a2c713d8a4
 Parents:
 cfd2e6ba5e3f99f46e106336cafaef280361dd61
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/tropical.lib
rcfd2e6 r6ea057 957 957 proc puiseuxExpansion (poly f,int n,list #) 958 958 "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; 959 ASSUME: f is a nonconstant 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; 962 964 the optional parameter # can be the string 'subst' 963 965 RETURN: list, where each entry of the list l describes the NewtonPuiseux … … 967 969 @* l[i][1] = is a ring 968 970 @* l[i][2] = int 969 @* l[i][3] = list971 @* l[i][3] = string 970 972 @* 971 973 @* WE NOW DESCRIBE THE LIST ENTRIES l[i] IN MORE DETAIL: … … 977 979 and LIFT[2] describes how the old parameter can be computed from the new one 978 980 @*  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 980 982 REMARK:  it is best to use the procedure displayPuiseuxExpansion to 981 983 display the result 982 984 @*  the procedure requires the @sc{Singular} procedure absPrimdecGTZ to be 983 985 present in the package primdec.lib 986 @*  if f is not squarefree it will be replaced by its squarefree part 984 987 EXAMPLE: example puiseuxExpansion; shows an example" 985 988 { 989 if (deg(f)<=0) 990 { 991 ERROR("The input polynomial is not allowed to be constant!"); 992 } 986 993 if (char(basering)!=0) 987 994 { … … 995 1002 { 996 1003 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."); 997 1008 } 998 1009 // if the ordering is not local, change to a local ordering … … 1000 1011 { 1001 1012 def GLOBALRING=basering; 1013 number mp=minpoly; 1002 1014 execute("ring LOCALRING=("+charstr(basering)+"),("+varstr(basering)+"),ds;"); 1003 1015 poly f=imap(GLOBALRING,f); 1016 minpoly=imap(GLOBALRING,mp); 1004 1017 } 1005 1018 // check if a substitution is necessary … … 1012 1025 } 1013 1026 } 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); 1050 1054 } 1051 1055 else 1052 1056 { 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; 1062 1143 } 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!"; 1073 1149 } 1074 1150 return(pexp); … … 1118 1194 variablen[j1]=string(var(j)); 1119 1195 } 1196 def BASERING=basering; 1120 1197 def LIFTRing=puiseux[1]; 1121 1198 int N=puiseux[2]; … … 1142 1219 ""; 1143 1220 "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 } 1153 1238 } 1154 1239 } … … 1157 1242 if (#[1]=="subst") 1158 1243 { 1244 setring BASERING; 1159 1245 ""; 1160 1246 "Substituting the expansion into the polynomial gives:"; 1161 "f="+puiseux[3] [size(puiseux[3])];1247 "f="+puiseux[3]; 1162 1248 } 1163 1249 } … … 5628 5714 { // we have to pass to a ring with two variables to compute the newtonpoly 5629 5715 def PRENPRing=basering; 5630 poly NPpoly=i[ size(i)];5716 poly NPpoly=i[1]; 5631 5717 ring NPRING=(0,@a),(x(1),t),ds; 5632 5718 poly NPpoly=imap(PRENPRing,NPpoly); … … 5860 5946 // ideal <t,x_1,...,[no x_lastvar],...,x_n> so that 5861 5947 // 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 5862 5951 /* 5863 5952 // DER NACHFOLGENDE TEIL IST MUELL UND WIRD NICHT MEHR GAMACHT … … 5896 5985 if (size(LI)==0) // if no power of t is in lead(I) (where @a is considered as a field element) 5897 5986 */ 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 } 5900 6001 LI=subst(I,var(nvars(basering)),0); 5901 6002 //size(deletedvariables)=anzahlvariablen(before elimination) … … 5949 6050 // all solutions with that component zero; 5950 6051 // 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 } 5952 6066 I=reduce(i,std(m)); // "remove" m from i 5953 6067 // test first, if i still is in the ideal <t,x_1,...,x_n>  if not, then … … 5959 6073 LI=subst(LI,var(j),0); 5960 6074 } 5961 if (size(LI)==0) // the entries of i have no tconstant part6075 if (size(LI)==0) // the entries of i have no constant part 5962 6076 { 5963 6077 // check now, if the leading ideal of i contains an element … … 6004 6118 for (j=anzahlvariablen1;j>=1;j) // the variable t is the last one !!! 6005 6119 { 6006 I=sat(ideal( var(j)+i),t)[1];6120 I=sat(ideal(std(var(j)+i)),t)[1]; 6007 6121 LI=subst(I,var(nvars(basering)),0); 6008 6122 //size(deletedvariables)=anzahlvariablen1(before elimination) … … 6225 6339 { 6226 6340 if (l!=1) // corresponds to x_(l1)  note t is the last variable 6227 { 6228 i[k]=subst itute(i[k],var(l1),(zeroneu[l]+var(l1))*t^(w[l]));6341 { 6342 i[k]=subst(i[k],var(l1),(zeroneu[l]+var(l1))*t^(w[l])); 6229 6343 } 6230 6344 else // corresponds to t 6231 6345 { 6232 i[k]=subst itute(i[k],var(nvars(basering)),var(nvars(basering))^(w[l]));6346 i[k]=subst(i[k],var(nvars(basering)),var(nvars(basering))^(w[l])); 6233 6347 } 6234 6348 }
Note: See TracChangeset
for help on using the changeset viewer.