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

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