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

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