source: git/Singular/LIB/assprimeszerodim.lib @ d008518

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