source: git/Singular/LIB/assprime.lib @ a15000c

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