source: git/Singular/LIB/numerDecom.lib @ ba64bb

spielwiese
Last change on this file since ba64bb was ba64bb, checked in by Andreas Steenpass <steenpass@…>, 5 years ago
do not load ring.lib for using create_ring()
  • Property mode set to 100644
File size: 45.5 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="version numerDecom.lib 4.1.2.0 Feb_2019 "; // $Id$
3category="Algebraic Geometry";
4info="
5LIBRARY:  NumDecom.lib    Numerical Decomposition of Ideals
6OVERVIEW:
7     The library contains procedures to compute
8      numerical irreducible decomposition, and
9      numerical primary decomposition of an algebraic variety defined by a
10      polynomial system. The use of the library requires to install Bertini
11           (@uref{http://www.nd.edu/~sommese/bertini}).
12AUTHOR: Shawki AlRashed, rashed@mathematik.uni-kl.de; sh.shawki@yahoo.de
13
14PROCEDURES:
15
16 re2squ(ideal I);               reduction to square system
17
18 UseBertini(ideal H,string sv); use Bertini to compute the solutions of the homotopy function
19
20 Singular2bertini(list L); adopt the list to be a read file in Bertini as a start solution set
21
22 bertini2Singular(string snp, int q); adopt the file of solutions of the homotopy function to be a list in SINGULAR
23
24 ReJunkUseHomo(ideal I, ideal L, list W, list w); remove junk points using the homotopy function
25
26 JuReTopDim(ideal J,list w,int tt, int d); remove junk points that are on top-dimensional component
27
28 JuReZeroDim(ideal J,list w, int d); remove junk points from 0-dimensional component
29
30 WitSupSet(ideal I);                 witness point super set
31
32 WitSet(ideal I);                    witness point set
33
34 NumIrrDecom(ideal I);               numerical irreducible decomposition
35
36 defl(ideal I, int d);               deflation of ideal I
37
38 NumPrimDecom(ideal I, int d);       numerical primary decomposition
39
40KEYWORDS: numerical irreducible decomposition; decomposition, numerical; primary decomposition, numerical;bertini; numerAlg_lib
41";
42
43
44LIB "solve.lib";
45LIB "matrix.lib";
46
47///////////////////////////////////////////////////////////////////////////////
48///////////////////////////////////////////////////////////////////////////////
49proc re2squ(ideal I)
50"USAGE:   re2squ(ideal I);I ideal
51RETURN:   ideal J defined by the polynomial system of the same number of polynomials and unknowns
52EXAMPLE: example re2squ;shows an example
53"
54{
55 def S=basering;
56 int n=nvars(basering);
57 ideal J;
58 poly p;
59 int N=size(I);
60 int i,j;
61 if(n==N)
62 {
63  J=I;
64 }
65 else
66 {
67  if(N<n)
68  {
69   for(i=1;i<=n;i++)
70   {
71    if(i<=N)
72    {
73     J[i]=I[i];
74    }
75    else
76    {
77     J[i]=0;
78    }
79   }
80  }
81  else
82  {
83   for(i=1;i<=n;i++)
84   {
85    p=0;
86    for(j=N-n;j<=N;j++)
87    {
88     p=p+random(1,101)*I[j];
89    }
90    J[i]=I[i]+p;
91   }
92  }
93 }
94 export(J);
95 setring S;
96 return(S);
97 }
98example
99{ "EXAMPLE:";echo = 2;
100   ring r=0,(x,y,z),dp;
101   ideal I= x3+y4,z4+yx,xz+3x,x2y+z;
102   def D=re2squ(I);
103   setring D;
104   J;
105}
106///////////////////////////////////////////////////////////////////////////////
107
108proc Singular2bertini(list L)
109"USAGE:   Singular2bertini(list L);L a list
110RETURN:  text file called start
111NOTE:    adopting the list L to be as a start solution of the homptopy function in Bertini
112EXAMPLE: Singular2bertini;shows an example
113"
114{
115 write("start",string(size(L)));
116 int i,j;
117 number a,b;
118 string s;
119 list LLL;
120 for(i=1;i<=size(L);i++)
121 {
122  LLL=L[i];
123  for(j=1;j<=size(LLL);j++)
124  {
125   a=repart(LLL[j]);
126   b=impart(LLL[j]);
127   s=string(a)+" "+string(b)+";";
128  write("start",s);
129  }
130 }
131 return(0);
132}
133example
134{ "EXAMPLE:";echo = 2;
135   ring r=(complex,16,I),(x,y,z),dp;
136   list L=list(1,2,3),list(4,5,6+I*2);
137   def D=Singular2bertini(L);
138}
139///////////////////////////////////////////////////////////////////////////////
140proc UseBertini(ideal H,string sv)
141"USAGE:   UseBertini(ideal H,string sv);
142          H ideal, sv string of the variable of ring
143RETURN:   text file called input used in Bertini to compute the solution of
144          the homotopy function H that existed in file input.text
145NOTE:     Need to define a start solution of H
146EXAMPLE:  example Use;shows an example
147"
148{
149 int ii,j,k, ph;
150 ph=size(H);
151 string sff,sf;
152 link l=":w ./input";
153 write(l,"");
154 write(l,"CONFIG");
155 write(l,"");
156 write(l,"USERHOMOTOPY: 1;");
157 write(l,"");
158 write(l,"END;");
159 write(l,"");
160 write(l,"INPUT");
161 write(l,"");
162 for( ii=1;ii<=size(sv);ii++)
163 {
164  if((sv[ii]=="(")||(sv[ii]==")"))
165  {
166   sv=sv[1,ii-1]+sv[ii+1,size(sv)];
167  }
168 }
169 write(l,"variable "+sv+";");
170  sff="function";
171 if(ph!=1)
172 {
173  for( ii=1;ii<=ph-1;ii++)
174  {
175   sff=sff+" f"+string(ii)+",";
176  }
177  sff=sff+"f"+string(ph)+";";
178 }
179 else
180 {
181  sff=sff+" f"+string(1)+",";
182 }
183 write(l,sff);
184 write(l,"pathvariable t;");
185 write(l,"parameter s;");
186 write(l,"constant gamma;");
187 write(l,"");
188 write(l,"gamma = 0.8 + 1.1*I;");
189 write(l,"");
190 write(l,"s=t;");
191 write(l,"");
192 short=0;
193 for( ii=1;ii<=ph;ii++)
194 {
195   sf=string(H[ii]);
196   k=find(sf,newline);
197  for( j=1;j<=size(sf);j++)
198  {
199   if(sf[j]=="(")
200   {
201    if(sf[j+2]==")")
202    {
203     sf[j]=" ";
204     sf=sf[1,j-1]+sf[j+1,size(sf)];
205     sf[j+1]=" ";
206     sf=sf[1,j]+sf[j+2,size(sf)];
207    }
208   }
209  }
210  write(l,"f"+string(ii)+"="+sf+";");
211 }
212 write(l,"END;");
213 system(("sh","bertini<./input"));
214 return(0);
215}
216example
217{ "EXAMPLE:";echo = 2;
218   ring r=0,(x,y,z),dp;
219   ideal I= x3+y4,z4+yx,xz+3x,x2y+z;
220   string sv=varstr(basering);
221   def A=UseBertini(I,sv);
222}
223///////////////////////////////////////////////////////////////////////////////
224proc bertini2Singular(string snp, int q)
225"USAGE:   bertini2Singular(string snp, int q);
226         snp string, q=nvars(basering) integer
227RETURN:  re list of the solutions of the homotopy function computed by Bertini
228EXAMPLE: example bertini2Singular;shows an example
229"
230{
231 def S=basering;
232 int nn=nvars(basering);
233 int n=q;
234 ring R = create_ring("(complex,18,I)", "("+varstr(S)+")", "dp");
235 number r1,r2;
236 list re,ru;
237 string sss=read(snp);
238 sss=sss+"";
239 int i,j,k,m,p;
240 string ss;
241 ss=sss[1];
242 i=2;
243 while(sss[i]!=" ")
244 {
245  ss=ss+sss[i];
246  i++;
247 }
248 execute("m="+ss+";");
249 for(i=1;i<=size(sss);i++)
250 {
251  if(sss[i]=="e")
252  {
253   if(!((sss[i+1]=="+")||(sss[i+1]=="-")))
254   {
255    ss=sss[i+1,size(sss)];
256    sss=sss[1,i];
257    sss=sss+"+"+ss;
258   }
259  }
260 }
261 j=1;
262 j=find(sss,newline,j)+1;
263 while(sss[j]==newline){j++;}
264 for(q=1;q<=m;q++)
265 {
266  for(p=1;p<=n;p++)
267  {
268   k=find(sss,newline,j);
269   ss=sss[j,k-j];
270   i=find (ss," ");
271   execute("r1="+ss[1,i-1]+";");
272   execute("r2="+ss[i+1,size(ss)-i+1]+";");
273   ru[p]=r1+I*r2;
274   j=k+1;
275  }
276  j=j+1;
277  re[size(re)+1]=ru;
278 }
279 export(re);
280 setring S;
281 return(R);
282}
283example
284{ "EXAMPLE:";echo = 2;
285   ring r = 0,(a,b,c),ds;
286   int q=nvars(basering);
287   def T=bertini2Singular("nonsingular_solutions",q);
288   re;
289}
290///////////////////////////////////////////////////////////////////////////////
291proc   WitSupSet(ideal I)
292"USAGE:  WitSupSet(ideal I);I ideal
293RETURN:  list of Witness point Super Sets W(i) for i=1,...,dim(V(I)),
294          L list of generic linear polynomials and N(0) list of a polynomial system of the same number of
295          polynomials and unknowns. // if W(i) = x, then V(I) has no component of dimension i
296EXAMPLE: example WitSupSet;shows an example
297"
298{
299 def S=basering;
300 int n=nvars(basering);
301 ideal II=I;
302 int dd=dim(std(I));
303 if(n==1)
304 {
305   ERROR("n=1");
306 }
307 else
308 {
309  if(dd==0)
310  {
311   ring R = create_ring(0, "("+varstr(S)+")", "dp");
312   int i,j;
313   ideal I=imap(S,I);
314   list V(dd),W(dd);
315   def T(dd+1)=solve(I,"nodisplay");
316   setring T(dd+1);
317   W(dd)=SOL;
318   ideal N(dd)=imap(S,I);
319   export(N(dd));
320   ideal LL;
321   export(LL);
322   int c(0);
323   c(0)=0;
324   export(c(0));
325   export(dd);
326   list w(1..size(W(dd)));
327   for(i=1;i<=size(W(dd));i++)
328   {
329    w(i)=W(dd)[i];
330    export(w(i));
331   }
332  "===========================================";
333  "===========================================";
334   "Dimension";
335    dd;
336   "Number of Components";
337   size(W(dd));
338   setring S;
339   return(T(dd+1));
340  }
341  else
342  {
343   matrix MJJ=jacob(I);
344   int rn=rank(MJJ);
345   I=imap(S,II);
346   def rs=re2squ(I);
347   setring rs;
348   I=J;
349   if((n-rn)!=dd)
350   {
351    execute("ring R=0,("+varstr(S)+",z(1..dd)),dp;");
352    ideal I=imap(rs,I);
353    ideal H(0..n),L,LL,L(1..dd),LL(1..dd),h(1..dd),N(0..dd);
354    poly p,p(0..n),e;
355   int i,j,k,kk,q,qq,t,m,d,jj,rii,c(0),ii;
356   for(i=1;i<=dd;i++)
357   {
358    p=0;
359    for(j=1;j<=n;j++)
360    {
361     p=p+random(1,2*n+7)*var(j);
362    }
363    if(i<dd)
364    {
365     LL[i]=random(2*n+7,4*n+1)+p;
366    }
367    else
368    {
369     c(0)=random(4*n+1,5*n+13);
370     LL[i]=c(0)+p;
371    }
372   }
373   export(c(0));
374   p(0)=0;
375   for(t=1;t<=n;t++)
376   {
377    for(j=1;j<=dd;j++)
378    {
379     p(j)=p(j-1)+random(1,2*n+10)*var(n+j);
380     h(j)[t]=I[t]+p(j);
381    }
382   }
383   for(q=1;q<=dd;q++)
384   {
385    for(i=1;i<=q;i++)
386    {
387     L(q)[i]=LL[i]+var(n+i);
388    }
389   }
390   for(i=1;i<=dd;i++)
391   {
392    N(i)=h(i),L(i);
393   }
394   for(i=1;i<=n;i++)
395   {
396    N(0)[i]=I[i];
397   }
398   ideal JJ=N(0);
399   if(dim(std(N(dd)))!=0)
400   {
401    "Try Again";
402   }
403   else
404   {
405    def T=solve(N(dd),100,"nodisplay");
406    setring T;
407    ring T(dd+1) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
408    list M,Y;
409    list W(dd),V(dd);
410    list SOL=imap(T,SOL);
411    Y=SOL;
412    number rp,ip,rip;
413    for( i=1;i<=size(SOL);i++)
414    {
415     M=Y[i];
416     for(j=dd;j>=1;j--)
417     {
418      rp=repart(M[n+j]);
419      ip=impart(M[n+j]);
420      rip=rp^2 + ip^2;
421      if(rip<0.0000000000000001)
422      {
423       M=delete(M,n+j);
424       Y[i]=M;
425      }
426     }
427    }
428    k=1;
429    kk=1;
430    for( i=1;i<=size(Y);i++)
431    {
432     if(size(Y[i])==n)
433     {
434      W(dd)[k]=Y[i];
435      k=k+1;
436     }
437     else
438     {
439      V(dd)[kk]=Y[i];
440      kk=kk+1;
441     }
442    }
443    ideal JJ=imap(S,II);
444    k=1;
445    number al,ar,ai,ri;
446    for(j=1;j<=size(W(dd));j++)
447    {
448     ri=0;
449     al=0;
450     ai=0;
451     ar=0;
452     for(ii=1;ii<=size(JJ);ii++)
453     {
454      for(i=1;i<=n;i++)
455      {
456       JJ[ii]=subst(JJ[ii],var(i),W(dd)[j][i]);
457      }
458      al=leadcoef(JJ[ii]);
459      ar=repart(al);
460      ai=impart(al);
461      ri=ar^2+ai^2+ri;
462     }
463     if(ri<=0.000000000000000001)
464     {
465      W(dd)[k]=W(dd)[j];
466      k=k+1;
467     }
468    }
469    ideal L(dd)=imap(R,L(dd));
470    export(L(dd));
471    export(W(dd));
472    export(V(dd));
473    string sff,sf,sv;
474    int nv(dd)=size(V(dd));
475    int nv(0..dd-1);
476    if(size(W(dd))<size(Y))
477    {
478     def SB(dd)=Singular2bertini(V(dd));
479    }
480    for( q=dd;q>=n-rn+1;q--)
481    {
482     if(nv(q)!=0)
483     {
484      int w(q-1)=0;
485      execute("ring D(q)=(0,s,gamma),("+varstr(S)+",z(1..q)),dp;");
486      string nonsin(q),stnonsin(q);
487      ideal H(1..q);
488      ideal N(q)=imap(R,N(q));
489      ideal N(q-1)=imap(R,N(q-1));
490      for(j=1;j<=n+q-1;j++)
491      {
492       H(q)[j]=s*gamma*N(q)[j]+(1-s)*N(q-1)[j];
493      }
494      H(q)[n+q]=s*gamma*N(q)[n+q]+(1-s)*var(n+q);
495      ideal H=H(q);
496      export(H(q));
497      string sv(q)=varstr(basering);
498      sv=sv(q);
499      def Q(q)=UseBertini(H,sv);
500      system("sh","rm start");
501      nonsin(q)=read("nonsingular_solutions");
502      if(size(nonsin(q))>=52)
503      {
504       def T(q)=bertini2Singular("nonsingular_solutions",nvars(basering));
505       setring T(q);
506       list C=re;
507       list B,X,A,G;
508       for(i=1;i<=size(C);i++)
509       {
510        B=re[i];
511        B=delete(B,n+q);
512        C[i]=B;
513       }
514       X=C;
515       if(q>=2)
516       {
517        for(j=q-1;j>=1;j--)
518        {
519         for(i=1;i<=size(X);i++)
520         {
521          A[i]=X[i];
522          G=A[i];
523          G=delete(G,n+j);
524          A[i]=G;
525         }
526         X=A;
527        }
528       }
529       else
530       {
531        X=C;
532       }
533       list W(q-1),V(q-1);
534       ideal JJ=imap(S,II);
535       k=1;
536       poly tj;
537       number al,ar,ai,ri;
538       for(j=1;j<=size(C);j++)
539       {
540        ri=0;
541        al=0;
542        ai=0;
543        ar=0;
544        for(i=1;i<=size(JJ);i++)
545        {
546         tj=JJ[i];
547         for(i=1;i<=n;i++)
548         {
549          tj=subst(tj,var(i),X[j][i]);
550         }
551         al=leadcoef(tj);
552         ar=repart(al);
553         ai=impart(al);
554         ri=ar^2+ai^2+ri;
555        }
556        if(ri<=0.000000000000000001)
557        {
558          W(q-1)[k]=X[j];
559          k=k+1;
560        }
561        else
562        {
563         nv(q-1)=nv(q-1)+1;
564         V(q-1)[nv(q-1)]=C[j];
565        }
566       }
567       if(nv(q-1)==size(C))
568       {
569        list W(q-1)=var(1);
570       }
571       if(q>=2)
572       {
573        if(nv(q-1)!=0)
574        {
575         def SB(qq-1)=Singular2bertini(V(q-1));
576        }
577        else
578        {
579         for(qq=q-1;qq>=1;qq--)
580         {
581          execute("ring T(qq)=(complex,16,I),("+varstr(S)+",z(1..qq)),dp;");
582          list W(qq-1)=var(1);
583         }
584         q=1;
585        }
586       }
587      }
588      else
589      {
590       for(qq=q;qq>=1;qq--)
591       {
592        int w(qq-1);
593        execute("ring T(qq)=(complex,16,I),("+varstr(S)+",z(1..qq)),dp;");
594        list W(qq-1)=var(1);
595       }
596      }
597     }
598     else
599     {
600      for(qq=q;qq>=1;qq--)
601      {
602       execute("ring T(qq)=(complex,16,I),("+varstr(S)+",z(1..qq)),dp;");
603       list W(qq-1)=var(1);
604      }
605     }
606    }
607    ring D = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
608    for(i=0;i<=dd;i++)
609    {
610     list W(i)=imap(T(i+1),W(i));
611     export(W(i));
612    }
613     ideal L=imap(R,LL);
614     export(L);
615     ideal N(0)=imap(R,N(0));
616     export(N(0));
617     setring S;
618     return(D);
619    }
620   }
621   else
622   {
623    ring R = create_ring(0, "("+varstr(S)+")", "dp");
624    int i,j,c(0);
625    poly p;
626    ideal LL;
627    for(i=1;i<=dd;i++)
628    {
629     p=0;
630     for(j=1;j<=n;j++)
631     {
632      p=p+random(1,100)*var(j);
633     }
634     if(i<dd)
635     {
636      LL[i]=random(101,200)+p;
637     }
638     else
639     {
640      c(0)=random(201,300);
641      LL[i]=c(0)+p;
642     }
643    }
644    ideal I=imap(S,I);
645    ideal N(dd)=I,LL;
646    def T=solve(N(dd),100,"nodisplay");
647    setring T;
648    list W(0..dd);
649    W(dd)=SOL;
650    export(W(dd));
651    for(i=0;i<=dd-1;i++)
652    {
653     W(i)=var(1);
654     export(W(i));
655    }
656    ideal L=imap(R,LL);
657    export(L);
658    ideal N(0)=imap(S,I);
659    export(N(0));
660    setring S;
661    return(T);
662   }
663  }
664 }
665}
666example
667{ "EXAMPLE:";echo = 2;
668    ring r=0,(x,y,z),dp;
669    poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
670    poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
671    poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
672    ideal I=f1,f2,f3;
673    def W=WitSupSet(I);
674    setring W;
675    W(2);
676        // witness point super set of a pure 2-dimensional component of V(I)
677    W(1);
678        // witness point super set of a pure 1-dimensional component of V(I)
679    W(0);
680        // witness point super set of a pure 0-dimensional component of V(I)
681    L;
682        // list of generic linear polynomials
683}
684///////////////////////////////////////////////////////////////////////////////
685proc ReJunkUseHomo(ideal I, ideal L, list W, list w)
686"USAGE:   ReJunkUseHomo(ideal I, ideal L, list W, list w);
687          I ideal, L list of generic linear polynomials {l_1,...,l_i},
688          W list of a subset of the solution set of the generic slicing V(L) with V(J),
689          w list of a point in V(J)
690RETURN:  t=1 if w on an i-dimensional component of V(I),
691         otherwise t=0. Where i=size(L)
692EXAMPLE: example ReJunkUseHomo;shows an example
693"
694{
695 def S=basering;
696 int n=nvars(basering);
697 int ii,i,in,j,jjj,jj,k,zi,a,kk,kkk;
698 string sf,sff,sv;
699 i=size(W);
700 in=size(w);
701 ideal LL;
702 jjj=size(L);
703 poly pp;
704 for(jj=1;jj<=jjj;jj++)
705 {
706  for(ii=1;ii<=in;ii++)
707  {
708   pp=random(1,3*n+1)*(var(ii)-w[ii])+pp;
709  }
710  LL[jj]=pp;
711 }
712 export(LL);
713 ring R = create_ring("(complex,16,I)", "("+varstr(S)+",gamma,s)", "dp");
714 ideal L=imap(S,L);
715 ideal LL=imap(S,LL);
716 ideal I=imap(S,I);
717 list w=imap(S,w);
718 zi=size(I);
719 ideal H;
720 for(a=1;a<=zi;a++)
721 {
722  H[a]=s*gamma*I[a]+(1-s)*I[a];
723 }
724 for(kk=1;kk<=jjj;kk++)
725 {
726  H[kk+zi]=s*gamma*L[kk]+(1-s)*LL[kk];
727 }
728 list W=imap(S,W);
729 def SB1=Singular2bertini(W);
730 sv=varstr(S);
731 def Q=UseBertini(H,sv);
732 system("sh","rm start");
733 string nonsin=read("nonsingular_solutions");
734 if(size(nonsin)>=52)
735 {
736  def TT=bertini2Singular("nonsingular_solutions",nvars(basering)-2);
737  setring TT;
738  list w=imap(S,w);
739  list C=re;
740  list ww,v;
741  number rp,ip,rp(1..size(w)),ip(1..size(w)),irp,t;
742  for(k=1;k<=size(C);k++)
743  {
744   ww=re[k];
745   for(jj=1;jj<=size(w);jj++)
746   {
747    rp(jj)=(repart(ww[jj])-repart(w[jj]))^2;
748    ip(jj)=(impart(ww[jj])-impart(w[jj]))^2;
749    rp=rp+rp(jj);
750    ip=ip+ip(jj);
751   }
752   irp=ip+rp;
753   if(irp<=0.000000000000000000000001)
754   {
755    t=1.0;
756   }
757   else
758   {
759    t=0.0;
760   }
761  }
762 }
763 else
764 {
765 ring TT = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
766 list w=imap(S,w);
767 number t=1.0;
768 }
769 export(t);
770 setring S;
771 return(TT);
772}
773example
774{ "EXAMPLE:";echo = 2;
775    ring r=(complex,16,I),(x,y,z),dp;
776    poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
777    poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
778    poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
779    ideal J=f1,f2,f3;
780    poly l1=15x+16y+6z+17;
781    poly l2=2x+14y+4z+18;
782    ideal L=l1,l2;
783    list W1=list(0.5372775295412116,-0.7105339291010922,-2.2817700129167831+I*0),list(0.09201175741935605,-1.7791717821935455,1.6810953589677311);
784    list w=list(2,2,-131666666/10000000);
785    def D=ReJunkUseHomo(J,L,W1,w);
786    setring D;
787    t;
788}
789///////////////////////////////////////////////////////////////////////////////
790proc JuReTopDim(ideal J,list w,int tt, int d);
791"USAGE:  JuReTopDim(ideal J,list w,int tt, int d);J ideal, w list of a point in V(J),
792         tt the degree of d-dimensional component of V(J), d dimension of V(J)
793RETURN:  t=1 if w on a d-dimensional component of V(I), otherwise t=0.
794EXAMPLE: example JuReTopDim;shows an example
795"
796{
797 def S=basering;
798 int n=nvars(basering);
799 int i,j,k;
800 list iw,rw;
801 for(k=1;k<=n;k++)
802 {
803  rw[k]=repart(w[k]);
804  iw[k]=impart(w[k]);
805 }
806 ring R = create_ring(list("real"), "("+varstr(S)+",I)", "dp");
807 list iw=imap(S,iw);
808 list rw=imap(S,rw);
809 ideal J=imap(S,J);
810 ring RR = create_ring(0, "("+varstr(S)+",I)", "dp");
811 ideal J=imap(R,J);
812 list iw=imap(R,iw);
813 list rw=imap(R,rw);
814 ideal L;
815 poly p;
816 for(i=1;i<=d;i++)
817 {
818  p=0;
819  for(j=1;j<=n;j++)
820  {
821   p=p+random(1,100)*(var(j)-rw[j]-I*iw[j]);
822  }
823  L[i]=p;
824 }
825 ideal JJ;
826 for(i=1;i<=size(J);i++)
827 {
828  p=J[i];
829  for(j=1;j<=n;j++)
830  {
831   p=subst(p,var(j),rw[j]+I*iw[j]);
832  }
833  JJ[i]=p;
834 }
835 poly pp;
836 pp=I^2 +1;
837 ideal T=L,J,pp;
838 int di=dim(std(T));
839 if(di==0)
840 {
841  def T(d)=solve(T,10,"nodisplay");
842  setring T(d);
843  number t,ie,re,rt;
844  int zi=size(SOL);
845  list iw=imap(S,iw);
846  list rw=imap(S,rw);
847  if(zi==2*tt)
848  {
849   t=1.0/1;
850  }
851  else
852  {
853   t=0.0/1;
854  }
855 }
856 else
857 {
858  ring T(d) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
859  "Try Again";
860   -----
861 }
862 export(t);
863 setring S;
864 return(T(d));
865}
866example
867{ "EXAMPLE:";echo = 2;
868    ring r=(complex,16,I),(x,y,z),dp;
869    poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
870    poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
871    poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
872    ideal J=f1,f2,f3;
873    list w=list(0.5372775295412116,-0.7105339291010922,-2.2817700129167831);
874    def D=JuReTopDim(J,w,2,2);
875    setring D;
876    t;
877}
878///////////////////////////////////////////////////////////////////////////////
879proc JuReZeroDim(ideal J,list w, int d);
880"USAGE:  JuReZeroDim(ideal J,list w, int d);J ideal,
881         w list of a point in V(J), d dimension of V(J)
882RETURN:  t=1 if w on a positive-dimensional component of V(I),
883         i.e w is not isolated point in V(J)
884EXAMPLE: example JuReZeroDim;shows an example
885"
886{
887 def S=basering;
888 int n=nvars(basering);
889 int i,j,k;
890 list iw,rw;
891 for(k=1;k<=n;k++)
892 {
893  rw[k]=repart(w[k]);
894  iw[k]=impart(w[k]);
895 }
896 ring R = create_ring(list("real"), "("+varstr(S)+",I)", "dp");
897 list iw=imap(S,iw);
898 ideal J=imap(S,J);
899 list rw=imap(S,rw);
900 ring RR = create_ring(0, "("+varstr(S)+",I)", "dp");
901 list iw=imap(R,iw);
902 ideal J=imap(R,J);
903 list rw=imap(R,rw);
904 ideal LL;
905 poly p;
906 for(i=1;i<=d;i++)
907 {
908  p=0;
909  for(j=1;j<=n;j++)
910  {
911   p=p+random(1,100)*(var(j)-rw[j]-I*iw[j]);
912  }
913  LL[i]=p;
914 }
915 poly pp;
916 pp=I^2 +1;
917 ideal TT=LL,J,pp;
918 def TT(d)=solve(TT,16,"nodisplay");
919 setring TT(d);
920 int zii=size(SOL);
921 ring RR1 = create_ring(0, "("+varstr(S)+",I)", "dp");
922 list iw=imap(R,iw);
923 ideal J=imap(R,J);
924 list rw=imap(R,rw);
925 ideal L;
926 poly p;
927 for(i=1;i<=d;i++)
928 {
929  p=0;
930  for(j=1;j<=n;j++)
931  {
932   p=p+random(1,100)*(var(j)-rw[j]-I*iw[j]-1/100000000000000);
933  }
934  L[i]=p;
935 }
936 poly pp;
937 pp=I^2 +1;
938 ideal T=L,J,pp;
939 int di=dim(std(T));
940 if(di==0)
941 {
942  def T(d)=solve(T,16,"nodisplay");
943  setring T(d);
944  number t;
945  int zi=size(SOL);
946  list iw=imap(S,iw);
947  list rw=imap(S,rw);
948  if(zi==zii)
949  {
950   t=1.0/1;
951  }
952  else
953  {
954   t=0.0/1;
955  }
956 }
957 else
958 {
959  ring T(d) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
960  "Try Again";
961   -----
962 }
963 export(t);
964 setring S;
965 return(T(d));
966}
967example
968{ "EXAMPLE:";echo = 2;
969    ring r=(complex,16,I),(x,y,z),dp;
970    poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
971    poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
972    poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
973    ideal J=f1,f2,f3;
974    list w1=list(0.5372775295412116,-0.7105339291010922,-2.2817700129167831);
975    def D1=JuReZeroDim(J,w1,2);
976    setring D1;
977    t;
978}
979///////////////////////////////////////////////////////////////////////////////
980proc    WitSet(ideal I)
981"USAGE:   WitSet(ideal I); I ideal
982RETURN:  lists W(0..d) of witness point sets of i-dimensional components
983          of V(J) for i=0,...d respectively, where d the dimension of V(J),
984          L list of generic linear polynomials
985NOTE:    if W(i)=x, then V(J) has no component of dimension i
986EXAMPLE: example WitSet;shows an example
987"
988{
989 def S=basering;
990 int n=nvars(basering);
991 int ii,i,j,b,bb,k,kk,dt;
992 def TJ(0)=WitSupSet(I);
993 setring TJ(0);
994 ideal LL=L;
995 int d=size(LL);
996 if(d==0)
997 {
998  setring S;
999  return(TJ(0));
1000 }
1001 else
1002 {
1003  for( i=0;i<=d;i++)
1004  {
1005   list Ww(i)=W(i);
1006   int z(i)=size(W(i));
1007   export(Ww(i));
1008  }
1009  for(i=d-1;i>=0;i--)
1010  {
1011   list W(i)=imap(TJ(0),Ww(i));
1012   if(size(W(i)[1])>1)
1013   {
1014    for(j=1;j<=z(i);j++)
1015    {
1016     ring Rr(j+i) = create_ring("(complex,106,I)", "("+varstr(S)+")", "ds");
1017     list W(i)=imap(TJ(0),Ww(i));
1018     list w=W(i)[j];
1019     ideal J=imap(TJ(0),N(0));
1020     ideal J(j),K(j);
1021     for(k=1;k<=size(J);k++)
1022     {
1023      J(j)[k]=J[k];
1024      for(kk=1;kk<=n;kk++)
1025      {
1026       J(j)[k]=subst(J(j)[k],var(kk),w[kk]);
1027      }
1028      K(j)[k]=J[k]-J(j)[k];
1029     }
1030     poly p(1..n);
1031     for(k=1;k<=n;k++)
1032     {
1033      p(k)=var(k)+w[k];
1034     }
1035     map f(j)=Rr(j+i),p(1..n);
1036     ideal JJ=f(j)(K(j));
1037     if(dim(std(JJ))>i)
1038     {
1039      ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1040      list W(i)=imap(TJ(0),Ww(i));
1041      list w(j)=var(1);
1042     }
1043     else
1044     {
1045      if(i==0)
1046      {
1047       ring RR(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1048       list W(i)=imap(TJ(0),Ww(i));
1049       list w=W(i)[j];
1050       ideal J=imap(TJ(0),N(0));
1051       def AA(j)=JuReZeroDim( J,w,d);
1052       setring AA(j);
1053       list W(i)=imap(TJ(0),Ww(i));
1054       if(t<1)
1055       {
1056        ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1057        list W(i)=imap(TJ(0),Ww(i));
1058        list w(j)=W(i)[j];
1059       }
1060       else
1061       {
1062        ring RRR(j) = create_ring("(complex,106,I)", "("+varstr(S)+")", "dp");
1063        list W(i)=imap(TJ(0),Ww(i));
1064        list w=W(i)[j];
1065        ideal J=imap(TJ(0),N(0));
1066        def AAA(j)=JuReTopDim( J,w, z(d),d);
1067        setring AAA(j);
1068        number ts=t;
1069        list W(i)=imap(TJ(0),Ww(i));
1070        if(ts<1)
1071        {
1072         dt=d-1;
1073        }
1074        else
1075        {
1076         dt=d;
1077        }
1078        if(dt>i)
1079        {
1080         for(ii=i+1;ii<=dt;ii++)
1081         {
1082          ring RRRR(ii+j) = create_ring("(complex,106,I)", "("+varstr(S)+")", "ds");
1083          list Ww(ii)=imap(TJ(0),Ww(ii));
1084          if(size(Ww(ii)[1])>1)
1085          {
1086           ring RRRRR(ii+j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1087           list W(i)=imap(TJ(0),Ww(i));
1088           list w=W(i)[j];
1089           list Ww(ii)=imap(TJ(0),Ww(ii));
1090           ideal J=imap(TJ(0),N(0));
1091           ideal L=imap(TJ(0),LL);
1092           ideal L(ii);
1093           for(k=1;k<=ii;k++)
1094           {
1095           L(ii)[k]=L[k];
1096           }
1097           def AAA(ii+j)=ReJunkUseHomo(J,L(ii),Ww(ii),w);
1098           setring AAA(ii+j);
1099           number ts=t;
1100           list W(i)=imap(TJ(0),Ww(i));
1101           if(ts>0)
1102           {
1103           ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1104           list W(i)=imap(TJ(0),Ww(i));
1105           list w(j)=var(1);
1106           ii=d+1;
1107           }
1108           else
1109           {
1110            if(ii==dt)
1111            {
1112             ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1113             list W(i)=imap(TJ(0),Ww(i));
1114            }
1115            list w(j)=W(i)[j];
1116           }
1117          }
1118         }
1119        }
1120        else
1121        {
1122         ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1123         list W(i)=imap(TJ(0),Ww(i));
1124         list w(j)=W(i)[j];
1125        }
1126       }
1127      }
1128      else
1129      {
1130       ring RRRRRRR(j) = create_ring("(complex,106,I)", "("+varstr(S)+")", "dp");
1131       list W(i)=imap(TJ(0),Ww(i));
1132       list w=W(i)[j];
1133       ideal J=imap(TJ(0),N(0));
1134       def Aaa(j)=JuReTopDim( J,w,z(d),d);
1135       setring Aaa(j);
1136       list W(i)=imap(TJ(0),Ww(i));
1137       number ts =t;
1138       if(ts<1)
1139       {
1140        dt=d-1;
1141       }
1142       else
1143       {
1144        dt=d;
1145       }
1146       if(dt>i)
1147       {
1148        for(ii=i+1;ii<=dt;ii++)
1149        {
1150         ring RRRRRRRR(ii+j) = create_ring("(complex,106,I)", "("+varstr(S)+")", "ds");
1151         list Ww(ii)=imap(TJ(0),Ww(ii));
1152         if(size(Ww(ii)[1])>1)
1153         {
1154          ring R1(ii+j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1155          list W(i)=imap(TJ(0),Ww(i));
1156          list w=W(i)[j];
1157          list Ww(ii)=imap(TJ(0),Ww(ii));
1158          ideal J=imap(TJ(0),N(0));
1159          ideal L=imap(TJ(0),LL);
1160          ideal L(ii);
1161          for(k=1;k<=ii;k++)
1162          {
1163           L(ii)[k]=L[k];
1164          }
1165          def AA(ii+j)=ReJunkUseHomo(J,L(ii),Ww(ii),w);
1166          setring AA(ii+j);
1167          number ts=t;
1168          list W(i)=imap(TJ(0),Ww(i));
1169          if(ts>0)
1170          {
1171           ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1172           list W(i)=imap(TJ(0),Ww(i));
1173           list w(j)=var(1);
1174           ii=d+1;
1175          }
1176          else
1177          {
1178           if(ii==dt)
1179           {
1180           ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1181           list W(i)=imap(TJ(0),Ww(i));
1182           }
1183           list w(j)=W(i)[j];
1184          }
1185         }
1186        }
1187       }
1188       else
1189       {
1190        ring A(j) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1191        list W(i)=imap(TJ(0),Ww(i));
1192        list w(j)=W(i)[j];
1193       }
1194      }
1195     }
1196     if(j==z(i))
1197     {
1198      ring R(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1199      list W(i),w;
1200      int k(i)=0;
1201      for(k=1;k<=z(i);k++)
1202      {
1203       w=imap(A(k),w(k));
1204       if(size(w)>1)
1205       {
1206        k(i)=k(i)+1;
1207        W(i)[k(i)]=w;
1208       }
1209      }
1210      if(k(i)==0)
1211      {
1212       W(i)=var(1);
1213      }
1214     }
1215    }
1216   }
1217  }
1218  ring T = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1219  int bt=0;
1220  for(i=0;i<=d-1;i++)
1221  {
1222   list Ww(i)=imap(TJ(0),W(i));
1223   if(size(Ww(i)[1])>1)
1224   {
1225    list W(i)=imap(R(i),W(i));
1226   }
1227   else
1228   {
1229    list W(i)=Ww(i);
1230   }
1231   export(W(i));
1232  }
1233  list W(d)=imap(TJ(0),W(d));
1234  export(W(d));
1235  ideal L=imap(TJ(0),LL);
1236  export(L);
1237  ideal N(0)=imap(TJ(0),N(0));
1238  export(N(0));
1239  setring S;
1240  return(T);
1241 }
1242}
1243example
1244{ "EXAMPLE:";echo = 2;
1245    ring r=0,(x,y,z),dp;
1246    poly f1=(x3+z)*(x2-y);
1247    poly f2=(x3+y)*(x2-z);
1248    poly f3=(x3+y)*(x3+z)*(z2-y);
1249    ideal I=f1,f2,f3;
1250    def W=WitSet(I);
1251    setring W;
1252    W(1);
1253         // witness point  set of a pure 1-dimensional component of V(I)
1254    W(0);
1255         // witness point  set of a pure 0-dimensional component of V(I)
1256    L;
1257         // list of generic linear polynomials
1258}
1259///////////////////////////////////////////////////////////////////////////////
1260static proc ZSR1(ideal I, ideal L, list W )
1261"USAGE:  ZSR1(ideal I, ideal L, list W );I ideal,
1262         L ideal defined by generic linear polynomials,
1263         W list of a point in the generic slicing of V(I) and V(L)
1264RETURN:  ts number;zero sum relation of W
1265EXAMPLE: example ZSR1;shows an example
1266"
1267{
1268 def S=basering;
1269 int n=nvars(basering);
1270 int c(1)=5*n;
1271 int c(2)=23*n;
1272 int iii=size(L);
1273 ring R = create_ring("(complex,16,I)", "("+varstr(S)+")", "ds");
1274 number c(0);
1275 ideal LL=imap(S,L);
1276 c(0)=leadcoef(LL[iii]);
1277 string sv=varstr(S);
1278 int j,ii,jj,k,a,b,te,zi,si;
1279 string sf,sff;
1280 list VV;
1281 list W=imap(S,W);
1282 VV[1]=W;
1283 def SB1=Singular2bertini(VV);
1284 ring R = create_ring("(complex,16,I)", "("+varstr(S)+",gamma,s)", "dp");
1285 ideal I=imap(S,I);
1286 zi=size(I);
1287 ideal LL=imap(S,L);
1288 ideal H, ll;
1289 for(a=1;a<=zi;a++)
1290 {
1291  H[a]=s*gamma*I[a]+(1-s)*I[a];
1292 }
1293 if(iii>1)
1294 {
1295  for(k=1;k<=iii-1;k++)
1296  {
1297   ll[k]=LL[k];
1298   H[k+zi]=s*gamma*LL[k]+(1-s)*ll[k];
1299  }
1300  ll[iii]=LL[iii]+c(1)-c(0);
1301  H[iii+zi]=s*gamma*LL[iii]+(1-s)*ll[iii];
1302 }
1303 else
1304 {
1305  ll[iii]=LL[iii]+c(1)-c(0);
1306   H[iii+zi]=s*gamma*LL[iii]+(1-s)*ll[iii];
1307 }
1308 def Q(1)=UseBertini(H,sv);
1309 string siaa=read("singular_solutions");
1310 string saa=read("nonsingular_solutions");
1311 def TT(1)=bertini2Singular("nonsingular_solutions",nvars(basering)-2);
1312 setring TT(1);
1313 list wr=re;
1314 if(size(wr)==0)
1315 {
1316  ring TT(2) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1317  number tte, ts;
1318  tte=11;
1319  ts=0;
1320  export(ts);
1321  export(tte);
1322 }
1323 else
1324 {
1325  ring R1 = create_ring("(complex,16,I)", "("+varstr(S)+",gamma,s)", "dp");
1326  ideal I=imap(S,I);
1327  si=size(I);
1328  ideal LL=imap(S,L);
1329  ideal H, ll;
1330  for(a=1;a<=si;a++)
1331  {
1332   H[a]=s*gamma*I[a]+(1-s)*I[a];
1333  }
1334  if(iii>1)
1335  {
1336   for(k=1;k<=iii-1;k++)
1337   {
1338    ll[k]=LL[k];
1339    H[k+si]=s*gamma*LL[k]+(1-s)*ll[k];
1340   }
1341   ll[iii]=LL[iii]+c(2)-c(0);
1342   H[iii+si]=s*gamma*LL[iii]+(1-s)*ll[iii];
1343  }
1344  else
1345  {
1346   ll[iii]=LL[iii]+c(2)-c(0);
1347   H[iii+si]=s*gamma*LL[iii]+(1-s)*ll[iii];
1348  }
1349  def Q(2)=UseBertini(H,sv);
1350  string saaa=read("nonsingular_solutions");
1351  string siaaa=read("singular_solutions");
1352  if(size(saaa)<52)
1353  {
1354   if(size(siaaa)<52)
1355   {
1356    "ERROR( Try again try);";
1357   }
1358  }
1359  if(size(saaa)>=52)
1360  {
1361   def TT(2)=bertini2Singular("nonsingular_solutions",nvars(basering)-2);
1362   setring TT(2);
1363   list wwr=re;
1364   list wr=imap(TT(1),wr);
1365   list W=imap(S,W);
1366   list w,ww,www;
1367   number s(0),s(1),s(2),ts;
1368   zi=size(W)/n;
1369   for(jj=1;jj<=zi;jj++)
1370   {
1371    s(0)=0;
1372    s(1)=0;
1373    s(2)=0;
1374    w=W;
1375    ww=wr[jj];
1376    www=wwr[jj];
1377    for(j=1;j<=n;j++)
1378    {
1379     s(0)=s(0)+j*w[j];
1380    }
1381    for(j=1;j<=n;j++)
1382    {
1383     s(1)=s(1)+j*ww[j];
1384    }
1385    for(j=1;j<=n;j++)
1386    {
1387     s(2)=s(2)+j*www[j];
1388    }
1389   }
1390   ts=s(0)*(c(1)-c(2))+s(1)*(c(2)-c(0))+s(2)*(c(0)-c(1));
1391  }
1392  else
1393  {
1394   def TT(2)=bertini2Singular("singular_solutions",nvars(basering)-2);
1395   setring TT(2);
1396   list wwr=re;
1397   list wr=imap(TT(1),wr);
1398   list W=imap(S,W);
1399   list w,ww,www;
1400   number s(0),s(1),s(2),ts;
1401   zi=size(W)/n;
1402   for(jj=1;jj<=zi;jj++)
1403   {
1404    s(0)=0;
1405    s(1)=0;
1406    s(2)=0;
1407    w=W;
1408    ww=wr[jj];
1409    www=wwr[jj];
1410    for(j=1;j<=n;j++)
1411    {
1412     s(0)=s(0)+j*w[j];
1413    }
1414    for(j=1;j<=n;j++)
1415    {
1416     s(1)=s(1)+j*ww[j];
1417    }
1418    for(j=1;j<=n;j++)
1419    {
1420     s(2)=s(2)+j*www[j];
1421    }
1422   }
1423   ts=s(0)*(c(1)-c(2))+s(1)*(c(2)-c(0))+s(2)*(c(0)-c(1));
1424  }
1425 }
1426 ring e = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1427 number ts=imap(TT(2),ts);
1428 export(ts);
1429 number tte;
1430 tte=11;
1431 export(tte);
1432 system("sh","rm start");
1433 setring S;
1434 return (e);
1435}
1436example
1437{ "EXAMPLE:";echo = 2;
1438    ring r=(complex,16,I),(x,y,z),dp;
1439    poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
1440    poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
1441    poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
1442    ideal J=f1,f2,f3;
1443    ideal L=2*x+7*y+3*z+29;
1444    list W=2,1.999999999999999,-15.6666666666664;
1445    def D=ZSR1(J,L,W );
1446    setring D;
1447    ts;
1448}
1449///////////////////////////////////////////////////////////////////////////////
1450static proc perSumZ(list A)
1451"USAGE:  perSumZ(list A);A list of different complex numbers
1452RETURN:  all subsets of A, whose sum of their elements is zero
1453EXAMPLE: example perSumZ;shows an example
1454"
1455{
1456 list B, C;
1457 int i, j;
1458 number t,tr;
1459 if(size(A)==0)
1460 {
1461  B[1]=A;
1462 }
1463 if(size(A)==1)
1464 {
1465  if(((repart(A[1]))^2+(impart(A[1]))^2)<=0.000000000000001)
1466  {
1467   B[1]=A;
1468  }
1469 }
1470 for(i=1;i<=size(A);i++)
1471 {
1472  t=t+A[i];
1473 }
1474 if(((repart(t))^2+(impart(t))^2)<=0.000000000000001)
1475 {
1476  B[1]=A;
1477 }
1478 for(i=1;i<=size(A);i++)
1479 {
1480  C=delete(A,i);
1481  C=perSumZ(C);
1482  for(j=1;j<=size(C);j++)
1483  {
1484   if(size(C[j])>0)
1485   {
1486    B[size(B)+1]=C[j];
1487   }
1488  }
1489 }
1490 return(B);
1491}
1492example
1493{ "EXAMPLE:";echo = 2;
1494   ring r=(complex,16,I),x,lp;
1495   list A=1,-1,2-I,I,-2;
1496   def D=perSumZ(A);
1497   D;
1498}
1499///////////////////////////////////////////////////////////////////////////////
1500static proc  ZSROFWitSet(ideal I)
1501"USAGE:  ZSROFWitSet(ideal I);I ideal
1502RETURN:  ZSR(i) lists of the zero sum relation of witness point
1503          sets W(i) for i=1,...dim(V(I))
1504EXAMPLE: example ZSROFWitSet;shows an example
1505"
1506{
1507 def S=basering;
1508 int n=nvars(basering);
1509 def T(0)=WitSet(I);
1510 setring T(0);
1511 ideal LL=L;
1512 int dd=size(LL);
1513 int a=c(0);
1514 if(a==0)
1515 {
1516  return(T(0));
1517 }
1518 else
1519 {
1520  int i,j,ii,jj,k,sv(0..dd),j(0..dd),kk;
1521  string sv;
1522  for(i=1;i<=dd;i++)
1523  {
1524   jj=0;
1525   list V(i)=imap(T(0),W(i));
1526   if(size(V(i)[1])>1)
1527   {
1528    if(size(V(i))==1)
1529    {
1530     ring L(i)(1) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1531     list W(i),ZSR(i);
1532     list V(i)=imap(T(0),W(i));
1533     W(i)=V(i);
1534     ZSR(i)[1]=0.000;
1535    }
1536    else
1537    {
1538     if(i>1)
1539     {
1540      ring ee(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1541      list V(i)=imap(T(0),W(i));
1542     }
1543     sv(i)=size(V(i));
1544     for(j=1;j<=sv(i);j++)
1545     {
1546      ideal N=imap(T(0),N(0));
1547      ideal LLL=imap(T(0),LL);
1548      ideal L;
1549      for(kk=1;kk<=i;kk++)
1550      {
1551       L[kk]=LLL[kk];
1552      }
1553      def L(i)(j)=ZSR1(N,L,V(i)[j]);
1554      setring L(i)(j);
1555      if(j==1)
1556      {
1557       list W(i),ZSR(i);
1558      }
1559      else
1560      {
1561       list W(i)=imap(L(i)(j-1),W(i));
1562       list ZSR(i)=imap(L(i)(j-1),ZSR(i));
1563       export(ZSR(i));
1564      }
1565      list V(i)=imap(T(0),W(i));
1566      jj=jj+1;
1567      ZSR(i)[jj]=ts;
1568      W(i)[jj]=V(i)[j];
1569     }
1570    }
1571   }
1572  }
1573  ring Q = create_ring("(complex,12,I)", "("+varstr(S)+")", "dp");
1574  list W(0)=imap(T(0),W(0));
1575  export(W(0));
1576  for(jj=1;jj<=dd;jj++)
1577  {
1578   number pt(jj);
1579   list V(jj)=imap(T(0),W(jj));
1580   if(size(V(jj)[1])>1)
1581   {
1582    sv(jj)=size(V(jj));
1583    if(jj>0)
1584    {
1585     list ZSR(jj)=imap(L(jj)(sv(jj)),ZSR(jj));
1586     export(ZSR(jj));
1587     list W(jj)=imap(L(jj)(sv(jj)),W(jj));
1588     export(W(jj));
1589    }
1590   }
1591   else
1592   {
1593    list ZSR(jj)=var(1);
1594    export(ZSR(jj));
1595    list W(jj)=var(1);
1596    export(W(jj));
1597    }
1598   }
1599   ideal L=imap(T(0),LL);
1600   export(L);
1601   export(dd);
1602   system("sh","rm singular_solutions");
1603   system("sh","rm nonsingular_solutions");
1604   system("sh","rm real_solutions");
1605   system("sh","rm raw_solutions");
1606   system("sh","rm raw_data");
1607   system("sh","rm output");
1608   system("sh","rm midpath_data");
1609   system("sh","rm main_data");
1610   system("sh","rm input");
1611   system("sh","rm failed_paths");
1612   setring S;
1613   return(Q);
1614  }
1615}
1616example
1617{ "EXAMPLE:";echo = 2;
1618    ring  r = 0,(x,y,z),dp;
1619    poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
1620    poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
1621    poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
1622    ideal J=f1,f2,f3;
1623    def D=ZSROFWitSet(J);
1624    setring D;
1625    ZSR(1);
1626    W(1);
1627    ZSR(2);
1628    W(2);
1629}
1630///////////////////////////////////////////////////////////////////////////////
1631static proc ReWitZSR(list A, list W, int di)
1632"USAGE:  ReWitZSR(list A, list W, int di);A ideal of complex numbers,
1633          W list of points on di-dimensional component,
1634         di integer
1635RETURN:  tw(di) integer, list Z(size(Z));
1636        if tw(di)>0, else Z(0), list Z1(tw1(di))
1637EXAMPLE: example ReWitZSR;shows an example
1638"
1639{
1640 def S=basering;
1641 ring e = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1642 list A=imap(S,A);
1643 list W=imap(S,W);
1644 list  D, B(1..size(A)),C(1..size(A)),D(1..size(A)),Z1(1..size(A)),Z(1..size(A)),Z,Y,ZY;
1645 int i,j,k,tw(di),tw1(di),tw2(di),tw3(di),tr,tc;
1646 list AA;
1647 list WW;
1648 for(i=1;i<=size(A);i++)
1649 {
1650  if(((repart(A[i]))^2+(impart(A[i]))^2)<=0.0000000000000001)
1651  {
1652   tc=tc+1;
1653   tw1(di)=tw1(di)+1;
1654   ZY[i]=A[i];
1655   Z1(tw1(di))=W[i];
1656   export(Z1(tw1(di)));
1657  }
1658  else
1659  {
1660   tr=tr+1;
1661   tw2(di)=tw2(di)+1;
1662   AA[tr]=A[i];
1663   WW[tr]=W[i];
1664  }
1665 }
1666 A=AA;
1667 W=WW;
1668 if(size(A)>0)
1669 {
1670  def B=perSumZ(A);
1671  for(i=1;i<=size(A);i++)
1672  {
1673   tc=0;
1674   B(i)=A[i];
1675   for(j=1;j<=size(B);j++)
1676   {
1677    tr=0;
1678    for(k=1;k<=size(B[j]);k++)
1679    {
1680     if(B(i)[1]==B[j][k])
1681     {
1682      tr=tr+1;
1683     }
1684    }
1685    if(tr>0)
1686    {
1687     tc=tc+1;
1688     C(i)[tc]=B[j];
1689    }
1690   }
1691   for(j=1;j<=size(C(i));j++)
1692   {
1693    D(i)=C(i)[j];
1694    for(k=1;k<=size(C(i));k++)
1695    {
1696     if(size(D(i))<size(C(i)[k]))
1697     {
1698      D(i)=D(i);
1699     }
1700     else
1701     {
1702      D(i)=C(i)[k];
1703     }
1704    }
1705   }
1706  }
1707  for(i=1;i<=size(A);i++)
1708  {
1709   Z[i]=D(i);
1710  }
1711  for(i=1;i<=size(Z);i++)
1712  {
1713   if(size(Z[i])>0)
1714   {
1715    D=Z[i];
1716    for(k=size(Z);k>0;k--)
1717    {
1718     if(size(Z[k])>0)
1719     {
1720      B=Z[k];
1721      if(i!=k)
1722      {
1723       if(D[1]==B[1])
1724       {
1725        Z=delete(Z,k);
1726       }
1727      }
1728     }
1729    }
1730   }
1731  }
1732  for(j=1;j<=size(Z);j++)
1733  {
1734    tr=0;
1735    D=Z[j];
1736    for(i=1;i<=size(A);i++)
1737    {
1738     for(k=1;k<=size(D);k++)
1739     {
1740      if(A[i]==D[k])
1741      {
1742       tr=tr+1;
1743       tw(di)=tw(di)+1;
1744       Z(j)[tr]=W[i];
1745      }
1746     }
1747    }
1748    export(Z(j));
1749  }
1750 export(Z);
1751 }
1752 if(tw1(di)==0)
1753 {
1754  list Z1(0);
1755  Z1(0)="Empty set";
1756  export(Z1(0));
1757 }
1758 if(tw(di)==0)
1759 {
1760  list Z(0);
1761  Z(0)="Empty set";
1762  export(Z(0));
1763 }
1764 export(tw1(di));
1765 export(tw(di));
1766  setring S;
1767  return (e);
1768}
1769example
1770{ "EXAMPLE:";echo = 2;
1771    ring r=(complex,16,I),(x,y,z),dp;
1772    list A= 3.7794571034732007+I*21.1724850800421247,
1773           -3.7794571034752664-I*21.1724850800419908;
1774    list W=list(-2.0738016397747976,1.29520655909919,-0.1476032795495952),
1775           list(-1.354769788796631,-1.5809208448134761,1.2904604224067381);
1776    int di=1;
1777    def D=ReWitZSR(A,W,di);
1778    setring D;
1779    tw(di);
1780    Z(size(Z));// if tw(di)>0, else Z(0);
1781    Z1(tw1(di));
1782}
1783///////////////////////////////////////////////////////////////////////////////
1784proc NumIrrDecom(ideal I) Numerical Irreducible Decomposition
1785"USAGE:  NumIrrDecom(ideal I);I ideal
1786RETURN:  w(1),..., w(t) lists of irreducible witness point sets of
1787         irreducible components of V(J)
1788EXAMPLE: example NumIrrDecom;shows an example
1789"
1790{
1791 def S=basering;
1792 int i,ii;
1793 def WW=ZSROFWitSet(I);
1794 setring WW;
1795 if(c(0)==0)
1796 {
1797  setring S;
1798  return(WW);
1799 }
1800 else
1801 {
1802  int d=size(L);
1803  for(i=0;i<=d;i++)
1804  {
1805   int co(i)=0;
1806   if(i==0)
1807   {
1808    ring q(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1809    list V(i)=imap(WW,W(i));
1810    list W(0..size(V(i)));
1811    if(size(V(i)[1])>1)
1812    {
1813     co(i)=size(V(i));
1814     for(ii=1;ii<=size(V(i));ii++)
1815     {
1816      list w(ii)=V(i)[ii];
1817      export(w(ii));
1818     }
1819    }
1820    else
1821    {
1822     W(1)[1]="Empty Set";
1823    }
1824   }
1825   else
1826   {
1827    list WW(i);
1828    list V(i)=imap(WW,W(i));
1829    list a(i)=imap(WW,ZSR(i));
1830    if(size(V(i)[1])>1)
1831    {
1832     def q(i)=ReWitZSR(a(i),V(i),i);
1833     setring q(i);
1834     if(tw1(i)>0)
1835     {
1836      for(ii=1;ii<=tw1(i);ii++)
1837      {
1838       WW(i)[ii]=Z1(ii);
1839      }
1840      co(i)=tw1(i);
1841     }
1842     if(tw(i)>0)
1843     {
1844      for(ii=1;ii<=size(Z);ii++)
1845      {
1846       if(size(Z[ii])>1)
1847       {
1848        co(i)=co(i)+1;
1849        WW(i)[ii+tw1(i)]=Z(ii);
1850       }
1851      }
1852     }
1853     for(ii=1;ii<=size(WW(i));ii++)
1854     {
1855      list w(ii);
1856      w(ii)=WW(i)[ii];
1857     }
1858    }
1859    else
1860    {
1861     ring q(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1862     WW(i)[1]="Empty Set";
1863    }
1864   }
1865  }
1866  for(i=0;i<=d;i++)
1867  {
1868   ring qq(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
1869   for(ii=1;ii<=co(i);ii++)
1870   {
1871    list w(ii)=imap(q(i),w(ii));
1872    export w(ii);
1873   }
1874  "===========================================";
1875  "===========================================";
1876   "Dimension";
1877    i;
1878   "Number of Components";
1879    co(i);
1880   number cco(i)=co(i)/1;
1881   export(cco(i));
1882  }
1883  ideal L=imap(WW,L);
1884  export(L);
1885  "The generic Linear Space L";
1886  L;
1887  return(qq(0..d));
1888 }
1889}
1890example
1891{ "EXAMPLE:";echo = 2;
1892    ring r=0,(x,y,z),dp;
1893    poly f1=(x2+y2+z2-6)*(x-y)*(x-1);
1894    poly f2=(x2+y2+z2-6)*(x-z)*(y-2);
1895    poly f3=(x2+y2+z2-6)*(x-y)*(x-z)*(z-3);
1896    ideal I=f1,f2,f3;
1897    list W=NumIrrDecom(I);
1898      ==>
1899         Dimension
1900         0
1901         Number of Components
1902         1
1903         Dimension
1904         1
1905         Number of Components
1906         3
1907         Dimension
1908         2
1909         Number of Components
1910         1
1911    def A(0)=W[1];
1912       // corresponding 0-dimensional components
1913    setring A(0);
1914    w(1);
1915        // corresponding 0-dimensional irreducible component
1916         ==> 0-Witness point set (one point)
1917    def A(1)=W[2];
1918           // corresponding 1-dimensional components
1919    setring A(1);
1920    w(1);
1921         // corresponding 1-dimensional irreducible component
1922         ==> 1-Witness point set (one point)
1923    w(2);
1924       // corresponding 1-dimensional irreducible component
1925         ==> 1-Witness point set (one point)
1926    w(3);
1927      // corresponding 1-dimensional irreducible component
1928        ==> 1-Witness point set (one point)
1929    def A(2)=W[3];
1930    // corresponding 2-dimensional components
1931    setring A(2);
1932    w(1);
1933      // corresponding 2-dimensional irreducible component
1934        ==> 1-Witness point set (two points)
1935}
1936///////////////////////////////////////////////////////////////////////////////
1937proc defl(ideal I, int d)
1938"USAGE:   defl(ideal I, int d);  I ideal, int d order of the deflation
1939RETURN:   deflation ideal DI of I
1940EXAMPLE:  example defl; shows an example
1941"
1942{
1943 def S=basering;
1944 int n=nvars(basering);
1945 int i,j;
1946 for(i=1;i<=d;i++)
1947 {
1948  def R(i)=symmetricBasis(n,i,"x");
1949  setring R(i);
1950  ideal J(i)=symBasis;
1951  export(J(i));
1952 }
1953 execute("ring RR=0,(x(1..n),"+varstr(S)+"),dp;");
1954 for(i=1;i<=d;i++)
1955 {
1956  ideal J(i)=imap(R(i),J(i));
1957  for(j=1;j<=n;j++)
1958  {
1959   J(i)=subst(J(i),x(j),var(n+j));
1960  }
1961 }
1962 ring R = create_ring(0, "("+varstr(S)+")", "dp");
1963 ideal I=imap(S,I);
1964 if(d>1)
1965 {
1966  for(i=1;i<=d-1;i++)
1967  {
1968   ideal J(i)=imap(RR,J(i));
1969   for(j=1;j<=size(I);j++)
1970   {
1971    ideal I(j);
1972    for(k=1;k<=size(J(i));k++)
1973    {
1974     I(j)[k]=J(i)[k]*I[j];
1975    }
1976     export(I(j));
1977   }
1978  }
1979  ideal J(d)=imap(RR,J(d));
1980  ideal D(d)=J(1..d);
1981  ideal II(d)=I,I(1..size(I));
1982  matrix T(d)=diff(D(d),II(d));
1983  matrix TT(d)=transpose(T(d));
1984  export(TT(d));
1985 }
1986 else
1987 {
1988  ideal J(d)=imap(RR,J(d));
1989  ideal D(d)=J(d);
1990  ideal II(d)=I;
1991  matrix T(d)=diff(D(d),II(d));
1992  matrix TT(d)=transpose(T(d));
1993  export(TT(d));
1994 }
1995 int zc=size(D(d));
1996 export(zc);
1997 execute("ring DR=0,("+varstr(S)+",x(1..zc)),dp;");
1998 matrix TT(d)=imap(R,TT(d));
1999 ideal I=imap(S,I);
2000 vector v=[x(1..zc)];
2001 ideal DI=I,TT(d)*v;
2002 export(DI);
2003 export(I);
2004 setring S;
2005 return(DR);
2006}
2007example
2008{ "EXAMPLE:"; echo = 2;
2009   ring r=0,(x,y,z),dp;
2010   poly f1=z^2;
2011   poly f2=z*(x^2+y);
2012   ideal I=f1,f2;
2013   def D=defl(I,1);
2014   setring D;
2015   DI;
2016}
2017///////////////////////////////////////////////////////////////////////////////
2018static proc     NIDofDI(ideal I)
2019"USAGE:  NIDofDI(ideal I);  I ideal
2020RETURN:  numerical irreducible decomposition of I
2021EXAMPLE: NIDofDI; shows an example
2022"
2023{
2024 def S=basering;
2025 int i,ii;
2026 def WW=ZSROFWitSet(I);
2027 setring WW;
2028 if(c(0)==0)
2029 {
2030  setring S;
2031  return(WW);
2032 }
2033 else
2034 {
2035  int d=size(L);
2036  for(i=0;i<=d;i++)
2037  {
2038   int co(i)=0;
2039   if(i==0)
2040   {
2041    ring q(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
2042    list V(i)=imap(WW,W(i));
2043    list W(0..size(V(i)));
2044    if(size(V(i)[1])>1)
2045    {
2046     co(i)=size(V(i));
2047     for(ii=1;ii<=size(V(i));ii++)
2048     {
2049      list w(ii)=V(i)[ii];
2050      export(w(ii));
2051     }
2052    }
2053    else
2054    {
2055     W(1)[1]="Empty Set";
2056    }
2057   }
2058   else
2059   {
2060    list WW(i);
2061    list V(i)=imap(WW,W(i));
2062    list a(i)=imap(WW,ZSR(i));
2063    if(size(V(i)[1])>1)
2064    {
2065     def q(i)=ReWitZSR(a(i),V(i),i);
2066     setring q(i);
2067     if(tw1(i)>0)
2068     {
2069      for(ii=1;ii<=tw1(i);ii++)
2070      {
2071       WW(i)[ii]=Z1(ii);
2072      }
2073      co(i)=tw1(i);
2074     }
2075     if(tw(i)>0)
2076     {
2077      for(ii=1;ii<=size(Z);ii++)
2078      {
2079       if(size(Z[ii])>1)
2080       {
2081        co(i)=co(i)+1;
2082        WW(i)[ii+tw1(i)]=Z(ii);
2083       }
2084      }
2085     }
2086     for(ii=1;ii<=size(WW(i));ii++)
2087     {
2088      list w(ii);
2089      w(ii)=WW(i)[ii];
2090     }
2091    }
2092    else
2093    {
2094     ring q(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
2095     WW(i)[1]="Empty Set";
2096    }
2097   }
2098  }
2099  for(i=0;i<=d;i++)
2100  {
2101   ring qq(i) = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
2102   list ww(i);
2103   if(co(i)>0)
2104   {
2105    for(ii=1;ii<=co(i);ii++)
2106    {
2107     list v(ii)=imap(q(i),w(ii));
2108     ww(i)=v(ii)[1];
2109     if(size(ww(i))==1)
2110     {
2111      list w(ii);
2112      w(ii)[1]=v(ii);
2113     }
2114     else
2115     {
2116      list w(ii)=v(ii);
2117     }
2118     export(w(ii));
2119    }
2120   }
2121   else
2122   {
2123    list w(1);
2124    w(1)[1]=var(1);
2125    export(w(1));
2126   }
2127   number cco(i)=co(i)/1;
2128   export(cco(i));
2129  }
2130  ideal L=imap(WW,L);
2131  export(L);
2132  return(qq(0..d));
2133 }
2134}
2135example
2136{ "EXAMPLE:"; echo = 2;
2137   ring r=0,(x,y,z),dp;
2138   poly f1=z^2;
2139   poly f2=z*(x^2+y);
2140   ideal I=f1,f2;
2141   list DD=NIDofDI(I);
2142   def D(0)=DD[1];
2143   setring D(0);
2144   w(1);          // w(1)= x, i.e. no components
2145   def D(1)=DD[2];
2146   setring D(1);
2147   w(1);
2148   def D(2)=DD[3];
2149   setring D(2);
2150   w(1);
2151}
2152///////////////////////////////////////////////////////////////////////////////
2153proc NumPrimDecom(ideal I,int d)
2154"USAGE:   NumPrimDecom(ideal I,int d);  I ideal, d order of the deflation
2155RETURN:   lists of the numerical primary decomposition
2156EXAMPLE:  example NumPrimDecom; shows an example
2157"
2158{
2159 def S=basering;
2160 int n=nvars(basering);
2161 int i,Dd,j,k,jj;
2162 def D=defl(I,d);
2163 setring D;
2164 ideal J=DI;
2165 Dd=dim(std(DI));
2166 list W=NIDofDI(J);
2167 for(i=0;i<=Dd;i++)
2168 {
2169  def A(i+1)=W[i+1];
2170  setring A(i+1);
2171  if(cco(i)>0)
2172  {
2173   for(j=1;j<=cco(i);j++)
2174   {
2175    list W(j)=w(j);
2176    for(k=1;k<=size(w(j));k++)
2177    {
2178     for(jj=size(W(j)[k]);jj>=n+1;jj--)
2179     {
2180      W(j)[k]=delete(W(j)[k],jj);
2181     }
2182     W(j)[k]=W(j)[k];
2183    }
2184   }
2185  }
2186  else
2187  {
2188   list W(1)=var(1);
2189  }
2190 }
2191 ring R = create_ring("(complex,16,I)", "("+varstr(S)+")", "dp");
2192 jj=0;
2193 for(i=0;i<=Dd;i++)
2194 {
2195  number cco(i)=imap(A(i+1),cco(i));
2196  if(cco(i)>0)
2197  {
2198   for(j=1;j<=cco(i);j++)
2199   {
2200    jj=jj+1;
2201    list w(jj)=imap(A(i+1),W(j));
2202    export(w(jj));
2203  "===========================================";
2204  "===========================================";
2205    "Numerical Primary Component";
2206     w(jj);
2207   }
2208  }
2209 }
2210 return(R);
2211}
2212example
2213{ "EXAMPLE:"; echo = 2;
2214   ring r=0,(x,y),dp;
2215   poly f1=yx;
2216   poly f2=x2;
2217   ideal I=f1,f2;
2218   def W=NumPrimDecom(I,1);
2219   setring W;
2220   w(1);
2221      ==> 1-Witness point set (one point)
2222   w(2);
2223      ==> 1-Witness point set (one point)
2224}
2225///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.