source: git/Singular/LIB/assprimeszerodim.lib @ 731e4df

spielwiese
Last change on this file since 731e4df was 731e4df, checked in by Hans Schoenemann <hannes@…>, 14 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13280 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 22.5 KB
Line 
1////////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category = "Commutative Algebra";
4info="
5LIBRARY:  assprimeszerodim.lib   associated primes of a zero-dimensional ideal
6
7AUTHORS:  N. Idrees       nazeranjawwad@gmail.com
8@*        G. Pfister      pfister@mathematik.uni-kl.de
9@*        S. Steidel      steidel@mathematik.uni-kl.de
10
11OVERVIEW:
12
13  A library for computing the associated primes and the radical of a
14  zero-dimensional ideal in the polynomial ring over the rational numbers,
15  Q[x_1,...,x_n], using modular computations.
16
17SEE ALSO: primdec_lib
18
19PROCEDURES:
20 zeroR(I);             computes the radical of I
21 assPrimes(I);         computes the associated primes of I
22";
23
24LIB "primdec.lib";
25LIB "modstd.lib";
26
27////////////////////////////////////////////////////////////////////////////////
28
29proc zeroR(ideal I, list #)
30"USAGE:  zeroR(I,[n]); I ideal, optional: n number of processors (for parallel computing)
31ASSUME:  I is zero-dimensional in Q[variables]
32NOTE:    Parallelization is just applicable using 32-bit Singular version since
33         MP-links are not compatible with 64-bit Singular version.
34RETURN:  the radical of I
35EXAMPLE: example zeroR; shows an example
36"
37{
38   attrib(I,"isSB",1);
39   int i, k;
40   int j = 1;
41   int index = 1;
42   int crit;
43
44   list CO1, CO2, P;
45   ideal G, F;
46   bigint N;
47   poly f;
48
49//---------------------  Initialize optional parameter  ------------------------
50   if(size(#) > 0)
51   {
52      int n = #[1];
53      if((n > 1) && (1 - system("with","MP")))
54      {
55         "========================================================================";
56         "There is no MP available on your system. Since this is necessary to     ";
57         "parallelize the algorithm, the computation will be done without forking.";
58         "========================================================================";
59         n = 1;
60      }
61   }
62   else
63   {
64      int n = 1;
65   }
66
67//--------------------  Initialize the list of primes  -------------------------
68   intvec L = primeList(I,10);
69   L[5] = prime(random(100000000,536870912));
70
71   if(n > 1)
72   {
73
74//-----  Create n links l(1),...,l(n), open all of them and compute  -----------
75//-----  polynomial F for the primes L[2],...,L[n + 1].              -----------
76
77      for(i = 1; i <= n; i++)
78      {
79         link l(i) = "MPtcp:fork";
80         open(l(i));
81         write(l(i), quote(zeroRadP(eval(I), eval(L[i + 1]))));
82      }
83
84      int t = timer;
85      P = zeroRadP(I, L[1]);
86      t = timer - t;
87      if(t > 60) { t = 60; }
88      int i_sleep = system("sh", "sleep "+string(t));
89      CO1[index] = P[1];
90      CO2[index] = bigint(P[2]);
91      index++;
92
93      j = j + n + 1;
94   }
95
96//---------  Main computations in positive characteristic start here  ----------
97
98   while(!crit)
99   {
100      if(n > 1)
101      {
102         while(j <= size(L) + 1)
103         {
104            for(i = 1; i <= n; i++)
105            {
106               if(status(l(i), "read", "ready"))
107               {
108                  P = read(l(i));                                        // read the result from l(i)
109                  CO1[index] = P[1];
110                  CO2[index] = bigint(P[2]);
111                  index++;
112
113                  if(j <= size(L))
114                  {
115                     write(l(i), quote(zeroRadP(eval(I), eval(L[j]))));
116                     j++;
117                  }
118                  else
119                  {
120                     k++;
121                     close(l(i));
122                  }
123               }
124            }
125            if(k == n)                                                   // k describes the number of closed links
126            {
127               j++;
128            }
129            i_sleep = system("sh", "sleep "+string(t));                  // sleep for t seconds
130         }
131      }
132      else
133      {
134         while(j<=size(L))
135         {
136            P = zeroRadP(I, L[j]);
137            CO1[index] = P[1];
138            CO2[index] = bigint(P[2]);
139            index++;
140            j++;
141         }
142      }
143
144      // insert deleteUnluckyPrimes
145      G = chinrem(CO1,CO2);
146      N = CO2[1];
147      for(i = 2; i <= size(CO2); i++){N = N*CO2[i];}
148      F = farey(G,N);
149
150      crit = 1;
151      for(i = 1; i <= nvars(basering); i++)
152      {
153         if(reduce(F[i],I) != 0) { crit = 0; break; }
154      }
155
156      if(!crit)
157      {
158         CO1 = G;
159         CO2 = N;
160         index = 2;
161
162         j = size(L) + 1;
163         L = primeList(I,10,L);
164
165         if(n > 1)
166         {
167            for(i = 1; i <= n; i++)
168            {
169               open(l(i));
170               write(l(i), quote(zeroRadP(eval(I), eval(L[j+i-1]))));
171            }
172            j = j + n;
173            k = 0;
174         }
175      }
176   }
177
178   k = 0;
179   for(i = 1; i <= nvars(basering); i++)
180   {
181      f = gcd(F[i],diff(F[i],var(i)));
182      k = k + deg(f);
183      F[i] = F[i]/f;
184   }
185
186   if(k == 0) { return(I); }
187   else { return(std(I + F)); }
188}
189example
190{ "EXAMPLE:";  echo = 2;
191   ring R = 0, (x,y), dp;
192   ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
193   zeroR(I);
194}
195
196////////////////////////////////////////////////////////////////////////////////
197
198proc assPrimes(list #)
199"USAGE:  assPrimes(I,[a],[n]); I ideal or module,
200         optional: n number of processors (for parallel computing), a
201         - a = 1: method of Eisenbud/Hunecke/Vasconcelos
202         - a = 2: method of Gianni/Trager/Zacharias
203         - a = 3: mathod of Monico
204         assPrimes(I) chooses n = a = 1
205ASSUME:  I is zero-dimensional over Q[variables]
206NOTE:    Parallelization is just applicable using 32-bit Singular version since
207         MP-links are not compatible with 64-bit Singular version.
208RETURN:  a list pr of associated primes of I:
209EXAMPLE: example assPrimes; shows an example
210"
211
212{
213   int T = timer;
214   int RT = rtimer;
215
216   ideal I;
217   if(typeof(#[1]) == "ideal")
218   {
219      I = #[1];
220   }
221   else
222   {
223      module M = #[1];
224      I = Ann(M);
225   }
226
227   //---------------------  Initialize optional parameter  ------------------------
228   if(size(#) > 1)
229   {
230      if(size(#) == 2)
231      {
232         int alg = #[2];
233         int n = 1;
234      }
235      if(size(#) == 3)
236      {
237         int alg = #[2];
238         int n = #[3];
239         if((n > 1) && (1 - system("with","MP")))
240         {
241            "========================================================================";
242            "There is no MP available on your system. Since this is necessary to     ";
243            "parallelize the algorithm, the computation will be done without forking.";
244            "========================================================================";
245            n = 1;
246         }
247      }
248   }
249   else
250   {
251      int alg = 1;
252      int n = 1;
253   }
254
255   def SPR = basering;
256   I = modStd(I,n);
257   int d = vdim(I);
258   if(d == -1) { ERROR("Ideal is not zero-dimensional."); }
259   if(homog(I) == 1){ return(list(maxideal(1))); }
260   poly f = findGen(I);
261   if(size(f) == nvars(SPR))
262   {
263      int spT = pTestRad(d,f,I);
264      if(printlevel>=10)
265      {
266         "pTestRad(d,f,I) = "+string(spT);
267      }
268      if(!spT)
269      {
270         if(typeof(attrib(#[1],"isRad")) == "int")
271         {
272            if(attrib(#[1],"isRad") == 0)
273            {
274               I = zeroR(I);
275               I = std(I);
276               d = vdim(I);
277               f = findGen(I);
278            }
279         }
280         else
281         {
282            I = zeroR(I);
283            I = std(I);
284            d = vdim(I);
285            f = findGen(I);
286         }
287      }
288   }
289   if(printlevel>=10)
290   {
291      "Real-time for radical-check is "+string(rtimer - RT)+" seconds.";
292      "CPU-time for radical-check is "+string(timer - T)+" seconds.";
293   }
294
295   export(SPR);
296   poly f_for_fork = f;
297   export(f_for_fork);           // f available for each link
298   ideal I_for_fork = I;
299   export(I_for_fork);           // I available for each link
300
301//--------------------  Initialize the list of primes  -------------------------
302   intvec L = primeList(I,10);
303   L[5] = prime(random(1000000000,2134567879));
304
305   list Re;
306
307   ring rHelp = 0,T,dp;
308   list CO1,CO2,P;
309   ideal F,G,testF;
310   bigint N;
311
312   list ringL = ringlist(SPR);
313   int i,k,e,int_break;
314   int j = 1;
315   int index = 1;
316
317//-----  If there is more than one processor available, we parallelize the  ----
318//-----  main standard basis computations in positive characteristic        ----
319
320   if(n > 1)
321   {
322
323//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
324//-----  standard basis for the primes L[2],...,L[n + 1].              ---------
325
326      for(i = 1; i <= n; i++)
327      {
328         link l(i) = "MPtcp:fork";
329         open(l(i));
330         write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[i + 1]), eval(alg))));
331      }
332
333      int t = timer;
334      P = modpSpecialAlgDep(ringL, L[1], alg);
335      t = timer - t;
336      if(t > 60) { t = 60; }
337      int i_sleep = system("sh", "sleep "+string(t));
338      CO1[index] = P[1];
339      CO2[index] = bigint(P[2]);
340      index++;
341
342      j = j + n + 1;
343   }
344
345//---------  Main computations in positive characteristic start here  ----------
346
347   int tt = timer;
348   int rt = rtimer;
349
350   while(1)
351   {
352      tt = timer;
353      rt = rtimer;
354
355      if(n > 1)
356      {
357         while(j <= size(L) + 1)
358         {
359            for(i = 1; i <= n; i++)
360            {
361               if(status(l(i), "read", "ready"))           // ask if link l(i) is ready otherwise sleep for t seconds
362               {
363                  P = read(l(i));                          // read the result from l(i)
364                  CO1[index] = P[1];
365                  CO2[index] = bigint(P[2]);
366                  index++;
367
368                  if(j <= size(L))
369                  {
370                     write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[j]), eval(alg))));
371                     j++;
372                  }
373                  else
374                  {
375                     k++;
376                     close(l(i));
377                  }
378               }
379            }
380            if(k == n)                                     // k describes the number of closed links
381            {
382               j++;
383            }
384            i_sleep = system("sh", "sleep "+string(t));
385         }
386      }
387      else
388      {
389         while(j<=size(L))
390         {
391            P = modpSpecialAlgDep(ringL, L[j], alg);
392            CO1[index] = P[1];
393            CO2[index] = bigint(P[2]);
394            index++;
395            j++;
396         }
397      }
398      if(printlevel>=10)
399      {
400         "Real-time for computing list in assPrimes is "+string(rtimer - rt)+" seconds.";
401         "CPU-time for computing list in assPrimes is "+string(timer - tt)+" seconds.";
402      }
403
404      // insert deleteUnluckyPrimes
405      G = chinrem(CO1,CO2);
406      N = CO2[1];
407      for(j = 2; j <= size(CO2); j++){N = N*CO2[j];}
408      F = farey(G,N);
409      if(F[1]-testF[1]==0)
410      {
411         if(printlevel>=10)
412         {
413            "size(L) = "+string(size(L));
414         }
415         F = cleardenom(F[1]);
416
417         e = deg(F[1]);
418         if(e == d)
419         {
420            list H = factorize(F[1]);
421
422            int s = size(H[1]);
423            for(i = 1; i <= s; i++)
424            {
425               if(H[2][i] != 1)
426               {
427                  int_break = 1;
428               }
429            }
430
431            if(int_break == 0)
432            {
433               setring SPR;
434               map phi = rHelp,var(nvars(SPR));
435               list H = phi(H);
436               if(printlevel>=10)
437               {
438                  "Real-time without test is "+string(rtimer - RT)+" seconds.";
439                  "CPU-time without test is "+string(timer - T)+" seconds.";
440               }
441               T = timer;
442               RT = rtimer;
443
444               ideal F = phi(F);
445               poly F1;
446
447               if(n > 1)
448               {
449                  open(l(1));
450                  write(l(1), quote(quickSubst(eval(F[1]), eval(f), eval(I))));
451                  int t_sleep = timer;
452               }
453               else
454               {
455                  F1 = quickSubst(F[1],f,I);
456                  if(F1 != 0) { int_break = 1; }
457               }
458
459               if(int_break == 0)
460               {
461                  for(i = 2; i <= s; i++)
462                  {
463                     H[1][i] = quickSubst(H[1][i],f,I);
464                     Re[i-1] = I + ideal(H[1][i]);
465                  }
466
467                  if(n > 1)
468                  {
469                     t_sleep = timer - t_sleep;
470                     if(t_sleep > 5) { t_sleep = 5; }
471
472                     while(1)
473                     {
474                        if(status(l(1), "read", "ready"))
475                        {
476                           F1 = read(l(1));
477                           close(l(1));
478                           break;
479                        }
480                        i_sleep = system("sh", "sleep "+string(t_sleep));
481                     }
482                     if(F1 != 0) { int_break = 1; }
483                  }
484                  if(printlevel>=10)
485                  {
486                     "Real-time for test is "+string(rtimer - RT)+" seconds.";
487                     "CPU-time for test is "+string(timer - T)+" seconds.";
488                  }
489                  if(int_break == 0)
490                  {
491                     kill f_for_fork;
492                     kill I_for_fork;
493                     kill SPR;
494                     return(Re);
495                  }
496               }
497            }
498         }
499      }
500
501      int_break = 0;
502      testF = F;
503      CO1 = G;
504      CO2 = N;
505      index = 2;
506
507      j = size(L) + 1;
508
509      setring(SPR);
510      L = primeList(I,10,L);
511      setring rHelp;
512
513      if(n > 1)
514      {
515         for(i = 1; i <= n; i++)
516         {
517            open(l(i));
518            write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[j+i-1]), eval(alg))));
519         }
520         j = j + n;
521         k = 0;
522      }
523   }
524}
525example
526{ "EXAMPLE:";  echo = 2;
527   ring R=0,(a,b,c,d,e,f),dp;
528   ideal I=
529   2fb+2ec+d2+a2+a,
530   2fc+2ed+2ba+b,
531   2fd+e2+2ca+c+b2,
532   2fe+2da+d+2cb,
533   f2+2ea+e+2db+c2,
534   2fa+f+2eb+2dc;
535   assPrimes(I);
536}
537////////////////////////////////////////////////////////////////////////////////
538
539static proc specialAlgDepEHV(poly p, ideal I)
540{
541//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
542//=== mapping T to p
543   def R = basering;
544   execute("ring Rhelp="+charstr(R)+",T,dp;");
545   setring R;
546   map phi = Rhelp,p;
547   setring Rhelp;
548   ideal F = preimage(R,phi,I); //corresponds to std(I,p-T) in dp(n),dp(1)
549   export(F);
550   setring R;
551   list L = Rhelp;
552   return(L);
553}
554
555////////////////////////////////////////////////////////////////////////////////
556
557static proc specialAlgDepGTZ(poly p, ideal I)
558{
559//=== assume I is zero-dimensional
560//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
561//=== mapping T to p
562   def R = basering;
563   execute("ring Rhelp = "+charstr(R)+",T,dp;");
564   setring R;
565   map phi = Rhelp,p;
566   def Rlp = changeord("dp("+string(nvars(R)-1)+"),dp(1)");
567   setring Rlp;
568   poly p = imap(R,p);
569   ideal K = maxideal(1);
570   K[nvars(R)] = 2*var(nvars(R))-p;
571   map phi = R,K;
572   ideal I = phi(I);
573   I = std(I);
574   poly q = subst(I[1],var(nvars(R)),var(1));
575   setring Rhelp;
576   map psi=Rlp,T;
577   ideal F=psi(q);
578   export(F);
579   setring R;
580   list L=Rhelp;
581   return(L);
582}
583
584////////////////////////////////////////////////////////////////////////////////
585
586static proc specialAlgDepMonico(poly p, ideal I)
587{
588//=== assume I is zero-dimensional
589//=== computes a poly F in Q[T], the characteristic polynomial of the map
590//=== basering/I ---> baserng/I  defined by the multiplication with p
591//=== in case I is radical it is the same poly as in specialAlgDepEHV
592   def R = basering;
593   execute("ring Rhelp = "+charstr(R)+",T,dp;");
594   setring R;
595   map phi = Rhelp,p;
596   poly q;
597   int j;
598   matrix m ;
599   poly va = var(1);
600   ideal J = std(I);
601   ideal ba = kbase(J);
602   int d = vdim(J);
603   matrix n[d][d];
604   for(j = 2; j <= nvars(R); j++)
605   {
606     va = va*var(j);
607   }
608   for(j = 1; j <= d; j++)
609   {
610     q = reduce(p*ba[j],J);
611     m = coeffs(q,ba,va);
612     n[j,1..d] = m[1..d,1];
613   }
614   setring Rhelp;
615   matrix n = imap(R,n);
616   ideal F = det(n-T*freemodule(d));
617   export(F);
618   setring R;
619   list L = Rhelp;
620   return(L);
621}
622
623////////////////////////////////////////////////////////////////////////////////
624
625static proc specialTest(int d, poly p, ideal I)
626{
627//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
628//=== mapping T to p and test if d=deg(F)
629   def R = basering;
630   execute("ring Rhelp="+charstr(R)+",T,dp;");
631   setring R;
632   map phi = Rhelp,p;
633   setring Rhelp;
634   ideal F = preimage(R,phi,I);
635   int e=deg(F[1]);
636   setring R;
637   return((e==d));
638}
639
640////////////////////////////////////////////////////////////////////////////////
641
642static proc findGen(ideal J, list #)
643{
644//=== try to find a sparse linear form r such that vector space dim(baserng/J)=deg(F),
645//=== F a poly in Q[T] such that <F>=kernel(Q[T]--->basering) mapping T to r
646//=== if not found returns a generic (randomly choosen) r
647   int d = vdim(J);
648   def R = basering;
649   int n = nvars(R);
650   list rl = ringlist(R);
651   if(size(#) > 0) { int p = #[1]; }
652   else { int p = prime(random(1000000000,2134567879)); }
653   rl[1] = p;
654   def @R = ring(rl);
655   setring @R;
656   ideal J = imap(R,J);
657   poly r = var(n);
658   int i,k;
659   k = specialTest(d,r,J);
660   if(!k)
661   {
662      for(i = 1; i < n; i++)
663      {
664         k = specialTest(d,r+var(i),J);
665         if(k){ r = r + var(i); break; }
666      }
667   }
668   for(i = 1; i < n; i++)
669   {
670      r = r + var(i);
671      k = specialTest(d,r,J);
672      if(k){break;}
673   }
674   setring R;
675   poly r = randomLast(100)[nvars(R)];
676   if(k){ r = imap(@R,r); }
677   return(r);
678}
679
680////////////////////////////////////////////////////////////////////////////////
681
682static proc pTestRad(int d, poly p1, ideal I)
683{
684//=== computes a poly F in Z/q1[T] such that <F>=kernel(Z/q1[T]--->Z/q1[vars(basering)])
685//=== mapping T to p1 and test if d=deg(squarefreepart(F)), q1 a prime randomly choosen
686//=== If not choose randomly another prime q2 and another linear form p2 and
687//=== computes a poly F in Z/q2[T] such that <F>=kernel(Z/q2[T]--->Z/q2[vars(basering)])
688//=== mapping T to p2 and test if d=deg(squarefreepart(F))
689//=== if the test is positive then I is radical
690   def R = basering;
691   list rl = ringlist(R);
692   int q1 = prime(random(100000000,536870912));
693   rl[1] = q1;
694   ring Shelp1 = q1,T,dp;
695   setring R;
696   def Rhelp1 = ring(rl);
697   setring Rhelp1;
698   poly p1 = imap(R,p1);
699   ideal I = imap(R,I);
700   map phi = Shelp1,p1;
701   setring Shelp1;
702   ideal F = preimage(Rhelp1,phi,I);
703   poly f = gcd(F[1],diff(F[1],var(1)));
704   int e = deg(F[1]/f);
705   setring R;
706   if(e != d)
707   {
708      poly p2 = findGen(I,q1);
709      setring Rhelp1;
710      poly p2 = imap(R,p2);
711      phi = Shelp1,p2;
712      setring Shelp1;
713      F = preimage(Rhelp1,phi,I);
714      f = gcd(F[1],diff(F[1],var(1)));
715      e = deg(F[1]/f);
716      setring R;
717      if(e == d){ return(1); }
718      if(e != d)
719      {
720         int q2 = prime(random(100000000,536870912));
721         rl[1] = q2;
722         ring Shelp2 = q2,T,dp;
723         setring R;
724         def Rhelp2 = ring(rl);
725         setring Rhelp2;
726         poly p1 = imap(R,p1);
727         ideal I = imap(R,I);
728         map phi = Shelp2,p1;
729         setring Shelp2;
730         ideal F = preimage(Rhelp2,phi,I);
731         poly f = gcd(F[1],diff(F[1],var(1)));
732         e = deg(F[1]/f);
733         setring R;
734         if(e == d){ return(1); }
735      }
736   }
737   return((e==d));
738}
739
740////////////////////////////////////////////////////////////////////////////////
741
742static proc zeroRadP(ideal I, int p)
743{
744//=== computes F=(F_1,...,F_n) such that <F_i>=IZ/p[x_1,...,x_n] intersected with Z/p[x_i], F_i monic
745   def R0 = basering;
746   list ringL = ringlist(R0);
747   ringL[1] = p;
748   def @r = ring(ringL);
749   setring @r;
750   ideal I = fetch(R0,I);
751   option(redSB);
752   I = std(I);
753   ideal F = finduni(I); //F[i] generates I intersected with K[var(i)]
754   int i;
755   for(i = 1; i <= size(F); i++){ F[i] = simplify(F[i],1); }
756   setring R0;
757   return(list(fetch(@r,F),p));
758}
759
760////////////////////////////////////////////////////////////////////////////////
761
762static proc quickSubst(poly h, poly r, ideal I)
763{
764//=== assume h is in Q[x_n], r is in Q[x_1,...,x_n], computes h(r) mod I
765   attrib(I,"isSB",1);
766   int n = nvars(basering);
767   poly q = 1;
768   int i,j,d;
769   intvec v;
770   list L;
771   for(i = 1; i <= size(h); i++)
772   {
773      L[i] = list(leadcoef(h[i]),leadexp(h[i])[n]);
774   }
775   d = L[1][2];
776   i = 0;
777   h = 0;
778
779   while(i <= d)
780   {
781      if(L[size(L)-j][2] == i)
782      {
783         h = reduce(h+L[size(L)-j][1]*q,I);
784         j++;
785      }
786      q = reduce(q*r,I);
787      i++;
788   }
789   return(h);
790}
791
792////////////////////////////////////////////////////////////////////////////////
793
794static proc modpSpecialAlgDep(list ringL, int p, list #)
795{
796//=== prepare parallel computing
797//===  #=1: method of Eisenbud/Hunecke/Vasconcelos
798//===  #=2: method of Gianni/Trager/Zacharias
799//===  #=3: mathod of Monico
800
801   def R0 = basering;
802
803   ringL[1] = p;
804   def @r = ring(ringL);
805   setring @r;
806   poly f = fetch(SPR,f_for_fork);
807   ideal I = fetch(SPR,I_for_fork);
808   if(size(#) > 0)
809   {
810      if(#[1] == 1) { list M = specialAlgDepEHV(f,I); }
811      if(#[1] == 2) { list M = specialAlgDepGTZ(f,I); }
812      if(#[1] == 3) { list M = specialAlgDepMonico(f,I); }
813   }
814   else
815   {
816      list M = specialAlgDepEHV(f,I);
817   }
818   def @S = M[1];
819
820   setring R0;
821   return(list(imap(@S,F),p));
822}
823
824////////////////////////////////////////////////////////////////////////////////
825
826/*
827Examples
828
829//Test for zeroR
830ring R = 0,(x,y),dp;
831ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
832
833//Cyclic 6
834ring R = 0,x(1..6),dp;
835ideal I = cyclic(6);
836
837//Amrhei
838ring R = 0,(a,b,c,d,e,f),dp;
839ideal I = 2fb+2ec+d2+a2+a,
840          2fc+2ed+2ba+b,
841          2fd+e2+2ca+c+b2,
842          2fe+2da+d+2cb,
843          f2+2ea+e+2db+c2,
844          2fa+f+2eb+2dc;
845
846
847//Becker-Niermann
848ring R = 0,(x,y,z),dp;
849ideal I = x2+xy2z-2xy+y4+y2+z2,
850          -x3y2+xy2z+xyz3-2xy+y4,
851          -2x2y+xy4+yz4-3;
852
853//Roczen
854ring R = 0,(a,b,c,d,e,f,g,h,k,o),dp;
855ideal I = o+1, k4+k, hk, h4+h, gk, gh, g3+h3+k3+1,
856          fk, f4+f, eh, ef, f3h3+e3k3+e3+f3+h3+k3+1,
857          e3g+f3g+g, e4+e, dh3+dk3+d, dg, df, de,
858          d3+e3+f3+1, e2g2+d2h2+c, f2g2+d2k2+b,
859          f2h2+e2k2+a;
860
861//4 bodies with equal masses, before symmetrisation.
862//We are looking for the real positive solutions
863ring R = 0,(B,b,D,d,F,f),dp;
864ideal I = (b-d)*(B-D)-2*F+2,
865          (b-d)*(B+D-2*F)+2*(B-D),
866          (b-d)^2-2*(b+d)+f+1,
867          B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
868
869ring R = 0,(B,b,D,d,F,f),dp;
870ideal I = (b-d)*(B-D)-F+1,
871          (b-d)*(B+D-F)+(B-D),
872          (b-d)^2-(b+d)+f+1,
873          B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
874
875
876//prime
877ring R = 0,(x,y,z,t,u),dp;
878ideal I = 2x2-2y2+2z2-2t2+2u2-1,
879          2x3-2y3+2z3-2t3+2u3-1,
880          2x4-2y4+2z4-2t4+2u4-1,
881          2x5-2y5+2z5-2t5+2u5-1,
882          2x6-2y6+2z6-2t6+2u6-1;
883
884*/
885
Note: See TracBrowser for help on using the repository browser.