source: git/Singular/LIB/lejeune.lib @ 1e1ec4

spielwiese
Last change on this file since 1e1ec4 was 1e1ec4, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Updated LIBs according to master add: new LIBs from master fix: updated LIBs due to minpoly/(de)numerator changes fix: -> $Id$ fix: Fixing wrong rebase of SW on master (LIBs)
  • Property mode set to 100644
File size: 17.1 KB
Line 
1//-*- mode:C++;-*-
2version="$Id$";
3category="??";
4
5info="
6LIBRARY: lejeune.lib     Arc space computations
7AUTHOR:  Nadine Cremer,  nadine.cremer@gmx.de
8PROCEDURES:
9    fgset(f,H);          returns simultaneously the sets F and G corresponding
10                         to H as described by M. Lejeune
11    trunc(f,i);          returns the set Tr(i) as described by M. Lejeune
12    ";
13
14
15LIB "ring.lib";
16LIB "general.lib";
17LIB "primdec.lib";
18LIB "standard.lib";
19LIB "sing.lib";
20
21
22
23//////////////////////////////////////////////////////////////////////
24proc trunc (poly f, int i)
25"USAGE:   trunc(f,i); (f polynomial, i integer)
26 CREATE:  list, contains lists, each one consisting of two ideals:
27          the first one giving equations, the second one giving
28          inequations for a part of Tr(i). All of them together give a
29          complete description of Tr(i)
30 RETURN:  ring, corresponds to f and i, i.e. plugging in a polynomial
31          in t of degree i for each variable of f yields a ring whose
32          variables are the original variables of f, t and the
33          according t-coefficients
34 EXAMPLE: example trunc; shows an example"
35{
36 def r=basering;
37
38 // we make sure that we obtain Tr(i), which requires
39 // computations up to m*i, where m is minimal s.th.
40 // x(j)^m in J(f) for each variable x(j) involved in f
41 int m=minpower(f);
42 int mi=m*i;
43 int k;
44 int l=order(f);
45 int s0=mi+1;
46 int s=l*(mi+1);
47 int z1=nvars(r)*(i+1)+2;
48 int z2=nvars(r)*k+1;
49 intvec H=l;
50
51 // initialization of an intvec H of size m*i
52 for(k=1;k<=mi;k++)
53 {
54   H[k+1]=1;
55 }
56 // this is the ring in which result lies
57 def R1=ringchange(i);
58 setring R1;
59 ideal I,J;
60 list intersec=I,J;
61
62 // will save the results:
63 list fresult,gresult;
64 list result;
65
66 // consider all possible H's
67 while(sum(H)<=s)
68 {
69   setring r;
70   def tmp=fgset(f,H);
71   setring R1;
72   intersec=imap(tmp,fgresult);
73   kill tmp;
74
75   // simplifications
76   intersec[1]=simplifymodd(interred(intersec[1]));
77   intersec[2]=simplifymodd(interred(intersec[2]));
78   //option(redSB);
79   //intersec[1]=std(intersec[1]);
80   //intersec[2]=std(intersec[2]);
81   intersec[1]=simplifymodd(intersec[1]);
82   intersec[2]=simplifymodd(intersec[2]);
83   intersec[1]=simplifymodd(intersec[1]);
84   intersec[2]=simplifymodd(intersec[2]);
85
86   // remove lists which contain
87   // the same ideal twice and
88   // therefore define the empty
89   // set
90   if(equalitytest(intersec[1],intersec[2])==1)
91     {
92       H=Hnew(H);
93       continue;
94   }
95   intersec[1]=radical(intersec[1]);
96   intersec[2]=radical(intersec[2]);
97
98   // remove lists which contain
99   // the same ideal twice and
100   // therefore define the empty
101   // set
102   if(equalitytest(intersec[1],intersec[2])==1)
103     {
104       H=Hnew(H);
105       continue;
106     }
107   result=insert(result,intersec);
108   H=Hnew(H);
109 }
110 // output:
111 int u=size(result);
112 newline;
113 " We obtain the following sets of equations
114 and inequations for Tr(" +string(i)+"). In order to be
115 contained in Tr(" +string(i)+"), a point has to fulfill
116 the conditions of one of those set, i.e. it has to be
117 in the zero-set of the first ideal without
118 the zero-set of the second ideal.";
119 newline;
120 string ending;
121
122 for(k=1;k<=u;k++)
123 {
124   if((k mod 10==1) && (k mod 100!=11))
125   {ending="st";}
126   else
127   {
128     if((k mod 10==2) && (k mod 100!=12))
129     {ending="nd";}
130     else
131     {
132       if((k mod 10==3) && (k mod 100!=13))
133       {ending="rd";}
134       else
135       {ending="th";}
136     }
137   }
138   print(string(k)+ending+" set of equations:"+newline);
139   print(result[k]);
140   newline;
141 }
142 // return the ring
143 return(R1);
144}
145
146example
147{
148  "EXAMPLE:"; echo=2;
149  def r=basering;
150  poly f=y2-x3;
151  def R=trunc(f,3);
152  setring R;
153  print(result);
154  setring r;
155  def R1=trunc(f,6);
156  setring R1;
157  print(result);
158}
159
160
161
162
163
164//////////////////////////////////////////////////////////////////////
165//                                                                  //
166//             COMPUTATION OF F AND G SIMULTANEOUSLY                //
167//                                                                  //
168//////////////////////////////////////////////////////////////////////
169
170//////////////////////////////////////////////////////////////////////
171proc fgset (poly f,intvec H)
172"USAGE:   fgset(f,H);  f polynomial, H integer vector
173 CREATE:  list, consists two ideals, the first one giving equations,
174          the second one giving inequations to be satisfied by the set
175          corresponding to H
176 RETURN:  ring, corresponds to f and size(H)-1, i.e. plugging in a
177          polynomial in t of degree size(H)-1 for each variable of f
178          yields a ring whose variables are the original variables of f,
179          t and the according t-coefficients
180 EXAMPLE: example fgset; shows an example"
181{
182  def r=basering;
183  int p;
184  int m0=order(f);
185  int b=size(H);
186  if(H[1]!=m0)                       // input admissible?!
187    {
188      ERROR("H[1]=ord(f) necessary");
189    }
190  for(p=1;p<b;p++)
191    {
192      if(H[p]<H[p+1])
193      {
194        ERROR("Inadmissible input, H[1]<=...<=H[b] necessary");
195      }
196    }
197
198  def R=ringchange(b-1);
199  setring R;
200  list l;
201  ideal fresult,gresult;
202  list fgresult;
203  for(p=2;p<=b;p++)
204  {
205    setring r;
206    def tmp=formaldiff(f,intvec(H[1..p]));
207    setring R;
208    l=imap(tmp,resultdiff);
209    kill tmp;
210    fresult=fresult,l[1];
211  }
212  gresult=fresult;                     // use computation of f for g
213
214  setring r;                           // last step, special for G
215  def tmp=formaldiff(f,H);;
216  setring R;
217  l=imap(tmp,resultdiff);
218  kill tmp;
219  gresult=gresult,l[2];
220  fresult=simplify(fresult,6);
221  gresult=simplify(gresult,6);
222  fgresult=fresult,gresult;
223  export(fgresult);
224  //fgresult;
225  return(R);
226}
227example
228{
229  "EXAMPLE:"; echo=2;
230  def r= basering;
231  poly f=y2-x3;
232  intvec H1=2,2,2;
233  intvec H2=2,2,2,1;
234  def R1=fgset(f,H1);
235  def R2=fgset(f,H2);
236  setring R1;
237  print(fgresult);
238  setring R2;
239  print(fgresult);
240}
241
242
243
244//////////////////////////////////////////////////////////////////////
245//                                                                  //
246//      PREPARATORY WORK: PLUGGING IN AND DIFFERENTIATING           //
247//                                                                  //
248//////////////////////////////////////////////////////////////////////
249
250
251//////////////////////////////////////////////////////////////////////
252static proc plugincoeffs (poly f,int i)
253"USAGE:   plugincoeffs(f,i); f polynomial, i integer
254 CREATE:  matrix, the t-coefficients obtained by plugging in a
255          polynomial in t of degree i in each variable of f
256 RETURN:  ring, corresponds to f and i, i.e. plugging in a polynomial
257          in t of degree i for each variable of f yields a ring whose
258          variables are the original variables of f, t and the
259          according t-coefficients"
260
261{
262  int startvar=nvars(basering);
263  def r=basering;
264  def R=ringchange(i);              // changes the ring
265  setring R;                        // makes it new basering;
266  ideal I=tpolys(i,startvar);
267  poly g=imap(r,f);                 // maps f to new basering
268  export(g);                        // export it
269  map h=r,I;                        // define map according to our purpose
270  ideal J=h(f);                     // gives f with power series plugged in
271  export(h);
272  matrix resultplug=coeffs(J[1],t); // gives the t-coefficients
273  export resultplug;                // export it i.o. to use it later on
274  return(R);                        // return ring (ring change!)
275}
276
277
278
279//////////////////////////////////////////////////////////////////////
280static proc tpolys (int i,int k)
281"USAGE:   tpolys(i,i); i,k integer
282 RETURN:  ideal, generated by k general polynomials in t of degree i
283          without constant term
284 NOTE:    called from plugincoeffs"
285
286{                                     // has to be called from pluin_coeffs
287  int s,t;                            // loop variables
288  int v;
289  poly sum;
290  ideal I;
291  for(t=1;t<=k;t++)
292   {
293     v=(t-1)*i;
294     for(s=1;s<=i;s++)
295      {
296        sum=sum+var(1+k+v+s)*var(k+1)^s;    // clumsy: working with "var(1)",
297      }                                     // depends on form of basering
298     I[t]=sum;
299     sum=0;
300   }
301  return(I);
302}
303
304
305
306//////////////////////////////////////////////////////////////////////
307static proc formaldiff (poly f,intvec H)
308"USAGE:   formaldiff(f,H); f polynomial, H integer vector
309 CREATE:  list, containing two ideals. Polynomials in t of
310          degree size(H)-1 are plugged into f. H defines, if and how
311          often we differentiate each t-coefficient. We distinguish
312          two different cases, the only difference being that we
313          diffentiate more often in the second case (this is still
314          defined by H). This leads to two systems of equations,
315          each one defining a Zariski-closed set and the second one
316          being contained in the first one
317 RETURN:  ring, corresponds to f and size(H)-1, i.e. plugging in a
318          polynomial in t of degree size(H)-1 for each variable of f
319          yields a ring whose variables are the original variables of f,
320          t and the according t-coefficients"
321{
322  int startvar=nvars(basering);
323  int s,t,v;                          // loop variables
324  int u;
325  int i=size(H)-1;
326  int c=sum(H,1..i);
327  int k=H[i+1];
328  def R=plugincoeffs(f,i);           // plugs the power series in...
329  setring R;                          // changes the ring
330  matrix coe=resultplug;            // gives the t-coeff. after plugging in
331  poly fkv;                           // need this stuff for the following
332  ideal step=diffidealstep(i,startvar);
333  list resultdiff;
334  ideal m,power,diffstep,J,gresultdiff,fresultdiff;
335  for(v=1;v<=k;v++)                   // consider the different t-coeff.
336    {
337     if(c+v>nrows(coe))
338     {
339       fkv=0;
340     }
341     else
342     {
343       fkv=coe[c+v,1];
344     }
345     m=fkv;
346     J=fkv;
347     for(s=1;s<=k-v+1;s++)           // "s<=k-v+1" special for G, ONLY DIFF.!
348       {
349         if(s==k-v+1)                // equations for F!
350         {
351           fresultdiff=fresultdiff,J;
352         }
353         power=step^s;
354         u=size(power);
355         for(t=1;t<=u;t++)
356          {
357            diffstep=contract(power[t],m);     // actual differentiation
358            J=J,diffstep;
359          }
360
361       }
362     gresultdiff=gresultdiff,J;
363   }
364  resultdiff=fresultdiff,gresultdiff;
365  export(resultdiff);                 // exports the result
366  return(R);                           // return the ring
367}
368
369
370
371
372//////////////////////////////////////////////////////////////////////
373//                                                                  //
374//                   CONSTRUCTING THE NEW RING                      //
375//                                                                  //
376//////////////////////////////////////////////////////////////////////
377
378//////////////////////////////////////////////////////////////////////
379static proc ringchange (int i)
380"USAGE:   ringchange(i); i integer
381 RETURN:  ring, extends basering by variables t and
382          #(variables of basering)*i new variables"
383
384{
385  def R=changevar(""+varstr(basering)+",t,"+variables_list(nvars(basering),i)+"");
386  return(R);
387}
388
389
390//////////////////////////////////////////////////////////////////////
391static proc variables_list (int k,int i)
392"USAGE:   variables_list(k,i); k,i integer
393 RETURN:  string of the names of the additional variables
394 NOTE:    called from ringchange, we use this procedure to obtain
395          a convenient shape of the ring created in ringchange"
396{
397  list l;
398  int s,u;                              // loop variables
399  string str;
400  for (u=1;u<=k;u++)
401   {
402     for (s=1;s<=i;s++)
403     {
404       str=""+atoz(u)+"("+string(s)+")"; // creates new variables
405       l[(u-1)*i+s]=str;                // saves them in a list
406     }
407   }
408  string str1=string(l);                // makes the list into a string,
409  return(str1);                         // (needed for ring change)
410}
411
412
413//////////////////////////////////////////////////////////////////////
414static proc atoz (int n)
415"USAGE:   atoz(n); n integer
416 RETURN:  string, the nth letter of the alphabet"
417
418{
419  if(1>n>26)
420   {
421     ERROR("n must range between 1 and 26!");
422   }
423  string s="ring r=0,("+A_Z("a",n)+"),ds;";
424  execute(s);
425  return (string(var(n)));
426}
427
428
429
430//////////////////////////////////////////////////////////////////////
431//                                                                  //
432//                  AUXILIARY PROCEDURES                            //
433//                                                                  //
434//////////////////////////////////////////////////////////////////////
435
436//////////////////////////////////////////////////////////////////////
437static proc diffidealstep (int i, int N)
438"USAGE:   diffidealstep(i,N); i,N integer
439 RETURN:  ideal, generated by variables specified by i,N
440 NOTE:    called from formaldiff, gives the variables by which is
441          differentiated in a certain step"
442{
443  ideal I=var(N+1+i);
444  int j;
445  for(j=2;j<=N;j++)
446   {
447     I=I,var(N+1+j*i);
448   }
449return(I);
450}
451
452
453//////////////////////////////////////////////////////////////////////
454static proc order (poly f)
455"USAGE:   order(f); f polynomial
456 RETURN:  int, the multiplicity of V(f) in 0
457 NOTE:    this order partly directs the differentiation in formaldiff
458          and, together with minpower, gives the size of the
459          integer vector in fgset"
460{
461  poly g=homog(f,var(1));
462  int k=deg(g);
463  int i;
464  for(i=1;i<=k;i++)
465    {
466      if(jet(f,i)!=0)
467        {
468          return(i);
469        }
470    }
471}
472
473
474//////////////////////////////////////////////////////////////////////
475static proc simplifyvar (ideal I)
476"USAGE:   simplifyvar(I); I ideal
477 RETURN:  ideal defining the same zeroset as I: if any generator
478          of I is a power of one single variable, replace it by the
479          variable
480 NOTE:    this procedure is supposed to simplify and clarify the
481          output of the calculations made in fgset, trunc a.o.
482          without using radical"
483{
484  int i,j;
485  int divisornumber=0;
486  int pos;
487  I=simplify(I,6);
488  for(j=1;j<=ncols(I);j++)
489  {
490    if(size(I[j])==1)
491    {
492       for(i=1;i<=nvars(basering);i++)
493       {
494          if(modd(I[j],var(i))==0)
495          {
496             divisornumber++;
497             pos=i;
498          }
499       }
500    }
501    if(divisornumber==1)
502      {
503        I[j]=var(pos);
504      }
505    divisornumber=0;
506  }
507  return(I);
508}
509
510
511//////////////////////////////////////////////////////////////////////
512static proc modd (poly f, poly g)
513"USAGE:   modd(f,g); f,g polynomials
514 RETURN:  poly, f mod g modulo operation in the polynomial ring
515 NOTE:    called from idealsimplify1 where it is used to simplify
516          a generating set of an ideal"
517
518{
519  poly result=f-(f/g)*g;
520  return(result);
521}
522
523
524//////////////////////////////////////////////////////////////////////
525static proc Hnew (intvec H)
526"USAGE:   Hnew(H); H integer vector
527 RETURN:  intvec, the vector needed in the following step of trunc"
528
529{
530  intvec H1=H;
531  int k;
532  int l=size(H);
533  for(k=0;k<=l-2;k++)
534  {
535    if(H[l-k]<H[l-k-1])
536    {
537      H[l-k]=H[l-k]+1;
538      break;
539    }
540  }
541  if(H==H1)
542  {
543    H[l]=H[l]+1;
544  }
545  return(intvec(H));
546}
547
548
549//////////////////////////////////////////////////////////////////////
550static proc simplifymodd (ideal I)
551"USAGE:   simplifymodd(I); I ideal
552 RETURN:  ideal defining the same zeroset as I: replace certain
553          generators of I by the generator modulo the other generators.
554 NOTE:    this procedure is supposed to simplify and clarify the
555          output of the calculations made in fgset, trunc a.o.
556          without using radical"
557{
558  int i,j;
559  I=simplify(I,6);
560  for(j=2;j<=ncols(I);j++)          // reduce with higher element
561  {
562    for(i=1;i<j;i++)
563    {
564      if(I[i]!=0)
565      {
566        I[j]=modd(I[j],I[i]);
567      }
568    }
569  }
570  for(j=ncols(I)-1;j>=1;j--)        // reduce with lower elements
571  {
572    for(i=ncols(I);i>j;i--)
573    {
574      if(I[i]!=0)
575      {
576        I[j]=modd(I[j],I[i]);
577      }
578    }
579  }
580  I=simplify(simplifyvar(I),6);
581  return(I);
582}
583
584
585//////////////////////////////////////////////////////////////////////
586static proc minpower (poly f)
587"USAGE:   minpower(f); f polynomial
588 RETURN:  int, the minimal z>=1 s.th. v^z in J(f) for each variable v
589          of the basering
590 NOTE:    called from trunc, gives; together with i, the size of the
591          integer vectors to be considered in trunc(f,i)"
592{
593 ideal J=jacob(f);
594 int s=ncols(J);
595 int control=0;                // control if conditions for ny are fulfilled
596 int control1=0;
597 int n=nvars(basering);
598 int ny=1;
599 int i,j;
600 while (control==0)           // while var(i)^ny not in J(f)...
601 {
602   for(i=1;i<=n;i++)          // consider all variables
603   {
604    control1=0;
605    for(j=1;j<=s;j++)         // consider all elements of J
606    {
607      if (modd(var(i)^ny,J[j])==0)  // var(i)^ny in J(f)?
608      {
609        control1=1;
610        break;
611      }
612    }
613    if (control1==0)          // increment ny if no var(i)^nt in J(f)
614    {
615      ny++;
616      break;
617    }
618   }
619   if (control1==1)           // if each step was successful...
620  {
621    control=1;
622  }
623 }
624 return(ny);
625}
626
627
628
629static proc equalitytest (ideal I,ideal J)
630"USAGE:   equalitytest(I,J); I,J ideals
631 RETURN:  1, if I=J, 0 else
632 NOTE:    we assume I contained in J already"
633{
634  int s=ncols(J);
635  int i;
636  int p=0;
637  I=std(I);
638  for(i=1;i<=s;i++)
639  {
640    if(reduce(J[i],I)!=0)
641    {
642      p=1;
643    }
644  }
645  if(p==1)
646  {
647    return(0);
648  }
649  else{return(1);}
650}
651
Note: See TracBrowser for help on using the repository browser.