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