1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id$"; |
---|
3 | category="Algebraic Geometry"; |
---|
4 | info=" |
---|
5 | LIBRARY: NumDecom.lib Numerical Decomposition of Ideals |
---|
6 | OVERVIEW: |
---|
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). |
---|
12 | AUTHOR: Shawki AlRashed, rashed@mathematik.uni-kl.de; sh.shawki@yahoo.de |
---|
13 | |
---|
14 | PROCEDURES: |
---|
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 | |
---|
42 | LIB "solve.lib"; |
---|
43 | LIB "matrix.lib"; |
---|
44 | |
---|
45 | /////////////////////////////////////////////////////////////////////////////// |
---|
46 | /////////////////////////////////////////////////////////////////////////////// |
---|
47 | proc re2squ(ideal I) |
---|
48 | "USAGE: re2squ(ideal I);I ideal |
---|
49 | RETURN: ideal J defined by the polynomial system of the same number of polynomials and unknowns |
---|
50 | EXAMPLE: 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 | } |
---|
96 | example |
---|
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 | |
---|
106 | proc Singular2bertini(list L) |
---|
107 | "USAGE: Singular2bertini(list L);L a list |
---|
108 | RETURN: text file called start |
---|
109 | NOTE: adopting the list L to be as a start solution of the homptopy function in Bertini |
---|
110 | EXAMPLE: 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 | } |
---|
131 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
138 | proc UseBertini(ideal H,string sv) |
---|
139 | "USAGE: UseBertini(ideal H,string sv); |
---|
140 | H ideal, sv string of the variable of ring |
---|
141 | RETURN: text file called input used in Bertini to compute the solution of |
---|
142 | the homotopy function H that existed in file input.text |
---|
143 | NOTE: Need to define a start solution of H |
---|
144 | EXAMPLE: 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 | } |
---|
214 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
222 | proc bertini2Singular(string snp, int q) |
---|
223 | "USAGE: bertini2Singular(string snp, int q); |
---|
224 | snp string, q=nvars(basering) integer |
---|
225 | RETURN: re list of the solutions of the homotopy function computed by Bertini |
---|
226 | EXAMPLE: 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 | } |
---|
281 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
289 | proc WitSupSet(ideal I) |
---|
290 | "USAGE: WitSupSet(ideal I);I ideal |
---|
291 | RETURN: 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 |
---|
294 | EXAMPLE: 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 | } |
---|
664 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
683 | proc 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) |
---|
688 | RETURN: t=1 if w on an i-dimensional component of V(I), |
---|
689 | otherwise t=0. Where i=size(L) |
---|
690 | EXAMPLE: 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 | } |
---|
771 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
788 | proc 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) |
---|
791 | RETURN: t=1 if w on a d-dimensional component of V(I), otherwise t=0. |
---|
792 | EXAMPLE: 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 | } |
---|
864 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
877 | proc 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) |
---|
880 | RETURN: t=1 if w on a positive-dimensional component of V(I), |
---|
881 | i.e w is not isolated point in V(J) |
---|
882 | EXAMPLE: 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 | } |
---|
965 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
978 | proc WitSet(ideal I) |
---|
979 | "USAGE: WitSet(ideal I); I ideal |
---|
980 | RETURN: 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 |
---|
983 | NOTE: if W(i)=x, then V(J) has no component of dimension i |
---|
984 | EXAMPLE: 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 | } |
---|
1241 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
1258 | static 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) |
---|
1262 | RETURN: ts number;zero sum relation of W |
---|
1263 | EXAMPLE: 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 | } |
---|
1434 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
1448 | static proc perSumZ(list A) |
---|
1449 | "USAGE: perSumZ(list A);A list of different complex numbers |
---|
1450 | RETURN: all subsets of A, whose sum of their elements is zero |
---|
1451 | EXAMPLE: 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 | } |
---|
1490 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
1498 | static proc ZSROFWitSet(ideal I) |
---|
1499 | "USAGE: ZSROFWitSet(ideal I);I ideal |
---|
1500 | RETURN: ZSR(i) lists of the zero sum relation of witness point |
---|
1501 | sets W(i) for i=1,...dim(V(I)) |
---|
1502 | EXAMPLE: 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 | } |
---|
1614 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
1629 | static 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 |
---|
1633 | RETURN: tw(di) integer, list Z(size(Z)); |
---|
1634 | if tw(di)>0, else Z(0), list Z1(tw1(di)) |
---|
1635 | EXAMPLE: 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 | } |
---|
1767 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
1782 | proc NumIrrDecom(ideal I) Numerical Irreducible Decomposition |
---|
1783 | "USAGE: NumIrrDecom(ideal I);I ideal |
---|
1784 | RETURN: w(1),..., w(t) lists of irreducible witness point sets of |
---|
1785 | irreducible components of V(J) |
---|
1786 | EXAMPLE: 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 | } |
---|
1888 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
1935 | proc defl(ideal I, int d) |
---|
1936 | "USAGE: defl(ideal I, int d); I ideal, int d order of the deflation |
---|
1937 | RETURN: deflation ideal DI of I |
---|
1938 | EXAMPLE: 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 | } |
---|
2005 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
2016 | static proc NIDofDI(ideal I) |
---|
2017 | "USAGE: NIDofDI(ideal I); I ideal |
---|
2018 | RETURN: numerical irreducible decomposition of I |
---|
2019 | EXAMPLE: 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 | } |
---|
2133 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|
2151 | proc NumPrimDecom(ideal I,int d) |
---|
2152 | "USAGE: NumPrimDecom(ideal I,int d); I ideal, d order of the deflation |
---|
2153 | RETURN: lists of the numerical primary decomposition |
---|
2154 | EXAMPLE: 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 | } |
---|
2210 | example |
---|
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 | /////////////////////////////////////////////////////////////////////////////// |
---|