source: git/Singular/LIB/numerAlg.lib @ 810381

spielwiese
Last change on this file since 810381 was 810381, checked in by Hans Schoenemann <hannes@…>, 6 years ago
chg: version/date of libs, p2
  • Property mode set to 100644
File size: 11.8 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="version numerAlg.lib 4.1.1.0 Dec_2017 "; // $Id$
3category="Algebraic Geometry";
4info="
5LIBRARY:  NumerAlg.lib    Numerical Algebraic Algorithm
6OVERVIEW:
7        The library contains procedures to
8        test the inclusion, the equality of two ideals defined by polynomial systems,
9        compute the degree of a pure i-dimensional component of an algebraic variety
10         defined by a polynomial system,
11        compute the local dimension of an algebraic variety defined by a polynomial
12         system at a point computed as an approximate value. The use of the library
13         requires to install Bertini (http://www.nd.edu/~sommese/bertini).
14
15AUTHOR: Shawki AlRashed, rashed@mathematik.uni-kl.de; sh.shawki@yahoo.de
16PROCEDURES:
17
18 Incl(ideal I, ideal J);   test if I containes J
19
20 Equal(ideal I, ideal J);  test if I equals to J
21
22 Degree(ideal I, int i);   computes the degree of a pure i-dimensional
23
24 NumLocalDim(ideal I, p);  numerical local dimension at a point computed as
25                                  an approximate value
26KEYWORDS: bertini; numerDecom_lib
27";
28
29LIB "numerDecom.lib";
30
31///////////////////////////////////////////////////////////////////////////////
32///////////////////////////////////////////////////////////////////////////////
33proc Degree(ideal I,int i)
34"USAGE:  Degree(ideal I,int i); I ideal,  i positive integer
35RETURN:  the degree of the pure i-dimensional component of the algebraic
36          variety defined by I
37EXAMPLE: example Degree; shows an example
38"
39{
40 def S=basering;
41 def W=WitSet(I);
42 setring W;
43 int j;
44 if(size(W(i)[1])>1)
45 {
46  j=size(W(i));
47 }
48 else
49 {
50  j=-1; // no component of dimension i
51 }
52 "The Degree of Component";
53 j;
54 setring S;
55 return (W);
56}
57example
58{ "EXAMPLE:"; echo = 2;
59   ring r=0,(x,y,z),dp;
60   poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
61   poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
62   poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
63   ideal I=f1,f2,f3;
64   def W=Degree(I,1);
65        ==>
66           The Degree of Component
67           3
68   def W=Degree(I,2);
69        ==>
70           The Degree of Component
71           2
72}
73///////////////////////////////////////////////////////////////////////////////
74
75proc Incl(ideal I, ideal J)
76"USAGE:  Incl(ideal I, ideal J); I, J ideals
77RETURN:  t=1 if the algebraic variety defined by I contains the algebraic
78           variety defined by J, otherwise t=0
79EXAMPLE: example Incl; shows an example
80"
81{
82 def S=basering;
83 int n=nvars(basering);
84 int i,j,ii,k,z,zi,dd;
85 if(dim(std(I))==0)
86 {
87  def W=solve(I,"nodisplay");
88  setring W;
89  ideal J=imap(S,J);
90  ideal I=imap(S,I);
91  list w;
92  poly tj;
93  number al,ar,ai,ri,jj;
94  zi=size(SOL);
95  for(j=1;j<=zi;j++)
96  {
97   w=SOL[j];
98   for(k=1;k<=size(J);k++)
99   {
100    tj=J[k];
101    for(ii=1;ii<=n;ii++)
102    {
103     tj=subst(tj,var(ii),w[ii]);
104    }
105    al=leadcoef(tj);
106    ar=repart(al);
107    ai=impart(al);
108    ri=ar^2+ai^2;
109    if(ri>0.000000000000001)
110    {
111     jj=0;
112     k=size(I)+1;
113     j=zi+1;
114    }
115    else
116    {
117     jj=1;
118     ri=0;
119    }
120   }
121  }
122 }
123 else
124 {
125  def W=WitSupSet(I);
126  setring W;
127  ideal J=imap(S,J);
128  ideal I=imap(S,I);
129  list w;
130  number al,ar,ai,ri,jj;
131  poly tj;
132  dd=size(L);
133  for(i=0;i<=dd;i++)
134  {
135   z=size(W(i)[1]);
136   zi=size(W(i));
137   if(z>1)
138   {
139    for(j=1;j<=zi;j++)
140    {
141     w=W(i)[j];
142     for(k=1;k<=size(J);k++)
143     {
144      tj=J[k];
145      for(ii=1;ii<=n;ii++)
146      {
147       tj=subst(tj,var(ii),w[ii]);
148      }
149      al=leadcoef(tj);
150      ar=repart(al);
151      ai=impart(al);
152      ri=ar^2+ai^2;
153      if(ri>0.000000000000001)
154      {
155       jj=-1;
156       k=size(J)+1;
157       j=zi+1;
158       z=0;
159       i=dd+1;
160      }
161      else
162      {
163       jj=1;
164       ri=0;
165      }
166     }
167    }
168   }
169  }
170 }
171 if(ri>0.000000000000001)
172 {
173  jj=0;
174 }
175 else
176 {
177  jj=1;
178 }
179"================================================";
180 "Inclusion:";
181 jj;
182"================================================";
183 export(jj);
184 export(J);
185 export(I);
186   system("sh","rm singular_solutions");
187   system("sh","rm nonsingular_solutions");
188   system("sh","rm real_solutions");
189   system("sh","rm raw_solutions");
190   system("sh","rm raw_data");
191   system("sh","rm output");
192   system("sh","rm midpath_data");
193   system("sh","rm main_data");
194   system("sh","rm input");
195   system("sh","rm failed_paths");
196 setring S;
197 return (W);
198}
199example
200{ "EXAMPLE:"; echo = 2;
201   ring r=0,(x,y,z),dp;
202   poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
203   poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
204   poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
205   ideal I=f1,f2,f3;
206   poly g1=(x2+y2+z2-6)*(x-1);
207   poly g2=(x2+y2+z2-6)*(y-2);
208   poly g3=(x2+y2+z2-6)*(z-3);
209   ideal J=g1,g2,g3;
210   def W=Incl(I,J);
211      ==>
212         Inclusion:
213         0
214 def W=Incl(J,I);
215      ==>
216         Inclusion:
217         1
218}
219///////////////////////////////////////////////////////////////////////////////
220proc Equal(ideal I, ideal J)
221"USAGE:  Equal(ideal I, ideal J); I, J ideals
222RETURN:  t=1 if the algebraic variety defined by I equals to the algebraic
223           variety defined by J, otherwise t=0
224EXAMPLE: example Equal; shows an example
225"
226{
227 def S=basering;
228 int n=nvars(basering);
229 def W1=Incl(J,I);
230 setring W1;
231 number j1=jj;
232 execute("ring q=(real,0),("+varstr(S)+"),dp;");
233 ideal I=imap(W1,I);
234 ideal J=imap(W1,J);
235 execute("ring qq=0,("+varstr(S)+"),dp;");
236 ideal I=imap(S,I);
237 ideal J=imap(S,J);
238 def W2=Incl(I,J);
239 setring W2;
240 number j2=jj;
241 number j;
242 number j1=imap(W1,j1);
243 if(j2==1)
244 {
245  if(j1==1)
246  {
247   j=1/1;
248  }
249  else
250  {
251   j=0/1;
252  }
253 }
254 else
255 {
256  j=0/1;
257 }
258"================================================";
259 "Equality:";
260 j;
261"================================================";
262 setring S;
263 return (W2);
264}
265example
266{ "EXAMPLE:"; echo = 2;
267   ring r=0,(x,y,z),dp;
268   poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
269   poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
270   poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
271   ideal I=f1,f2,f3;
272   poly g1=(x2+y2+z2-6)*(x-1);
273   poly g2=(x2+y2+z2-6)*(y-2);
274   poly g3=(x2+y2+z2-6)*(z-3);
275   ideal J=g1,g2,g3;
276   def W=Equal(I,J);
277        ==>
278           Equality:
279           0
280
281
282  def W=Equal(J,J);
283        ==>
284           Equality:
285           1
286}
287///////////////////////////////////////////////////////////////////////////////
288proc NumLocalDim(ideal J, list w, int e)
289"USAGE:  NumLocalDim(ideal J, list w, int e); J ideal,
290           w list of an approximate value of a point v in the algebraic variety defined by J,
291               e integer
292RETURN: the local dimension of the algebraic variety defined by J at v
293EXAMPLE: example NumLocalDim; shows an example
294"
295{
296 def S=basering;
297 int n=nvars(basering);
298 int sI=size(J);
299 int i,j,jj,t,tt,sz1,sz2,ii,ph,ci,k;
300 poly p,pp;
301 list rw,iw;
302 for(i=1;i<=sI;i++)
303 {
304  p=J[i];
305  for(j=1;j<=n;j++)
306  {
307   w[j]=w[j]+I*0;
308   rw[j]=repart(w[j]);
309   iw[j]=impart(w[j]);
310   p=subst(p,var(j),w[j]);
311  }
312  pp=pp+p;
313 }
314 number u=leadcoef(pp);
315 if((u^2)==0)
316 {
317  execute("ring A=(real,e-1),("+varstr(S)+",I),ds;");
318  ideal II=imap(S,J);
319  list rw=imap(S,rw);
320  list iw=imap(S,iw);
321  poly p(1..n);
322  for(j=1;j<=n;j++)
323  {
324   p(j)=var(j)+rw[j]+I*iw[j];
325  }
326  map f=A,p(1..n);
327  ideal T=f(II);
328  tt=dim(std(T));
329  t=tt-1;
330 }
331 else
332 {
333  int d=dim(std(J));
334  execute("ring R=(complex,e-1,I),("+varstr(S)+"),ds;");
335  list w=imap(S,w);
336  ideal II=imap(S,J);
337  ideal JJ;
338  poly p, p(1..n);
339  for(i=1;i<=sI;i++)
340  {
341   p=II[i];
342   for(j=1;j<=n;j++)
343   {
344    p=subst(p,var(j),w[j]);
345   }
346   JJ[i]=II[i]-p;
347  }
348  for(j=1;j<=n;j++)
349  {
350   p(j)=var(j)+w[j];
351  }
352  map f=R,p(1..n);
353  ideal T=f(JJ);
354  tt=dim(std(T));
355  if(tt==d)
356  {
357   execute("ring A=(complex,e,I),("+varstr(S)+"),dp;");
358   t=tt;
359  }
360  else
361  {
362   execute("ring RR=(real,e-2),("+varstr(S)+",I),dp;");
363   ideal II=imap(S,J);
364   list rw=imap(S,rw);
365   list iw=imap(S,iw);
366   ideal L,LL,H,HH;
367   poly l(1..d),ll(1..d);
368   int c;
369   for(i=1;i<=d;i++)
370   {
371    for(j=1;j<=n;j++)
372    {
373     c=random(1,100);
374     l(i)=l(i)+c*(var(j));
375     ll(i)=ll(i)+c*(var(j)-rw[j]-I*iw[j]);
376    }
377    l(i)=l(i)+random(101,200);
378    L[i]=l(i);
379    LL[i]=ll(i);
380   }
381   poly pi=I^2+1;
382   H=L,II,pi;
383   ideal JJ;
384   poly p, p(1..n);
385   for(i=1;i<=sI;i++)
386   {
387    p=II[i];
388    for(j=1;j<=n;j++)
389    {
390     p=subst(p,var(j),rw[j]+I*iw[j]);
391    }
392    JJ[i]=II[i]-p;
393   }
394   HH=LL,JJ,pi;
395   if(dim(std(H))==0)
396   {
397    def M=solve(H,100,"nodisplay");
398    setring M;
399    sz1=size(SOL);
400    execute("ring RRRQ=(real,e-1),("+varstr(S)+",I),dp;");
401    ideal HH=imap(RR,HH);
402    if(dim(std(HH))==0)
403    {
404     def MM=solve(HH,100,"nodisplay");
405     setring MM;
406     sz2=size(SOL);
407    }
408   }
409   else
410   {
411    sz1=1;
412   }
413   if(sz1==sz2)
414   {
415    execute("ring A=(complex,e,I),("+varstr(S)+"),dp;");
416    t=d;
417   }
418   else
419   {
420    execute("ring RQ=(real,e-1),("+varstr(S)+"),dp;");
421    ideal II=imap(S,J);
422    def RW=WitSet(II);
423    setring RW;
424    list v;
425    list w=imap(S,w);
426    number nr,ni;
427    if(tt<0)
428    {
429     tt=0;
430    }
431    for(ii=tt;ii<=d;ii++)
432    {
433     list W(ii)=imap(RW,W(ii));
434     if(size(W(ii)[1])>1)
435     {
436      if(ii==0)
437      {
438       for(i=1;i<=size(W(0));i++)
439       {
440        v=W(ii)[i];
441        nr=0;
442        ni=0;
443        for(j=1;j<=n;j++)
444        {
445         nr=nr+(repart(v[j])-repart(w[j]))^2;
446         ni=ni+(impart(v[j])-impart(w[j]))^2;
447        }
448        if((ni+nr)<1/10^(2*e-3))
449        {
450         execute("ring A=(complex,e,I),("+varstr(S)+"),dp;");
451         list W(ii)=imap(RW,W(ii));
452         t=0;
453         i=size(W(ii))+1;
454         ii=d+1;
455        }
456       }
457      }
458      else
459      {
460       def SS=Singular2bertini(W(ii));
461       execute("ring D=(complex,e,I),("+varstr(S)+",s,gamma),dp;");
462       string nonsin;
463       ideal H,L;
464       ideal J=imap(RW,N(0));
465       ideal LL=imap(RW,L);
466       list w=imap(S,w);
467       poly p;
468       for(j=1;j<=ii;j++)
469       {
470        p=0;
471        for(jj=1;jj<=n;jj++)
472        {
473         p=p+random(1,100)*(var(jj)-w[jj]);
474        }
475        L[j]=p;
476       }
477       for(jj=1;jj<=size(J);jj++)
478       {
479        H[jj]=s*gamma*J[jj]+(1-s)*J[jj];
480       }
481       for(jj=1;jj<=ii;jj++)
482       {
483        H[size(J)+jj]=s*gamma*LL[jj]+(1-s)*L[jj];
484       }
485       string sv=varstr(S);
486       def Q(ii)=UseBertini(H,sv);
487       system("sh","rm start");
488       nonsin=read("nonsingular_solutions");
489       if(size(nonsin)>=52)
490       {
491        def T(ii)=bertini2Singular("nonsingular_solutions",nvars(basering)-2);
492        setring T(ii);
493        list C=re;
494        ci=size(C);
495        number tr;
496        list w=imap(S,w);
497        for(jj=1;jj<=ci;jj++)
498        {
499         tr=0;
500         for(k=1;k<=n;k++)
501         {
502          tr=tr+(repart(w[k])-repart(C[jj][k]))^2+(impart(w[k])-impart(C[jj][k]))^2;
503         }
504         if(tr<=1/10^(2*e-3))
505         {
506          execute("ring A=(complex,e,I),("+varstr(S)+"),dp;");
507          t=ii;
508          ii=d+1;
509          jj=ci+1;
510         }
511        }
512       }
513      }
514     }
515    }
516    system("sh","rm singular_solutions");
517    system("sh","rm nonsingular_solutions");
518    system("sh","rm real_solutions");
519    system("sh","rm raw_solutions");
520    system("sh","rm raw_data");
521    system("sh","rm output");
522    system("sh","rm midpath_data");
523    system("sh","rm main_data");
524    system("sh","rm input");
525    system("sh","rm failed_paths");
526   }
527  }
528 }
529 "=============================================";
530 "The Local Dimension:";
531 t;
532 setring S;
533 return(A);
534}
535example
536{ "EXAMPLE:"; echo = 2;
537   int e=14;
538   ring r=(complex,e,I),(x,y,z),dp;
539   poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
540   poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
541   poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
542   ideal J=f1,f2,f3;
543   list p0=0.99999999999999+I*0.00000000000001,2,3+I*0.00000000000001;
544   list p2=1,0.99999999999998,2;
545   list p1=5+I,4.999999999999998+I,5+I;
546   def D=NumLocalDim(J,p0,e);
547             ==>
548               The Local Dimension:
549                0
550   def D=NumLocalDim(J,p1,e);
551             ==>
552               The Local Dimension:
553                1
554   def D=NumLocalDim(J,p2,e);
555             ==>
556               The Local Dimension:
557                2
558}
559
560///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.