# source:git/Singular/LIB/arcAtPoint.lib@c1d97b

spielwiese
Last change on this file since c1d97b was c1d97b, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes/cremer: initial commit git-svn-id: file:///usr/local/Singular/svn/trunk@9126 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 16.2 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id: arcAtPoint.lib,v 1.1 2006-05-11 13:03:59 Singular Exp \$";
3//-*- mode:C++;-*-
4
5category="Singularities";
6info="
7LIBRARY:  arcAtPoint.lib        Truncations of arcs at a singular point
9
10OVERVIEW: An arc is given by a power series in one variable, say t, and
11          truncating it at a positive integer i means cutting
12          the t-powers > i. The set of arcs truncated at order
13          <bound> is denoted Tr(i). An algorithm for computing
14          these sets (which happen to be constructible) is given in
15          [Lejeune-Jalabert, M.: Courbes trac\'ees sur un germe
16          d'hypersurface, American Journal of Mathematics, 112 (1990)].
17          Our procedures for computing the locally closed sets contributing
18          to the set of truncations rely on this algorithm.
19
20PROCEDURES:
21    nashmult(f,bound);         determines locally closed sets relevant
22                               for computing truncations of arcs over a
23                               hypersurface with isolated singularity
24                               defined by f. The sets are given by two
25                               ideals specifying relations between
26                               coefficients of power series in t. One
27                               of the ideals defines an open set, the
28                               other one the complement of a closed
29                               set within the open one.
30                               We consider only coefficients up to
31                               t^<bound>.
32                               Moreover, the sequence of Nash
33                               Multiplicities of each set is
34                               displayed
35    removepower(I);            modifies the ideal I such that the
36                               algebraic set defined by it remains
37                               the same: removes powers of variables
38    idealsimplify(I,maxiter);  further simplification of I in the
39                               above sense: reduction with other
40                               elements of I. The positive integer
41                               <maxiter> gives a bound to the number of
42                               repetition steps
43    equalJinI(I,J);            tests if two ideals I and J are equal
44                               under the assumption that J is
45                               contained in I. Returns 1 if this is
46                               true and 0 otherwise
47    ";
48
49LIB "ring.lib";
50LIB "general.lib";
51LIB "standard.lib";
52LIB "sing.lib";
53
54
55//////////////////////////////////////////////////////////////////////
56//                                                                  //
57//                  RELEVANT LOCALLY CLOSED SETS                    //
58//                                                                  //
59//////////////////////////////////////////////////////////////////////
60
61
62proc nashmult (poly f, int bound)
63"USAGE:   nashmult(f,bound); f polynomial, bound positive integer
64CREATE:  allsteps:
65@format
66         a list containing all relevant locally closed sets
67         up to order <bound> and their sequences of
68         Nash Multiplicities
69@end format
70         setstep:
71@format
72         list of relevant locally closed sets
73         obtained from sequences of length bound+1
74@end format
75RETURN:  ring, original basering with additional
76         variables t and coeffiecients up to t^<bound>
77EXAMPLE: example nashmult; shows an example"
78
79{
80 // Make sure that only admissible  parameters are entered
81
82 if (bound<1)
83 {
84   ERROR("Integer parameter must be positive!");
85 }
86
87 // General setup, declarations and initialization...
88
89 int k,s,step,loop;                   // loop variables
90 int pos;                             // position parameter
91 int countall;                        // position parameter
92 list allsteps;                       // saves results from all
93                                      // steps
94 def r=basering;
95 int startvar=nvars(basering);
96 intvec control=order(f);             // initialize
97 def R=ringchange(bound+1);           // ring in which result lies
98
99
100 setring R;                           // make it basering
101 ideal I0;
102 list init=control,I0,I0;
103 list setstep=insert(setstep,init);   // stores Nash multiplicities
104 kill I0;                             // and underlying sets (given
105 kill init;                           // that the sets are not
106                                       // empty),will be a list of
107                                      // lists, each of which
108                                      // containing an intvec and
109                                      // two ideals
110
111
112  // consider all possible sequences of multiplicities<=<bound>:
113
114 for(step=2;step<=bound+1;step++)
115 {
116   list setsteptmp;                    // temporary variable, local
117                                       // to this loop
118   int count;                          // position parameter
119   int setsize=size(setstep);
120
121   setring r;
122   def rplug=pluginCoeffs(f,step);     // substitute variables in f
123                                       // by polynomials in t of
124                                       // degree <step>
125   setring R;
126   ideal coe=imap(rplug,resultcoe);    // gives the t-coefficients
127   kill rplug;
128
129
130   // consider all sequences of lenght <step-1> giving rise to a
131   // family...
132
133   for(loop=1;loop<=setsize;loop++)
134   {
135     control=setstep[loop][1];       // initialization. <control>
136     int sizecontrol=size(control);  // describes Nash mutiplicities
137     ideal gprev=setstep[loop][3];   // in this step, <fprev> and
138     ideal fprev=setstep[loop][2];   // <gprev> the already obtained
139                                     // relations
140
141     ideal actcoe=reduce(coe,slimgb(fprev));  // take into account
142                                              // existing relations
143
144     pos=1;                          // determine first nonzero
145     while(actcoe[pos]==0)           // t-coefficient
146     {
147       pos++;
148     }
149
150     ideal fset;                     // will store relations
151     ideal gset;                     // defining the
152                                     // constructible sets
153
154     int m0=control[sizecontrol];    // bounds the computations
155
156     // consider all possible sequences arising from <control>...
157
158     control=control,0;
159     for (s=1;s<=m0;s++)
160     {
161       control[step]=control[step]+1;  // the next possible sequence
162                                       // of multiplicities
163
164       fset=fset,actcoe[pos+s-1],gset; // add new condition
165       gset=gset,fset;
166
167       for(k=0;k<startvar;k++)         // additional conditions for
168       {                               // <gset>:
169         ideal coevar=coeffs(actcoe[pos+s],
170         var(startvar+1+step+k*(bound+1)));
171         int coesize=ncols(coevar);
172         if (coesize>=2)               // add coeff. of nonconstant
173         {                             // terms in "highest"
174           gset=gset,coevar[2..coesize]; // variables
175         }
176         kill coevar;
177         kill coesize;
178       }
179       fset=fprev,fset;                // add obtained conditions
180       gset=fprev,gset;                // to the existing ones...
181       gset=idealsimplify(gset,1000);  // ...and simplify
182       fset=idealsimplify(fset,1000);
183
184
185        // if we have found a nontrivial component...
186
187       if (control[step-1]==1)
188       {
189         list comp=control,fset,gset;              // ...add it and
190         setsteptmp=insert(setsteptmp,comp,count); // multiplicity
191         count++;
192         kill comp;
193       }
194       else
195       {
196         if (equalJinI(gset,fset)==0)
197         {
198           list comp=control,fset,gset;             // ...add it and
199           setsteptmp=insert(setsteptmp,comp,count);// multiplicity
200           count++;
201           kill comp;
202         }
203       }
204     }
205     kill fset,gset,actcoe,sizecontrol,fprev,gprev,m0;
206  }
207  setstep=setsteptmp;
208
209  allsteps=insert(allsteps,setstep,countall);   // add results from
210  countall++;                                   // this step
211
212  kill setsteptmp,count,coe,setsize;
213 }
214 export(setstep);
215 export(allsteps);
216 return(R);
217}
218example
219{
220  "EXAMPLE:"; echo=2;
221  ring r=0,(x,y,z),dp;
222  poly f=z4+y3-x2;
223  def R=nashmult(f,2);
224  setring R;
225  allsteps;
226}
227
228
229
230//////////////////////////////////////////////////////////////////////
231//                                                                  //
232//               SUBSTITUE VARIABLES BY POLYNOMIALS                 //
233//                                                                  //
234//////////////////////////////////////////////////////////////////////
235
236
237static proc pluginCoeffs (poly f,int i)
238"USAGE:   pluginCoeffs(f,i); f polynomial, i integer
239CREATE:  matrix, the t-coefficients obtained by replacing each
240         variable of f by a polynomial in t
241RETURN:  ring, corresponds to f and i in the sense that it is
242         the original ring with added variables t and
243         t-coefficients up to t^<bound> "
244{
245  int startvar=nvars(basering);
246  def r=basering;
247  def R=ringchange(i);              // changes the ring
248  setring R;                        // makes it new basering;
249  ideal I=tpolys(i,startvar);
250  map h=r,I;                        // define map
251  ideal resultplug=h(f);
252  matrix resultcoe=coeffs(resultplug[1],t);
253  export(resultcoe);                // keep accessible
254  export(resultplug);
255  return(R);                        //
256}
257
258//////////////////////////////////////////////////////////////////////
259static proc tpolys (int i,int k)
260"USAGE:   tpolys(i,k); i,k integer
261RETURN:  ideal, generated by k polynomials in t of degree i
262         of the form a(1)*t+..+a(i)*t^i
263NOTE:    called from pluginCoeffs"
264{
265  int s,t;
266  int v;
267  poly sum;
268  ideal I;
269  for(t=1;t<=k;t++)
270  {
271    v=(t-1)*i;
272    for(s=1;s<=i;s++)
273    {
274      sum=sum+var(1+k+v+s)*var(k+1)^s;
275    }
276    I[t]=sum;
277    sum=0;
278  }
279  return(I);
280}
281
282
283//////////////////////////////////////////////////////////////////////
284//                                                                  //
285//                  CONSTRUCTING THE RESULT RING                    //
286//                                                                  //
287//////////////////////////////////////////////////////////////////////
288
289
290static proc ringchange (int i)
291"USAGE:   ringchange(i); i integer
292RETURN:  ring, extends basering by variables t and
293          #(variables of basering)*i new variables"
294{
295  def R=changevar(""+varstr(basering)+",t,"
296  +variables(nvars(basering),i)+"");
297  return(R);
298}
299
300/////////////////////////////////////////////////////////////////////
301static proc variables (int k,int i)
302"USAGE:   variables(k,i); k,i integer
303ASSUME:  1<=k<=26, i>=1
304RETURN:  string of names of variables added in ringchange
305NOTE:    called from ringchange, we use this procedure to
306          obtain a nice shape of the ring created "
307{
308  list l;
309  int s,u;
310  string str;
311  for (u=1;u<=k;u++)
312  {
313    for (s=1;s<=i;s++)
314    {
315      str=""+atoz(u)+"("+string(s)+")"; // creates new variables
316      l[(u-1)*i+s]=str;                 // saves them in a list
317    }
318  }
319  string str1=string(l);               // makes a string of the
320  return(str1);                        // list (needed for change
321}                                      // of ring)
322
323//////////////////////////////////////////////////////////////////////
324static proc atoz (int n)
325"USAGE:   atoz(n); n integer
326ASSUME:  1<=n<=26
327RETURN:  string, the nth letter of the alphabet"
328
329 string s="ring r=0,("+A_Z("a",n)+"),ds;";
330 execute(s);
331 return (string(var(n)));
332}
333
334
335//////////////////////////////////////////////////////////////////////
336//                                                                  //
337//                  AUXILIARY PROCEDURES                            //
338//                                                                  //
339//////////////////////////////////////////////////////////////////////
340
341
342static proc order (poly f)
343"USAGE:   order(f); f polynomial
344RETURN:  int i, such that f_i is the smallest (in terms of degree)
345         non-zero homogeneous part
346NOTE:    is designed for ordering dp"
347{
348 int k=deg(f);
349 int i;
350 for(i=1;i<=k;i++)
351 {
352  if(jet(f,i)!=0)
353  {
354    return(i);
355  }
356 }
357}
358
359//////////////////////////////////////////////////////////////////////
360static proc modd (poly f, poly g)
361"USAGE:   modd(f,g); f,g polynomials
362RETURN:  poly, f mod g division with remainder
363NOTE:    called from idealsimplify where it is used to modify
364          the generating set of an ideal"
365{
366  poly result=f-(f/g)*g;
367  return(result);
368}
369
370//////////////////////////////////////////////////////////////////////
371proc removepower (ideal I)
372"USAGE:   removepower(I); I ideal
374RETURN:  ideal defining the same zeroset as I: if any  generator
375         of I is a power of one single variable, replace it by the
376         respective variable
377EXAMPLE: example removepower; shows an example"
378{
379  int i,j;
380  int divisornumber=0;
381  int pos;
382  I=simplify(I,6);            // remove 0 and multiple generators
383  for(j=1;j<=ncols(I);j++)
384  {
385    if(size(I[j])==1)         // test if generators are powers
386    {                         // of variables...
387      for(i=1;i<=nvars(basering);i++)
388      {
389        if(modd(I[j],var(i))==0)
390        {
391          divisornumber++;
392          pos=i;
393        }
394      }
395    }
396    if(divisornumber==1)      // ...if so, replace by variable
397    {
398      I[j]=var(pos);
399    }
400    divisornumber=0;
401  }
402  return(I);
403}
404example
405{
406  "EXAMPLE:"; echo=2;
407   ring r=0,(x,y,z),dp;
408   ideal I = x3,y+z2-x2;
409   I;
410
411   removepower(I);
412}
413
414//////////////////////////////////////////////////////////////////////
415proc idealsimplify (ideal I, int maxiter)
416"USAGE:   idealsimplify(I); I ideal
417ASSUME:  procedure is stable for sufficiently large maxiter
418RETURN:  ideal defining the same zeroset as I: replace   generators
419         of I by the generator modulo other generating elements
420EXAMPLE: example idealsimplify; shows an example "
421{
422  if(maxiter<1)
423    {ERROR("The integer argument has to be positive!")}
424  ideal comp;
425  int iteration;
426  int changed=0;
427  int i,j,ci,n,cols;
428  for(iteration=0;iteration<maxiter;iteration++)
429  {
430    comp=I;
431    n=ncols(I);
432    for(j=2;j<=n;j++)         // reduce with lower elements
433    {
434      for(i=1;i<j;i++)
435      {
436        if(I[i]!=0)
437        {
438          I[j]=modd(I[j],I[i]);
439        }
440      }
441    }
442    I=simplify(removepower(I),7);
443
444    kill n;
445    int n=ncols(I);
446    for(j=n-1;j>=1;j--)       // reduce with higher elements
447    {
448      for(i=n;i>j;i--)
449      {
450        if(I[i]!=0)
451        {
452          I[j]=modd(I[j],I[i]);
453        }
454      }
455     }
456    I=simplify(removepower(I),7);
457
458    if (ncols(I)==ncols(comp))      //check if I has changed
459    {
460      cols=ncols(I);
461      changed=0;
462      for(ci=1;ci<=cols;ci++)
463      {
464        if (I[ci]!=comp[ci])
465        {
466          changed=1;
467          break;
468        }
469      }
470      if (changed==0)
471       break;
472    }
473  }
474  return(I);
475}
476example
477{
478  "EXAMPLE:"; echo=2;
479  ring r=0,(x,y,z),dp;
480  ideal I = x3,y+z2-x2;
481  I;
482
483  idealsimplify(I,10);
484}
485
486//////////////////////////////////////////////////////////////////////
487proc equalJinI (ideal I, ideal J)
488"USAGE:   equalJinI(I,J); (I,J ideals)
489ASSUME:  J contained in I and both I and J have been processed
490         with idealsimplify before
492RETURN:  1, if I=J, 0 otherwise
493EXAMPLE: example equalJinI; shows an example"
494{
495  int col=ncols(I);
496  J=slimgb(J);
497  int k;
498  for(k=1;k<=col;k++)
499  {
500    if(reduce(I[k],J)!=0)
501    { return(0);}
502  }
503  return(1);
504}
505example
506{
507  "EXAMPLE:"; echo=2;
508  ring r=0,(x,y,z),dp;
509  ideal I = x,y+z2;
510  ideal J1= x;
511  ideal J2= x,y+z2;
512
513  equalJinI(I,J1);
514  equalJinI(I,J2);
515}
Note: See TracBrowser for help on using the repository browser.