source: git/Singular/LIB/modstd.lib @ 8d1432e

spielwiese
Last change on this file since 8d1432e was 9761b8, checked in by Andreas Steenpass <steenpass@…>, 9 years ago
chg: use alias in Modstd::pTest_std()
  • Property mode set to 100644
File size: 18.9 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version modstd.lib 4.0.2.0 Aug_2015 "; // $Id$
3category="Commutative Algebra";
4info="
5LIBRARY:  modstd.lib      Groebner bases of ideals using modular methods
6
7AUTHORS:  A. Hashemi      Amir.Hashemi@lip6.fr
8          G. Pfister      pfister@mathematik.uni-kl.de
9          H. Schoenemann  hannes@mathematik.uni-kl.de
10          A. Steenpass    steenpass@mathematik.uni-kl.de
11          S. Steidel      steidel@mathematik.uni-kl.de
12
13OVERVIEW:
14  A library for computing Groebner bases of ideals in the polynomial ring over
15  the rational numbers using modular methods.
16
17REFERENCES:
18  E. A. Arnold: Modular algorithms for computing Groebner bases.
19  J. Symb. Comp. 35, 403-419 (2003).
20
21  N. Idrees, G. Pfister, S. Steidel: Parallelization of Modular Algorithms.
22  J. Symb. Comp. 46, 672-684 (2011).
23
24PROCEDURES:
25  modStd(I);    standard basis of I using modular methods
26";
27
28LIB "poly.lib";
29LIB "modular.lib";
30
31proc modStd(ideal I, list #)
32"USAGE:   modStd(I[, exactness]); I ideal, exactness int
33RETURN:   a standard basis of I
34NOTE:     The procedure computes a standard basis of I (over the rational
35          numbers) by using modular methods.
36       @* An optional parameter 'exactness' can be provided.
37          If exactness = 1, the procedure computes a standard basis of I for
38          sure; if exactness = 0, it computes a standard basis of I
39          with high probability.
40SEE ALSO: modular
41EXAMPLE:  example modStd; shows an example"
42{
43    /* read optional parameter */
44    int exactness = 1;
45    if (size(#) > 0) {
46        /* For compatibility, we only test size(#) > 4. This can be changed to
47         * size(#) > 1 in the future. */
48        if (size(#) > 4 || typeof(#[1]) != "int") {
49            ERROR("wrong optional parameter");
50        }
51        exactness = #[1];
52    }
53
54    /* save options */
55    intvec opt = option(get);
56    option(redSB);
57
58    /* choose the right command */
59    string command = "groebner";
60    if (npars(basering) > 0) {
61        command = "Modstd::groebner_norm";
62    }
63
64    /* call modular() */
65    if (exactness) {
66        I = modular(command, list(I), primeTest_std,
67            deleteUnluckyPrimes_std, pTest_std, finalTest_std);
68    }
69    else {
70        I = modular(command, list(I), primeTest_std,
71            deleteUnluckyPrimes_std, pTest_std);
72    }
73
74    /* return the result */
75    attrib(I, "isSB", 1);
76    option(set, opt);
77    return(I);
78}
79example
80{
81    "EXAMPLE:";
82    echo = 2;
83    ring R1 = 0, (x,y,z,t), dp;
84    ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4;
85    ideal J = modStd(I);
86    J;
87    I = homog(I, t);
88    J = modStd(I);
89    J;
90
91    ring R2 = 0, (x,y,z), ds;
92    ideal I = jacob(x5+y6+z7+xyz);
93    ideal J = modStd(I, 0);
94    J;
95
96    ring R3 = 0, x(1..4), lp;
97    ideal I = cyclic(4);
98    ideal J1 = modStd(I, 1);   // default
99    ideal J2 = modStd(I, 0);
100    size(reduce(J1, J2));
101    size(reduce(J2, J1));
102}
103
104/* compute a normalized GB via groebner() */
105static proc groebner_norm(ideal I)
106{
107    I = simplify(groebner(I), 1);
108    attrib(I, "isSB", 1);
109    return(I);
110}
111
112/* test if the prime p is suitable for the input, i.e. it does not divide
113 * the numerator or denominator of any of the coefficients */
114static proc primeTest_std(int p, alias list args)
115{
116    /* erase zero generators */
117    ideal I = simplify(args[1], 2);
118
119    /* clear denominators and count the terms */
120    ideal J;
121    ideal K;
122    int n = ncols(I);
123    intvec sizes;
124    number cnt;
125    int i;
126    for(i = n; i > 0; i--) {
127        J[i] = cleardenom(I[i]);
128        cnt = leadcoef(J[i])/leadcoef(I[i]);
129        K[i] = numerator(cnt)*var(1)+denominator(cnt);
130    }
131    sizes = size(J[1..n]);
132
133    /* change to characteristic p */
134    def br = basering;
135    list lbr = ringlist(br);
136    if (typeof(lbr[1]) == "int") {
137        lbr[1] = p;
138    }
139    else {
140        lbr[1][1] = p;
141    }
142    def rp = ring(lbr);
143    setring(rp);
144    ideal Jp = fetch(br, J);
145    ideal Kp = fetch(br, K);
146
147    /* test if any coefficient is missing */
148    if (intvec(size(Kp[1..n])) != 2:n) {
149        setring(br);
150        return(0);
151    }
152    if (intvec(size(Jp[1..n])) != sizes) {
153        setring(br);
154        return(0);
155    }
156    setring(br);
157    return(1);
158}
159
160/* find entries in modresults which come from unlucky primes.
161 * For this, sort the entries into categories depending on their leading
162 * ideal and return the indices in all but the biggest category. */
163static proc deleteUnluckyPrimes_std(alias list modresults)
164{
165    int size_modresults = size(modresults);
166
167    /* sort results into categories.
168     * each category is represented by three entries:
169     * - the corresponding leading ideal
170     * - the number of elements
171     * - the indices of the elements
172     */
173    list cat;
174    int size_cat;
175    ideal L;
176    int i;
177    int j;
178    for (i = 1; i <= size_modresults; i++) {
179        L = lead(modresults[i]);
180        attrib(L, "isSB", 1);
181        for (j = 1; j <= size_cat; j++) {
182            if (size(L) == size(cat[j][1])
183                && size(reduce(L, cat[j][1])) == 0
184                && size(reduce(cat[j][1], L)) == 0) {
185                cat[j][2] = cat[j][2]+1;
186                cat[j][3][cat[j][2]] = i;
187                break;
188            }
189        }
190        if (j > size_cat) {
191            size_cat++;
192            cat[size_cat] = list();
193            cat[size_cat][1] = L;
194            cat[size_cat][2] = 1;
195            cat[size_cat][3] = list(i);
196        }
197    }
198
199    /* find the biggest categories */
200    int cat_max = 1;
201    int max = cat[1][2];
202    for (i = 2; i <= size_cat; i++) {
203        if (cat[i][2] > max) {
204            cat_max = i;
205            max = cat[i][2];
206        }
207    }
208
209    /* return all other indices */
210    list unluckyIndices;
211    for (i = 1; i <= size_cat; i++) {
212        if (i != cat_max) {
213            unluckyIndices = unluckyIndices + cat[i][3];
214        }
215    }
216    return(unluckyIndices);
217}
218
219/* test if 'command' applied to 'args' in characteristic p is the same as
220   'result' mapped to characteristic p */
221static proc pTest_std(string command, alias list args, alias ideal result,
222    int p)
223{
224    /* change to characteristic p */
225    def br = basering;
226    list lbr = ringlist(br);
227    if (typeof(lbr[1]) == "int") {
228        lbr[1] = p;
229    }
230    else {
231        lbr[1][1] = p;
232    }
233    def rp = ring(lbr);
234    setring(rp);
235    ideal Ip = fetch(br, args)[1];
236    ideal Gp = fetch(br, result);
237    attrib(Gp, "isSB", 1);
238
239    /* test if Ip is in Gp */
240    int i;
241    for (i = ncols(Ip); i > 0; i--) {
242        if (reduce(Ip[i], Gp, 1) != 0) {
243            setring(br);
244            return(0);
245        }
246    }
247
248    /* compute command(args) */
249    execute("Ip = "+command+"(Ip);");
250
251    /* test if Gp is in Ip */
252    for (i = ncols(Gp); i > 0; i--) {
253        if (reduce(Gp[i], Ip, 1) != 0) {
254            setring(br);
255            return(0);
256        }
257    }
258    setring(br);
259    return(1);
260}
261
262/* test if 'result' is a GB of the input ideal */
263static proc finalTest_std(string command, alias list args, ideal result)
264{
265    /* test if args[1] is in result */
266    attrib(result, "isSB", 1);
267    int i;
268    for (i = ncols(args[1]); i > 0; i--) {
269        if (reduce(args[1][i], result, 1) != 0) {
270            return(0);
271        }
272    }
273
274    /* test if result is a GB */
275    ideal G = std(result);
276    if (reduce_parallel(G, result)) {
277        return(0);
278    }
279    return(1);
280}
281
282/* return 1, if I_reduce is _not_ in G_reduce,
283 *        0, otherwise
284 * (same as size(reduce(I_reduce, G_reduce))).
285 * Uses parallelization. */
286static proc reduce_parallel(def I_reduce, def G_reduce)
287{
288    exportto(Modstd, I_reduce);
289    exportto(Modstd, G_reduce);
290    int size_I = ncols(I_reduce);
291    int chunks = Modular::par_range(size_I);
292    intvec range;
293    int i;
294    for (i = chunks; i > 0; i--) {
295        range = Modular::par_range(size_I, i);
296        task t(i) = "Modstd::reduce_task", list(range);
297    }
298    startTasks(t(1..chunks));
299    waitAllTasks(t(1..chunks));
300    int result = 0;
301    for (i = chunks; i > 0; i--) {
302        if (getResult(t(i))) {
303            result = 1;
304            break;
305        }
306    }
307    kill I_reduce;
308    kill G_reduce;
309    return(result);
310}
311
312/* compute a chunk of reductions for reduce_parallel */
313static proc reduce_task(intvec range)
314{
315    int result = 0;
316    int i;
317    for (i = range[1]; i <= range[2]; i++) {
318        if (reduce(I_reduce[i], G_reduce, 1) != 0) {
319            result = 1;
320            break;
321        }
322    }
323    return(result);
324}
325
326////////////////////////////////////////////////////////////////////////////////
327/*
328 * The following procedures are kept for backward compatibility with the old
329 * version of modstd.lib. As of now (May 2014), they are still needed in
330 * modnormal.lib, modwalk.lib, and symodstd.lib. They can be removed here as
331 * soon as they are not longer needed in these libraries.
332 */
333
334LIB "parallel.lib";
335
336static proc mod_init()
337{
338   newstruct("idealPrimeTest", "ideal Ideal");
339}
340
341static proc redFork(ideal I, ideal J, int n)
342{
343   attrib(J,"isSB",1);
344   return(reduce(I,J,1));
345}
346
347proc isIncluded(ideal I, ideal J, list #)
348"USAGE:  isIncluded(I,J); I,J ideals
349RETURN:  1 if J includes I,
350@*       0 if there is an element f in I which does not reduce to 0 w.r.t. J.
351EXAMPLE: example isIncluded; shows an example
352"
353{
354   def R = basering;
355   setring R;
356
357   attrib(J,"isSB",1);
358   int i,j,k;
359
360   if(size(#) > 0)
361   {
362      int n = #[1];
363      if(n >= ncols(I)) { n = ncols(I); }
364      if(n > 1)
365      {
366         for(i = 1; i <= n - 1; i++)
367         {
368            //link l(i) = "MPtcp:fork";
369            link l(i) = "ssi:fork";
370            open(l(i));
371
372            write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1)));
373         }
374
375         int t = timer;
376         if(reduce(I[ncols(I)], J, 1) != 0)
377         {
378            for(i = 1; i <= n - 1; i++)
379            {
380               close(l(i));
381            }
382            return(0);
383         }
384         t = timer - t;
385         if(t > 60) { t = 60; }
386         int i_sleep = system("sh", "sleep "+string(t));
387
388         j = ncols(I) - n;
389
390         while(j >= 0)
391         {
392            for(i = 1; i <= n - 1; i++)
393            {
394               if(status(l(i), "read", "ready"))
395               {
396                  if(read(l(i)) != 0)
397                  {
398                     for(i = 1; i <= n - 1; i++)
399                     {
400                        close(l(i));
401                     }
402                     return(0);
403                  }
404                  else
405                  {
406                     if(j >= 1)
407                     {
408                        write(l(i), quote(redFork(eval(I[j]), eval(J), 1)));
409                        j--;
410                     }
411                     else
412                     {
413                        k++;
414                        close(l(i));
415                     }
416                  }
417               }
418            }
419            if(k == n - 1)
420            {
421               j--;
422            }
423            i_sleep = system("sh", "sleep "+string(t));
424         }
425         return(1);
426      }
427   }
428
429   for(i = ncols(I); i >= 1; i--)
430   {
431      if(reduce(I[i],J,1) != 0){ return(0); }
432   }
433   return(1);
434}
435example
436{ "EXAMPLE:"; echo = 2;
437   ring r=0,(x,y,z),dp;
438   ideal I = x+1,x+y+1;
439   ideal J = x+1,y;
440   isIncluded(I,J);
441   isIncluded(J,I);
442   isIncluded(I,J,4);
443
444   ring R = 0, x(1..5), dp;
445   ideal I1 = cyclic(4);
446   ideal I2 = I1,x(5)^2;
447   isIncluded(I1,I2,4);
448}
449
450proc deleteUnluckyPrimes(list T, list L, int ho, list #)
451"USAGE:  deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer
452RETURN:  lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #),
453         lT ideal
454NOTE:    - if ho = 1, the polynomials in T are homogeneous, else ho = 0,
455@*       - lT is prevalent, i.e. the most appearing leading ideal in T
456EXAMPLE: example deleteUnluckyPrimes; shows an example
457"
458{
459   ho = ((ho)||(ord_test(basering) == -1));
460   int j,k,c;
461   intvec hl,hc;
462   ideal cT,lT,cK;
463   lT = lead(T[size(T)]);
464   attrib(lT,"isSB",1);
465   if(!ho)
466   {
467      for(j = 1; j < size(T); j++)
468      {
469         cT = lead(T[j]);
470         attrib(cT,"isSB",1);
471         if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
472         {
473            cK = cT;
474            c++;
475         }
476      }
477      if(c > size(T) div 2){ lT = cK; }
478   }
479   else
480   {
481      hl = hilb(lT,1);
482      for(j = 1; j < size(T); j++)
483      {
484         cT = lead(T[j]);
485         attrib(cT,"isSB",1);
486         hc = hilb(cT,1);
487         if(hl == hc)
488         {
489            for(k = 1; k <= size(lT); k++)
490            {
491               if(lT[k] < cT[k]) { lT = cT; c++; break; }
492               if(lT[k] > cT[k]) { c++; break; }
493            }
494         }
495         else
496         {
497            if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; }
498         }
499      }
500   }
501
502   int addList;
503   if(size(#) > 0) { list M = #; addList = 1; }
504   j = 1;
505   attrib(lT,"isSB",1);
506   while((j <= size(T))&&(c > 0))
507   {
508      cT = lead(T[j]);
509      attrib(cT,"isSB",1);
510      if((size(reduce(cT,lT)) != 0)||(size(reduce(lT,cT)) != 0))
511      {
512         T = delete(T,j);
513         if(j == 1)
514         {
515            L = L[2..size(L)];
516            if(addList == 1) { M = M[2..size(M)]; }
517         }
518         else
519         {
520            if(j == size(L))
521            {
522               L = L[1..size(L)-1];
523               if(addList == 1) { M = M[1..size(M)-1]; }
524            }
525            else
526            {
527               L = L[1..j-1],L[j+1..size(L)];
528               if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; }
529            }
530         }
531         j--;
532      }
533      j++;
534   }
535
536   for(j = 1; j <= size(L); j++)
537   {
538      L[j] = bigint(L[j]);
539   }
540
541   if(addList == 0) { return(list(T,L,lT)); }
542   if(addList == 1) { return(list(T,L,M,lT)); }
543}
544example
545{ "EXAMPLE:"; echo = 2;
546   list L = 2,3,5,7,11;
547   ring r = 0,(y,x),Dp;
548   ideal I1 = 2y2x,y6;
549   ideal I2 = yx2,y3x,x5,y6;
550   ideal I3 = y2x,x3y,x5,y6;
551   ideal I4 = y2x,11x3y,x5;
552   ideal I5 = y2x,yx3,x5,7y6;
553   list T = I1,I2,I3,I4,I5;
554   deleteUnluckyPrimes(T,L,1);
555   list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5);
556   deleteUnluckyPrimes(T,L,1,P);
557}
558
559proc primeTest(def II, bigint p)
560{
561   if(typeof(II) == "string")
562   {
563      ideal I = `II`.Ideal;
564   }
565   else
566   {
567      ideal I = II;
568   }
569
570   I = simplify(I, 2);   // erase zero generators
571
572   int i,j;
573   poly f;
574   number cnt;
575   for(i = 1; i <= size(I); i++)
576   {
577      f = cleardenom(I[i]);
578      if(f == 0) { return(0); }
579      cnt = leadcoef(I[i])/leadcoef(f);
580      if((bigint(numerator(cnt)) mod p) == 0) { return(0); }
581      if((bigint(denominator(cnt)) mod p) == 0) { return(0); }
582      for(j = size(f); j > 0; j--)
583      {
584         if((bigint(leadcoef(f[j])) mod p) == 0) { return(0); }
585      }
586   }
587   return(1);
588}
589
590proc primeList(ideal I, int n, list #)
591"USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
592         n integer
593RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
594         < L[size(L)] union with L) such that none of these primes divides any
595         coefficient occuring in I
596NOTE:    The number of cores to use can be defined by ncores, default is 1.
597EXAMPLE: example primeList; shows an example
598"
599{
600   intvec L;
601   int i,p;
602   int ncores = 1;
603
604//-----------------  Initialize optional parameter ncores  ---------------------
605   if(size(#) > 0)
606   {
607      if(size(#) == 1)
608      {
609         if(typeof(#[1]) == "int")
610         {
611            ncores = #[1];
612            # = list();
613         }
614      }
615      else
616      {
617         ncores = #[2];
618      }
619   }
620
621   if(size(#) == 0)
622   {
623      p = 2147483647;
624      while(!primeTest(I,p))
625      {
626         p = prime(p-1);
627         if(p == 2) { ERROR("no more primes"); }
628      }
629      L[1] = p;
630   }
631   else
632   {
633      L = #[1];
634      p = prime(L[size(L)]-1);
635      while(!primeTest(I,p))
636      {
637         p = prime(p-1);
638         if(p == 2) { ERROR("no more primes"); }
639      }
640      L[size(L)+1] = p;
641   }
642   if(p == 2) { ERROR("no more primes"); }
643   if(ncores == 1)
644   {
645      for(i = 2; i <= n; i++)
646      {
647         p = prime(p-1);
648         while(!primeTest(I,p))
649         {
650            p = prime(p-1);
651            if(p == 2) { ERROR("no more primes"); }
652         }
653         L[size(L)+1] = p;
654      }
655   }
656   else
657   {
658      int neededSize = size(L)+n-1;;
659      list parallelResults;
660      list arguments;
661      int neededPrimes = neededSize-size(L);
662      idealPrimeTest Id;
663      Id.Ideal = I;
664      export(Id);
665      while(neededPrimes > 0)
666      {
667         arguments = list();
668         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
669            *ncores; i > 0; i--)
670         {
671            p = prime(p-1);
672            if(p == 2) { ERROR("no more primes"); }
673            arguments[i] = list("Id", p);
674         }
675         parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores);
676         for(i = size(arguments); i > 0; i--)
677         {
678            if(parallelResults[i])
679            {
680               L[size(L)+1] = arguments[i][2];
681            }
682         }
683         neededPrimes = neededSize-size(L);
684      }
685      kill Id;
686      if(size(L) > neededSize)
687      {
688         L = L[1..neededSize];
689      }
690   }
691   return(L);
692}
693example
694{ "EXAMPLE:"; echo = 2;
695   ring r = 0,(x,y,z),dp;
696   ideal I = 2147483647x+y, z-181;
697   intvec L = primeList(I,10);
698   size(L);
699   L[1];
700   L[size(L)];
701   L = primeList(I,5,L);
702   size(L);
703   L[size(L)];
704}
705
706////////////////////////////// further examples ////////////////////////////////
707
708/*
709ring r = 0, (x,y,z), lp;
710poly s1 = 5x3y2z+3y3x2z+7xy2z2;
711poly s2 = 3xy2z2+x5+11y2z2;
712poly s3 = 4xyz+7x3+12y3+1;
713poly s4 = 3x3-4y3+yz2;
714ideal i =  s1, s2, s3, s4;
715
716ring r = 0, (x,y,z), lp;
717poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
718poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
719poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
720poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
721ideal i =  s1, s2, s3, s4;
722
723ring r = 0, (x,y,z), lp;
724poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
725poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
726poly s3 = 8x3 + 12y3 + xz2 + 3;
727poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
728ideal i = s1, s2, s3, s4;
729
730int n = 6;
731ring r = 0,(x(1..n)),lp;
732ideal i = cyclic(n);
733ring s = 0, (x(1..n),t), lp;
734ideal i = imap(r,i);
735i = homog(i,t);
736
737ring r = 0, (x(1..4),s), (dp(4),dp);
738poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
739poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
740poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
741poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4);
742poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
743ideal i =  s1, s2, s3, s4, s5;
744
745ring r = 0, (x,y,z), ds;
746int a = 16;
747int b = 15;
748int c = 4;
749int t = 1;
750poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
751         +x^(c-2)*y^c*(y2+t*x)^2;
752ideal i = jacob(f);
753
754ring r = 0, (x,y,z), ds;
755int a = 25;
756int b = 25;
757int c = 5;
758int t = 1;
759poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3
760         +x^(c-2)*y^c*(y2+t*x)^2;
761ideal i = jacob(f),f;
762
763ring r = 0, (x,y,z), ds;
764int a = 10;
765poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
766ideal i = jacob(f);
767
768ring r = 0, (x,y,z), ds;
769int a = 6;
770int b = 8;
771int c = 10;
772int alpha = 5;
773int beta = 5;
774int t = 1;
775poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
776         +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
777ideal i = jacob(f);
778*/
779
Note: See TracBrowser for help on using the repository browser.