source: git/Singular/LIB/assprimeszerodim.lib @ 0c2f6d

fieker-DuValspielwiese
Last change on this file since 0c2f6d was 8f1d129, checked in by Stefan Steidel <steidel@…>, 13 years ago
Changes in procedures zeroR resp. findGen and attachment of further examples git-svn-id: file:///usr/local/Singular/svn/trunk@13777 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 27.7 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 zeroRadical(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 zeroRadical(ideal I, list #)
30"USAGE:  zeroRadical(I,[n]); I ideal, optional: n number of processors (for
31         parallel computing)
32ASSUME:  I is zero-dimensional in Q[variables]
33NOTE:    Parallelization is just applicable using 32-bit Singular version since
34         MP-links are not compatible with 64-bit Singular version.
35RETURN:  the radical of I
36EXAMPLE: example zeroRadical; shows an example
37"
38{
39   return(zeroR(modStd(I,#),#));
40}
41example
42{ "EXAMPLE:";  echo = 2;
43   ring R = 0, (x,y), dp;
44   ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
45   zeroRadical(I);
46}
47
48////////////////////////////////////////////////////////////////////////////////
49
50static proc zeroR(ideal I, list #)
51// compute the radical of I provided that I is zero-dimensional in Q[variables]
52// and a standard basis
53{
54   attrib(I,"isSB",1);
55   int i, k;
56   int j = 1;
57   int index = 1;
58   int crit;
59
60   list CO1, CO2, P;
61   ideal G, F;
62   bigint N;
63   poly f;
64
65//---------------------  Initialize optional parameter  ------------------------
66   if(size(#) > 0)
67   {
68      int n = #[1];
69      if((n > 1) && (1 - system("with","MP")))
70      {
71     "========================================================================";
72     "There is no MP available on your system. Since this is necessary to     ";
73     "parallelize the algorithm, the computation will be done without forking.";
74     "========================================================================";
75         n = 1;
76      }
77   }
78   else
79   {
80      int n = 1;
81   }
82
83//--------------------  Initialize the list of primes  -------------------------
84   intvec L = primeList(I,10);
85   L[5] = prime(random(100000000,536870912));
86
87   if(n > 1)
88   {
89
90//-----  Create n links l(1),...,l(n), open all of them and compute  -----------
91//-----  polynomial F for the primes L[2],...,L[n + 1].              -----------
92
93      for(i = 1; i <= n; i++)
94      {
95         link l(i) = "MPtcp:fork";
96         open(l(i));
97         write(l(i), quote(zeroRadP(eval(I), eval(L[i + 1]))));
98      }
99
100      int t = timer;
101      P = zeroRadP(I, L[1]);
102      t = timer - t;
103      if(t > 60) { t = 60; }
104      int i_sleep = system("sh", "sleep "+string(t));
105      CO1[index] = P[1];
106      CO2[index] = bigint(P[2]);
107      index++;
108
109      j = j + n + 1;
110   }
111
112//---------  Main computations in positive characteristic start here  ----------
113
114   while(!crit)
115   {
116      if(n > 1)
117      {
118         while(j <= size(L) + 1)
119         {
120            for(i = 1; i <= n; i++)
121            {
122               if(status(l(i), "read", "ready"))
123               {
124                  //--- read the result from l(i) ---
125                  P = read(l(i));
126                  CO1[index] = P[1];
127                  CO2[index] = bigint(P[2]);
128                  index++;
129
130                  if(j <= size(L))
131                  {
132                     write(l(i), quote(zeroRadP(eval(I), eval(L[j]))));
133                     j++;
134                  }
135                  else
136                  {
137                     k++;
138                     close(l(i));
139                  }
140               }
141            }
142            //--- k describes the number of closed links ---
143            if(k == n)
144            {
145               j++;
146            }
147            //--- sleep for t seconds ---
148            i_sleep = system("sh", "sleep "+string(t));
149         }
150      }
151      else
152      {
153         while(j<=size(L))
154         {
155            P = zeroRadP(I, L[j]);
156            CO1[index] = P[1];
157            CO2[index] = bigint(P[2]);
158            index++;
159            j++;
160         }
161      }
162
163      // insert deleteUnluckyPrimes
164      G = chinrem(CO1,CO2);
165      N = CO2[1];
166      for(i = 2; i <= size(CO2); i++){N = N*CO2[i];}
167      F = farey(G,N);
168
169      crit = 1;
170      for(i = 1; i <= nvars(basering); i++)
171      {
172         if(reduce(F[i],I) != 0) { crit = 0; break; }
173      }
174
175      if(!crit)
176      {
177         CO1 = G;
178         CO2 = N;
179         index = 2;
180
181         j = size(L) + 1;
182         L = primeList(I,10,L);
183
184         if(n > 1)
185         {
186            for(i = 1; i <= n; i++)
187            {
188               open(l(i));
189               write(l(i), quote(zeroRadP(eval(I), eval(L[j+i-1]))));
190            }
191            j = j + n;
192            k = 0;
193         }
194      }
195   }
196
197   k = 0;
198   for(i = 1; i <= nvars(basering); i++)
199   {
200      f = gcd(F[i],diff(F[i],var(i)));
201      k = k + deg(f);
202      F[i] = F[i]/f;
203   }
204
205   if(k == 0) { return(I); }
206   else { return(modStd(I + F)); }
207}
208
209////////////////////////////////////////////////////////////////////////////////
210
211proc assPrimes(list #)
212"USAGE:  assPrimes(I,[a],[n]); I ideal or module,
213         optional: int n: number of processors (for parallel computing), int a:
214         - a = 1: method of Eisenbud/Hunecke/Vasconcelos
215         - a = 2: method of Gianni/Trager/Zacharias
216         - a = 3: method of Monico
217         assPrimes(I) chooses n = a = 1
218ASSUME:  I is zero-dimensional over Q[variables]
219NOTE:    Parallelization is just applicable using 32-bit Singular version since
220         MP-links are not compatible with 64-bit Singular version.
221RETURN:  a list Re of associated primes of I:
222EXAMPLE: example assPrimes; shows an example
223"
224{
225   ideal I;
226   if(typeof(#[1]) == "ideal")
227   {
228      I = #[1];
229   }
230   else
231   {
232      module M = #[1];
233      I = Ann(M);
234   }
235
236//---------------------  Initialize optional parameter  ------------------------
237   if(size(#) > 1)
238   {
239      if(size(#) == 2)
240      {
241         int alg = #[2];
242         int n = 1;
243      }
244      if(size(#) == 3)
245      {
246         int alg = #[2];
247         int n = #[3];
248         if((n > 1) && (1 - system("with","MP")))
249         {
250     "========================================================================";
251     "There is no MP available on your system. Since this is necessary to     ";
252     "parallelize the algorithm, the computation will be done without forking.";
253     "========================================================================";
254            n = 1;
255         }
256      }
257   }
258   else
259   {
260      int alg = 1;
261      int n = 1;
262   }
263
264   def SPR = basering;
265   I = modStd(I,n);
266   int d = vdim(I);
267   if(d == -1) { ERROR("Ideal is not zero-dimensional."); }
268   if(homog(I) == 1){ return(list(maxideal(1))); }
269   poly f = findGen(I);
270
271   int T = timer;
272   int RT = rtimer;
273   int TT;
274   if(size(f) == nvars(SPR))
275   {
276      TT = timer;
277      int spT = pTestRad(d,f,I);
278      if(printlevel>=10)
279      {
280         "pTestRad(d,f,I) = "+string(spT)+" and takes
281         "+string(timer-TT)+" seconds.";
282      }
283      if(!spT)
284      {
285         if(typeof(attrib(#[1],"isRad")) == "int")
286         {
287            if(attrib(#[1],"isRad") == 0)
288            {
289               TT = timer;
290               I = zeroR(I,n);
291               if(printlevel>=10)
292               {
293                  "zeroR(I,n) takes "+string(timer-TT)+" seconds.";
294               }
295               TT = timer;
296               I = std(I);
297               if(printlevel>=10)
298               {
299                  "std(I) takes "+string(timer-TT)+" seconds.";
300               }
301               d = vdim(I);
302               f = findGen(I);
303            }
304         }
305         else
306         {
307            TT = timer;
308            I = zeroR(I,n);
309            if(printlevel>=10)
310            {
311               "zeroR(I,n) takes "+string(timer-TT)+" seconds.";
312            }
313            TT = timer;
314            I = std(I);
315            if(printlevel>=10)
316            {
317               "std(I) takes "+string(timer-TT)+" seconds.";
318            }
319            d = vdim(I);
320            f = findGen(I);
321         }
322      }
323   }
324   if(printlevel>=10)
325   {
326      "Real-time for radical-check is "+string(rtimer - RT)+" seconds.";
327      "CPU-time for radical-check is "+string(timer - T)+" seconds.";
328   }
329
330   export(SPR);
331   poly f_for_fork = f;
332   export(f_for_fork);           // f available for each link
333   ideal I_for_fork = I;
334   export(I_for_fork);           // I available for each link
335
336//--------------------  Initialize the list of primes  -------------------------
337   intvec L = primeList(I,10);
338   L[5] = prime(random(1000000000,2134567879));
339
340   list Re;
341
342   ring rHelp = 0,T,dp;
343   list CO1,CO2,P;
344   ideal F,G,testF;
345   bigint N;
346
347   list ringL = ringlist(SPR);
348   int i,k,e,int_break;
349   int j = 1;
350   int index = 1;
351
352//-----  If there is more than one processor available, we parallelize the  ----
353//-----  main standard basis computations in positive characteristic        ----
354
355   if(n > 1)
356   {
357
358//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
359//-----  standard basis for the primes L[2],...,L[n + 1].              ---------
360
361      for(i = 1; i <= n; i++)
362      {
363         link l(i) = "MPtcp:fork";
364         open(l(i));
365         write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[i + 1]),
366                                                          eval(alg))));
367      }
368
369      int t = timer;
370      P = modpSpecialAlgDep(ringL, L[1], alg);
371      t = timer - t;
372      if(t > 60) { t = 60; }
373      int i_sleep = system("sh", "sleep "+string(t));
374      CO1[index] = P[1];
375      CO2[index] = bigint(P[2]);
376      index++;
377
378      j = j + n + 1;
379   }
380
381//---------  Main computations in positive characteristic start here  ----------
382
383   int tt = timer;
384   int rt = rtimer;
385
386   while(1)
387   {
388      tt = timer;
389      rt = rtimer;
390
391      if(n > 1)
392      {
393         while(j <= size(L) + 1)
394         {
395            for(i = 1; i <= n; i++)
396            {
397               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
398               if(status(l(i), "read", "ready"))
399               {
400                  //--- read the result from l(i) ---
401                  P = read(l(i));
402                  CO1[index] = P[1];
403                  CO2[index] = bigint(P[2]);
404                  index++;
405
406                  if(j <= size(L))
407                  {
408                     write(l(i), quote(modpSpecialAlgDep(eval(ringL),
409                                                         eval(L[j]),
410                                                         eval(alg))));
411                     j++;
412                  }
413                  else
414                  {
415                     k++;
416                     close(l(i));
417                  }
418               }
419            }
420            //--- k describes the number of closed links ---
421            if(k == n)
422            {
423               j++;
424            }
425            i_sleep = system("sh", "sleep "+string(t));
426         }
427      }
428      else
429      {
430         while(j<=size(L))
431         {
432            P = modpSpecialAlgDep(ringL, L[j], alg);
433            CO1[index] = P[1];
434            CO2[index] = bigint(P[2]);
435            index++;
436            j++;
437         }
438      }
439      if(printlevel>=10)
440      {
441         "Real-time for computing list in assPrimes is "+string(rtimer - rt)+
442         " seconds.";
443         "CPU-time for computing list in assPrimes is "+string(timer - tt)+
444         " seconds.";
445      }
446
447      // insert deleteUnluckyPrimes
448      G = chinrem(CO1,CO2);
449      N = CO2[1];
450      for(j = 2; j <= size(CO2); j++){N = N*CO2[j];}
451      F = farey(G,N);
452      if(F[1]-testF[1]==0)
453      {
454         if(printlevel>=10)
455         {
456            "size(L) = "+string(size(L));
457         }
458         F = cleardenom(F[1]);
459
460         e = deg(F[1]);
461         if(e == d)
462         {
463            list H = factorize(F[1]);
464
465            int s = size(H[1]);
466            for(i = 1; i <= s; i++)
467            {
468               if(H[2][i] != 1)
469               {
470                  int_break = 1;
471               }
472            }
473
474            if(int_break == 0)
475            {
476               setring SPR;
477               map phi = rHelp,var(nvars(SPR));
478               list H = phi(H);
479               if(printlevel>=10)
480               {
481                  "Real-time without test is "+string(rtimer - RT)+" seconds.";
482                  "CPU-time without test is "+string(timer - T)+" seconds.";
483               }
484               T = timer;
485               RT = rtimer;
486
487               ideal F = phi(F);
488               poly F1;
489
490               if(n > 1)
491               {
492                  open(l(1));
493                  write(l(1), quote(quickSubst(eval(F[1]), eval(f), eval(I))));
494                  int t_sleep = timer;
495               }
496               else
497               {
498                  F1 = quickSubst(F[1],f,I);
499                  if(F1 != 0) { int_break = 1; }
500               }
501
502               if(int_break == 0)
503               {
504                  for(i = 2; i <= s; i++)
505                  {
506                     H[1][i] = quickSubst(H[1][i],f,I);
507                     Re[i-1] = I + ideal(H[1][i]);
508                  }
509
510                  if(n > 1)
511                  {
512                     t_sleep = timer - t_sleep;
513                     if(t_sleep > 5) { t_sleep = 5; }
514
515                     while(1)
516                     {
517                        if(status(l(1), "read", "ready"))
518                        {
519                           F1 = read(l(1));
520                           close(l(1));
521                           break;
522                        }
523                        i_sleep = system("sh", "sleep "+string(t_sleep));
524                     }
525                     if(F1 != 0) { int_break = 1; }
526                  }
527                  if(printlevel>=10)
528                  {
529                     "Real-time for test is "+string(rtimer - RT)+" seconds.";
530                     "CPU-time for test is "+string(timer - T)+" seconds.";
531                  }
532                  if(int_break == 0)
533                  {
534                     kill f_for_fork;
535                     kill I_for_fork;
536                     kill SPR;
537                     return(Re);
538                  }
539               }
540            }
541         }
542      }
543
544      int_break = 0;
545      testF = F;
546      CO1 = G;
547      CO2 = N;
548      index = 2;
549
550      j = size(L) + 1;
551
552      setring(SPR);
553      L = primeList(I,10,L);
554      setring rHelp;
555
556      if(n > 1)
557      {
558         for(i = 1; i <= n; i++)
559         {
560            open(l(i));
561            write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[j+i-1]),
562                                                             eval(alg))));
563         }
564         j = j + n;
565         k = 0;
566      }
567   }
568}
569example
570{ "EXAMPLE:";  echo = 2;
571   ring R=0,(a,b,c,d,e,f),dp;
572   ideal I=
573   2fb+2ec+d2+a2+a,
574   2fc+2ed+2ba+b,
575   2fd+e2+2ca+c+b2,
576   2fe+2da+d+2cb,
577   f2+2ea+e+2db+c2,
578   2fa+f+2eb+2dc;
579   assPrimes(I);
580}
581////////////////////////////////////////////////////////////////////////////////
582
583static proc specialAlgDepEHV(poly p, ideal I)
584{
585//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
586//=== mapping T to p
587   def R = basering;
588   execute("ring Rhelp="+charstr(R)+",T,dp;");
589   setring R;
590   map phi = Rhelp,p;
591   setring Rhelp;
592   ideal F = preimage(R,phi,I); //corresponds to std(I,p-T) in dp(n),dp(1)
593   export(F);
594   setring R;
595   list L = Rhelp;
596   return(L);
597}
598
599////////////////////////////////////////////////////////////////////////////////
600
601static proc specialAlgDepGTZ(poly p, ideal I)
602{
603//=== assume I is zero-dimensional
604//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
605//=== mapping T to p
606   def R = basering;
607   execute("ring Rhelp = "+charstr(R)+",T,dp;");
608   setring R;
609   map phi = Rhelp,p;
610   def Rlp = changeord("dp("+string(nvars(R)-1)+"),dp(1)");
611   setring Rlp;
612   poly p = imap(R,p);
613   ideal K = maxideal(1);
614   K[nvars(R)] = 2*var(nvars(R))-p;
615   map phi = R,K;
616   ideal I = phi(I);
617   I = std(I);
618   poly q = subst(I[1],var(nvars(R)),var(1));
619   setring Rhelp;
620   map psi=Rlp,T;
621   ideal F=psi(q);
622   export(F);
623   setring R;
624   list L=Rhelp;
625   return(L);
626}
627
628////////////////////////////////////////////////////////////////////////////////
629
630static proc specialAlgDepMonico(poly p, ideal I)
631{
632//=== assume I is zero-dimensional
633//=== computes a poly F in Q[T], the characteristic polynomial of the map
634//=== basering/I ---> baserng/I  defined by the multiplication with p
635//=== in case I is radical it is the same poly as in specialAlgDepEHV
636   def R = basering;
637   execute("ring Rhelp = "+charstr(R)+",T,dp;");
638   setring R;
639   map phi = Rhelp,p;
640   poly q;
641   int j;
642   matrix m ;
643   poly va = var(1);
644   ideal J = std(I);
645   ideal ba = kbase(J);
646   int d = vdim(J);
647   matrix n[d][d];
648   for(j = 2; j <= nvars(R); j++)
649   {
650     va = va*var(j);
651   }
652   for(j = 1; j <= d; j++)
653   {
654     q = reduce(p*ba[j],J);
655     m = coeffs(q,ba,va);
656     n[j,1..d] = m[1..d,1];
657   }
658   setring Rhelp;
659   matrix n = imap(R,n);
660   ideal F = det(n-T*freemodule(d));
661   export(F);
662   setring R;
663   list L = Rhelp;
664   return(L);
665}
666
667////////////////////////////////////////////////////////////////////////////////
668
669static proc specialTest(int d, poly p, ideal I)
670{
671//=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
672//=== mapping T to p and test if d=deg(F)
673   def R = basering;
674   execute("ring Rhelp="+charstr(R)+",T,dp;");
675   setring R;
676   map phi = Rhelp,p;
677   setring Rhelp;
678   ideal F = preimage(R,phi,I);
679   int e=deg(F[1]);
680   setring R;
681   return((e==d));
682}
683
684////////////////////////////////////////////////////////////////////////////////
685
686static proc findGen(ideal J, list #)
687{
688//=== try to find a sparse linear form r such that
689//=== vector space dim(basering/J)=deg(F),
690//=== F a poly in Q[T] such that <F>=kernel(Q[T]--->basering) mapping T to r
691//=== if not found returns a generic (randomly choosen) r
692   int d = vdim(J);
693   def R = basering;
694   int n = nvars(R);
695   list rl = ringlist(R);
696   if(size(#) > 0) { int p = #[1]; }
697   else { int p = prime(random(1000000000,2134567879)); }
698   rl[1] = p;
699   def @R = ring(rl);
700   setring @R;
701   ideal J = imap(R,J);
702   poly r = var(n);
703   int i,k;
704   k = specialTest(d,r,J);
705   if(!k)
706   {
707      for(i = 1; i < n; i++)
708      {
709         k = specialTest(d,r+var(i),J);
710         if(k){ r = r + var(i); break; }
711      }
712   }
713   if((!k) && (n > 2))
714   {
715      for(i = 1; i < n; i++)
716      {
717         r = r + var(i);
718         k = specialTest(d,r,J);
719         if(k){ break; }
720      }
721   }
722   setring R;
723   poly r = randomLast(100)[nvars(R)];
724   if(k){ r = imap(@R,r); }
725   return(r);
726}
727
728////////////////////////////////////////////////////////////////////////////////
729
730static proc pTestRad(int d, poly p1, ideal I)
731{
732//=== computes a poly F in Z/q1[T] such that
733//===                    <F> = kernel(Z/q1[T]--->Z/q1[vars(basering)])
734//=== mapping T to p1 and test if d=deg(squarefreepart(F)), q1 a prime randomly
735//=== chosen
736//=== If not choose randomly another prime q2 and another linear form p2 and
737//=== computes a poly F in Z/q2[T] such that
738//===                    <F> = kernel(Z/q2[T]--->Z/q2[vars(basering)])
739//=== mapping T to p2 and test if d=deg(squarefreepart(F))
740//=== if the test is positive then I is radical
741   def R = basering;
742   list rl = ringlist(R);
743   int q1 = prime(random(100000000,536870912));
744   rl[1] = q1;
745   ring Shelp1 = q1,T,dp;
746   setring R;
747   def Rhelp1 = ring(rl);
748   setring Rhelp1;
749   poly p1 = imap(R,p1);
750   ideal I = imap(R,I);
751   map phi = Shelp1,p1;
752   setring Shelp1;
753   ideal F = preimage(Rhelp1,phi,I);
754   poly f = gcd(F[1],diff(F[1],var(1)));
755   int e = deg(F[1]/f);
756   setring R;
757   if(e != d)
758   {
759      poly p2 = findGen(I,q1);
760      setring Rhelp1;
761      poly p2 = imap(R,p2);
762      phi = Shelp1,p2;
763      setring Shelp1;
764      F = preimage(Rhelp1,phi,I);
765      f = gcd(F[1],diff(F[1],var(1)));
766      e = deg(F[1]/f);
767      setring R;
768      if(e == d){ return(1); }
769      if(e != d)
770      {
771         int q2 = prime(random(100000000,536870912));
772         rl[1] = q2;
773         ring Shelp2 = q2,T,dp;
774         setring R;
775         def Rhelp2 = ring(rl);
776         setring Rhelp2;
777         poly p1 = imap(R,p1);
778         ideal I = imap(R,I);
779         map phi = Shelp2,p1;
780         setring Shelp2;
781         ideal F = preimage(Rhelp2,phi,I);
782         poly f = gcd(F[1],diff(F[1],var(1)));
783         e = deg(F[1]/f);
784         setring R;
785         if(e == d){ return(1); }
786      }
787   }
788   return((e==d));
789}
790
791////////////////////////////////////////////////////////////////////////////////
792
793static proc zeroRadP(ideal I, int p)
794{
795//=== computes F=(F_1,...,F_n) such that <F_i>=IZ/p[x_1,...,x_n] intersected
796//=== with Z/p[x_i], F_i monic
797   def R0 = basering;
798   list ringL = ringlist(R0);
799   ringL[1] = p;
800   def @r = ring(ringL);
801   setring @r;
802   ideal I = fetch(R0,I);
803   option(redSB);
804   I = std(I);
805   ideal F = finduni(I); //F[i] generates I intersected with K[var(i)]
806   int i;
807   for(i = 1; i <= size(F); i++){ F[i] = simplify(F[i],1); }
808   setring R0;
809   return(list(fetch(@r,F),p));
810}
811
812////////////////////////////////////////////////////////////////////////////////
813
814static proc quickSubst(poly h, poly r, ideal I)
815{
816//=== assume h is in Q[x_n], r is in Q[x_1,...,x_n], computes h(r) mod I
817   attrib(I,"isSB",1);
818   int n = nvars(basering);
819   poly q = 1;
820   int i,j,d;
821   intvec v;
822   list L;
823   for(i = 1; i <= size(h); i++)
824   {
825      L[i] = list(leadcoef(h[i]),leadexp(h[i])[n]);
826   }
827   d = L[1][2];
828   i = 0;
829   h = 0;
830
831   while(i <= d)
832   {
833      if(L[size(L)-j][2] == i)
834      {
835         h = reduce(h+L[size(L)-j][1]*q,I);
836         j++;
837      }
838      q = reduce(q*r,I);
839      i++;
840   }
841   return(h);
842}
843
844////////////////////////////////////////////////////////////////////////////////
845
846static proc modpSpecialAlgDep(list ringL, int p, list #)
847{
848//=== prepare parallel computing
849//===  #=1: method of Eisenbud/Hunecke/Vasconcelos
850//===  #=2: method of Gianni/Trager/Zacharias
851//===  #=3: method of Monico
852
853   def R0 = basering;
854
855   ringL[1] = p;
856   def @r = ring(ringL);
857   setring @r;
858   poly f = fetch(SPR,f_for_fork);
859   ideal I = fetch(SPR,I_for_fork);
860   if(size(#) > 0)
861   {
862      if(#[1] == 1) { list M = specialAlgDepEHV(f,I); }
863      if(#[1] == 2) { list M = specialAlgDepGTZ(f,I); }
864      if(#[1] == 3) { list M = specialAlgDepMonico(f,I); }
865   }
866   else
867   {
868      list M = specialAlgDepEHV(f,I);
869   }
870   def @S = M[1];
871
872   setring R0;
873   return(list(imap(@S,F),p));
874}
875
876////////////////////////////////////////////////////////////////////////////////
877
878/*
879Examples:
880=========
881
882//=== Test for zeroR
883ring R = 0,(x,y),dp;
884ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
885
886//=== Cyclic_6
887ring R = 0,x(1..6),dp;
888ideal I = cyclic(6);
889
890//=== Amrhein
891ring R = 0,(a,b,c,d,e,f),dp;
892ideal I = 2fb+2ec+d2+a2+a,
893          2fc+2ed+2ba+b,
894          2fd+e2+2ca+c+b2,
895          2fe+2da+d+2cb,
896          f2+2ea+e+2db+c2,
897          2fa+f+2eb+2dc;
898
899//=== Becker-Niermann
900ring R = 0,(x,y,z),dp;
901ideal I = x2+xy2z-2xy+y4+y2+z2,
902          -x3y2+xy2z+xyz3-2xy+y4,
903          -2x2y+xy4+yz4-3;
904
905//=== Roczen
906ring R = 0,(a,b,c,d,e,f,g,h,k,o),dp;
907ideal I = o+1, k4+k, hk, h4+h, gk, gh, g3+h3+k3+1,
908          fk, f4+f, eh, ef, f3h3+e3k3+e3+f3+h3+k3+1,
909          e3g+f3g+g, e4+e, dh3+dk3+d, dg, df, de,
910          d3+e3+f3+1, e2g2+d2h2+c, f2g2+d2k2+b,
911          f2h2+e2k2+a;
912
913//=== FourBodyProblem
914//=== 4 bodies with equal masses, before symmetrisation.
915//=== We are looking for the real positive solutions
916ring R = 0,(B,b,D,d,F,f),dp;
917ideal I = (b-d)*(B-D)-2*F+2,
918          (b-d)*(B+D-2*F)+2*(B-D),
919          (b-d)^2-2*(b+d)+f+1,
920          B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
921
922//=== Reimer_5
923ring R = 0,(x,y,z,t,u),dp;
924ideal I = 2x2-2y2+2z2-2t2+2u2-1,
925          2x3-2y3+2z3-2t3+2u3-1,
926          2x4-2y4+2z4-2t4+2u4-1,
927          2x5-2y5+2z5-2t5+2u5-1,
928          2x6-2y6+2z6-2t6+2u6-1;
929
930//=== ZeroDim.example_12
931ring R = 0, (x, y, z), lp;
932ideal I = 7xy+x+3yz+4y+2z+10,
933          x3+x2y+3xyz+5xy+3x+2y3+6y2z+yz+1,
934          3x4+2x2y2+3x2y+4x2z2+xyz+xz2+6y2z+5z4;
935
936//=== ZeroDim.example_27
937ring R = 0, (w, x, y, z), lp;
938ideal I = -2w2+9wx+9wy-7wz-4w+8x2+9xy-3xz+8x+6y2-7yz+4y-6z2+8z+2,
939          3w2-5wx-3wy-6wz+9w+4x2+2xy-2xz+7x+9y2+6yz+5y+7z2+7z+5,
940          7w2+5wx+3wy-5wz-5w+2x2+9xy-7xz+4x-4y2-5yz+6y-4z2-9z+2,
941          8w2+5wx-4wy+2wz+3w+5x2+2xy-7xz-7x+7y2-8yz-7y+7z2-8z+8;
942
943//=== Cassou_1
944ring R = 0, (b,c,d,e), dp;
945ideal I = 6b4c3+21b4c2d+15b4cd2+9b4d3+36b4c2+84b4cd+30b4d2+72b4c+84b4d-8b2c2e
946          -28b2cde+36b2d2e+48b4-32b2ce-56b2de-144b2c-648b2d-32b2e-288b2-120,
947          9b4c4+30b4c3d+39b4c2d2+18b4cd3+72b4c3+180b4c2d+156b4cd2+36b4d3+216b4c2
948          +360b4cd+156b4d2-24b2c3e-16b2c2de+16b2cd2e+24b2d3e+288b4c+240b4d
949          -144b2c2e+32b2d2e+144b4-432b2c2-720b2cd-432b2d2-288b2ce+16c2e2-32cde2
950          +16d2e2-1728b2c-1440b2d-192b2e+64ce2-64de2-1728b2+576ce-576de+64e2
951          -240c+1152e+4704,
952          -15b2c3e+15b2c2de-90b2c2e+60b2cde-81b2c2+216b2cd-162b2d2-180b2ce
953          +60b2de+40c2e2-80cde2+40d2e2-324b2c+432b2d-120b2e+160ce2-160de2-324b2
954          +1008ce-1008de+160e2+2016e+5184,
955          -4b2c2+4b2cd-3b2d2-16b2c+8b2d-16b2+22ce-22de+44e+261;
956
957================================================================================
958
959The following timings are conducted on an Intel Xeon X5460 with 4 CPUs, 3.16 GHz
960each, 64 GB RAM under the Gentoo Linux operating system by using the 32-bit
961version of Singular 3-1-1.
962The results of the timinings are summarized in the following table where
963assPrimes* denotes the parallelized version of the algorithm on 4 CPUs and
964 (1) approach of Eisenbud, Hunecke, Vasconcelos (cf. specialAlgDepEHV),
965 (2) approach of Gianni, Trager, Zacharias      (cf. specialAlgDepGTZ),
966 (3) approach of Monico                         (cf. specialAlgDepMonico).
967
968Example             | minAssGTZ |     assPrimes         assPrimes*
969                    |           |  (1)   (2)   (3)   (1)   (2)   (3)
970--------------------------------------------------------------------
971Cyclic_6            |      5    |    5    10     7     4     7     6
972Amrhein             |      1    |    3     3     5     1     2    21
973Becker-Niermann     |      -    |    0     0     1     0     0     0
974Roczen              |      0    |    3     2     0     2     4     1
975FourBodyProblem     |      -    |  139   139   148    96    83    96
976Reimer_5            |      -    |  132   128   175    97    70   103
977ZeroDim.example_12  |    170    |  125   125   125    67    68    63
978ZeroDim.example_27  |     27    |  215   226   215   113   117   108
979Cassou_1            |    525    |  112   112   112    56    56    57
980
981minAssGTZ (cf. primdec_lib) runs out of memory for Becker-Niermann,
982FourBodyProblem and Reimer_5.
983
984================================================================================
985
986//=== One component at the origin
987
988ring R = 0, (x,y), dp;
989poly f1 = (y5 + y4x7 + 2x8);
990poly f2 = (y3 + 7x4);
991poly f3 = (y7 + 2x12);
992poly f = f1*f2*f3 + y19;
993ideal I = f, diff(f, x), diff(f, y);
994
995ring R = 0, (x,y), dp;
996poly f1 = (y5 + y4x7 + 2x8);
997poly f2 = (y3 + 7x4);
998poly f3 = (y7 + 2x12);
999poly f4 = (y11 + 2x18);
1000poly f = f1*f2*f3*f4 + y30;
1001ideal I = f, diff(f, x), diff(f, y);
1002
1003ring R = 0, (x,y), dp;
1004poly f1 = (y15 + y14x7 + 2x18);
1005poly f2 = (y13 + 7x14);
1006poly f3 = (y17 + 2x22);
1007poly f = f1*f2*f3 + y49;
1008ideal I = f, diff(f, x), diff(f, y);
1009
1010ring R = 0, (x,y), dp;
1011poly f1 = (y15 + y14x20 + 2x38);
1012poly f2 = (y19 + 3y17x50 + 7x52);
1013poly f = f1*f2 + y36;
1014ideal I = f, diff(f, x), diff(f, y);
1015
1016//=== Several components
1017
1018ring R = 0, (x,y), dp;
1019poly f1 = (y5 + y4x7 + 2x8);
1020poly f2 = (y13 + 7x14);
1021poly f = f1*subst(f2, x, x-3, y, y+5);
1022ideal I = f, diff(f, x), diff(f, y);
1023
1024ring R = 0, (x,y), dp;
1025poly f1 = (y5  + y4x7 + 2x8);
1026poly f2 = (y3 + 7x4);
1027poly f3 = (y7 + 2x12);
1028poly f = f1*f2*subst(f3, y, y+5);
1029ideal I = f, diff(f, x), diff(f, y);
1030
1031ring R = 0, (x,y), dp;
1032poly f1 = (y5  + 2x8);
1033poly f2 = (y3 + 7x4);
1034poly f3 = (y7 + 2x12);
1035poly f = f1*subst(f2,x,x-1)*subst(f3, y, y+5);
1036ideal I = f, diff(f, x), diff(f, y);
1037*/
1038
Note: See TracBrowser for help on using the repository browser.