source: git/Singular/LIB/nfmodsyz.lib @ a3c7bd

spielwiese
Last change on this file since a3c7bd was a3c7bd, checked in by dere <dkifle16@…>, 8 years ago
fix: typo and delete commented procedures
  • Property mode set to 100644
File size: 18.0 KB
Line 
1////////////////////////////////////////////////////////////////////////////////
2version="Id";  // $Id$
3category="Commutative Algebra";
4info="
5
6LIBRARY:   nfmodsyz.lib    Syzygy modules of submodules of free modules over
7                           algebraic number fields
8AUTHORS:   D.K. Boku       boku@mathematik.uni-kl.de
9@*         W. Decker       decker@mathematik.uni-kl.de
10@*         C. Fieker       fieker@mathematik.uni-kl.de
11
12OVERVIEW:
13  A library for computing the syzygy module of a given submodule I in a polynomial ring
14  over an algebraic number field Q(t), where t is an algebraic number, using modular methods.
15  For the case Q(t)=Q, that is, where t is an element of Q, we compute, following
16  [1], the syzygy module of I as follows: For a submodule I of A^m with A = Q[X], we first
17  choose a sufficiently large set of primes P and compute the reduced Groebner basis of the
18  syzygy module of I_p, for each p in P, in parallel. We then use the Chinese remainder
19  algorithm and rational reconstruction to obtain the syzygy module of I over Q.
20  For the case where t is not in Q, we compute, following [2], the syzygy module of I as
21  follows:
22  Let f be the minimal polynomial of t. For a submodule I in A^m with A = Q(t)[X], we map I
23  to a submodule I' in A^m with A = (Q[t]/<f>)[X] via the map sending t to t + <f>. We first
24  choose a prime p such that f has at least two factors in characteristic p. For each
25  factor f_{i,p} of f_p:= (f mod p), we set I'_{i,p} := (I'_p mod f_{i,p}). We then
26  compute the reduced Groebner bases G'_i of the syzygy modules of I'_{i,p} over
27  F_p[t]/<f_{i,p}> and combine the G'_i to G_p (the syzygy module of I'_p) using chinese
28  remaindering for polynomials. As described in [2], the procedure is repeated for many primes
29  p, where we compute the G_p in parallel until the number of primes is sufficiently large to
30  recover the correct generating set for the syzygy module G' of I' which is, considered over
31  Q(t), also a generating set for the syzygy module of I.
32
33REFERENCES:
34  [1] E. A. Arnold: Modular algorithms for computing Groebner bases.
35      J. Symb. Comp. 35, 403-419 (2003).
36  [2] D. Boku, W. Decker, C. Fieker, and A. Steenpass. Groebner bases over algebraic
37      number fields. In: Proceedings of the 2015 International Workshop on Parallel
38      Symb. Comp. PASCO'15, pages 16-24 (2015).
39
40PROCEDURES:
41  nfmodSyz(I);          syzygy module of I over algebraic number field using modular
42                        methods
43";
44
45LIB "nfmodstd.lib";
46
47////////////////////////////////////////////////////////////////////////////////
48
49static proc testPrime(int p, list args)
50{
51    /*
52     * test whether a prime p divides the denominator(s)
53     * and leading coefficients of generating set of ideal
54     */
55    int i,j,k;
56    vector f;
57    number num;
58    module I = args[1];
59    bigint d1,d2,d3;
60    for(i = 1; i <= ncols(I); i++)
61    {
62        f = cleardenom(I[i]);
63        if(f == 0)
64        {
65            return(0);
66        }
67        num = leadcoef(I[i])/leadcoef(f);
68        d1 = bigint(numerator(num));
69        d2 = bigint(denominator(num));
70        if( (d1 mod p) == 0)
71        {
72            return(0);
73        }
74        if((d2 mod p) == 0)
75        {
76            return(0);
77        }
78        for(j = nrows(f); j > 0; j--)
79        {
80            for(k=1;k<=size(f[j]);k++)
81            {
82                d3 = bigint(leadcoef(f[j][k]));
83                if((d3 mod p) == 0)
84                {
85                    return(0);
86                }
87            }
88        }
89    }
90    return(1);
91}
92
93////////////////////////////////////////////////////////////////////////////////
94/* return 1 if the number of factors are in the required bound , 0 else */
95
96static proc minpolyTask(poly f,int p)
97{
98    /*
99     * bound for irreducible factor(s) of (f mod p)
100     * see testfact()
101     */
102    int nr,k,ur;
103    ur=deg(f);
104    list L=factmodp(f,p);
105    if(degtest(L[2])==1)
106    {
107        // now each factor is squarefree
108        if(ur<=3)
109        {
110            return(1);
111        }
112        else
113        {
114            nr = testfact(ur);
115            k=ncols(L[1]);
116            if(nr < k && k < (ur-nr)) // set a bound for k
117            {
118                return(1);
119            }
120        }
121    }
122    return(0);
123}
124
125////////////////////////////////////////////////////////////////////////////////
126/* return 1 if both testPrime(p,J) and minpolyTask(f,p) is true, 0 else */
127
128static proc PrimeTestTask_syz(int p, list L)
129{
130    /* L=list(I), I=J,f; J ideal , f minpoly */
131    int sz,nr;
132    module J = L[1];
133    sz=ncols(J);
134    def f=J[sz];
135    poly g = f[1];
136    if(!testPrime(p,list(J)) or !minpolyTask(g,p))
137    {
138        return(0);
139    }
140    return(1);
141}
142
143////////////////////////////////////////////////////////////////////////////////
144/* compute factors of f mod p with multiplicity */
145
146static proc factmodp(poly f, int p)
147{
148    def R=basering;
149    list l=ringlist(R);
150    l[1]=p;
151    def S=ring(l);
152    setring S;
153    list L=factorize(imap(R,f),2);
154    ideal J=L[1];
155    intvec v=L[2];
156    list scx=J,v;
157    setring R;
158    return(imap(S,scx));
159    kill S;
160}
161
162////////////////////////////////////////////////////////////////////////////////
163/* set a bound for number of factors w.r.t degree nr*/
164
165static proc testfact(int nr)
166{
167    // nr must be greater than 3
168    int i;
169    if(nr>3 and nr<=5)
170    {
171        i=1;
172    }
173    if(nr>5 and nr<=10)
174    {
175        i=2;
176    }
177    if(nr>10 and nr<=15)
178    {
179        i=3;
180    }
181    if(nr>15 and nr<=20)
182    {
183        i=4;
184    }
185    if(nr>20 and nr<=25)
186    {
187        i=5;
188    }
189    if(nr>25 and nr<=30)
190    {
191        i=6;
192    }
193    if(nr>30)
194    {
195        i=10;
196    }
197    return(i);
198}
199
200///////////////////////////////////////////////////////////////////////////////
201// return 1 if v[i]>1 , 0 else
202
203static proc degtest(intvec v)
204{
205    for(int j=1;j<=nrows(v);j++)
206    {
207        if(v[j]>1)
208        {
209            return(0);
210        }
211    }
212    return(1);
213}
214
215////////////////////////////////////////////////////////////////////////////////
216
217static proc check_leadmonom_and_size(list L)
218{
219    /*
220     * compare the size of ideals in the list and
221     * check the corresponding leading monomials
222     * size(L)>=2
223     */
224    def J=L[1];
225    int i=size(L);
226    int sc=ncols(J);
227    int j,k;
228    def g=leadmonom(J[1]);
229    for(j=1;j<=i;j++)
230    {
231        if(ncols(L[j])!=sc)
232        {
233            return(0);
234        }
235    }
236    for(k=2;k<=i;k++)
237    {
238        for(j=1;j<=sc;j++)
239        {
240            if(leadmonom(J[j])!=leadmonom(L[k][j]))
241            {
242                return(0);
243            }
244        }
245    }
246    return(1);
247}
248
249////////////////////////////////////////////////////////////////////////////////
250
251static proc LiftPolyCRT_syz(def I)
252{
253    /*
254     * compute syz for each factor and combine this result
255     * to modulo minpoly via CRT for poly over char p>0
256     */
257    def sl;
258    int u,in,j;
259    list LL,Lk,T2;
260    module J,II;
261    vector f;
262    u=ncols(I);
263    J=I[1..u-1];
264    f=I[u];
265    poly ff = f[1];
266    ideal K=factorize(ff,1);
267    in=ncols(K);
268    def Ls = basering;
269    list l = ringlist(Ls);
270    if(l[3][1][1]=="c")
271    {
272        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
273        list(list(l[3][size(l[3])]))+list(ideal(0));
274        l[2] = delete(l[2],size(l[2]));
275        l[3] = delete(l[3],size(l[3]));
276    }
277    else
278    {
279        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
280        list(list(l[3][size(l[3])-1]))+list(ideal(0));
281        l[2] = delete(l[2],size(l[2]));
282        l[3] = delete(l[3],size(l[3])-1);
283    }
284
285    def S1 = ring(l);
286    setring S1;
287    number Num= number(imap(Ls,ff));
288    list l = ringlist(S1);
289    l[1][4][1] = Num;
290    S1 = ring(l);
291    setring S1;
292    ideal K = imap(Ls,K);
293    def S2;
294    module II;
295    number Num;
296    /* ++++++ if minpoly is irreducible then K will be the zero ideal +++ */
297    if(size(K)==0)
298    {
299        module M = syz(imap(Ls,J));
300        if(size(M)==0)
301        {
302            setring Ls;
303            return(module([0]));
304        }
305        II = normalize(M);
306    }
307    else
308    {
309        for(j=1;j<=in;j++)
310        {
311            LL[j]=K[j];
312            Num = number(K[j]);
313            T2 = ringlist(S1);
314            T2[1][4][1] = Num;
315            S2 = ring(T2);
316            setring S2;
317            module M = syz(imap(Ls,J));
318            if(size(M)==0)
319            {
320                setring Ls;
321                return(module([0]));
322                break;
323            }
324            setring S1;
325            Lk[j] = imap(S2,M);
326        }
327
328        if(check_leadmonom_and_size(Lk))
329        {
330            // apply CRT for polynomials
331            setring Ls;
332            II =chinrempoly(imap(S1,Lk),imap(S1,LL));
333            setring S1;
334            II = normalize(imap(Ls,II));
335        }
336        else
337        {
338            setring S1;
339            II=[0];
340        }
341     }
342     setring Ls;
343     return(imap(S1,II));
344}
345
346////////////////////////////////////////////////////////////////////////////////
347
348static proc final_Test_syz(string command, alias list args, def result)
349{
350    /*
351     * test if the set generating 'result' also generates the syzygy module
352     * of args[1] in characteristic zero
353     */
354    def Ls = basering;
355    def Ip = args[1];
356    vector f;
357    int u=ncols(Ip);
358    module J=Ip[1..u-1];
359    f=Ip[u];
360    poly ff = f[1];
361    list l = ringlist(Ls);
362
363    if(l[3][1][1]=="c")
364    {
365        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
366        list(list(l[3][size(l[3])]))+list(ideal(0));
367        l[2] = delete(l[2],size(l[2]));
368        l[3] = delete(l[3],size(l[3]));
369    }
370    else
371    {
372        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
373        list(list(l[3][size(l[3])-1]))+list(ideal(0));
374        l[2] = delete(l[2],size(l[2]));
375        l[3] = delete(l[3],size(l[3])-1);
376    }
377
378    def S1 = ring(l);
379    setring S1;
380    number Num= number(imap(Ls,ff));
381    list l = ringlist(S1);
382    l[1][4][1] = Num;
383    S1 = ring(l);
384    setring S1;
385    def result2 = imap(Ls,result);
386    def M = imap(Ls,J);
387    if(size(result2)==0)
388    {
389        return(1);
390    }
391    else
392    {
393        if(size(module(matrix(M)*matrix(result2)))!=0)
394        {
395            return(0);
396        }
397        return(1);
398    }
399}
400
401////////////////////////////////////////////////////////////////////////////////
402
403static proc final_test(string command, alias list args, def result)
404{
405    /*
406     * test if the set  generating 'result' also generates the syzygy module
407     * of args[1] in characteristic zero
408     */
409    module M=args[1];
410    if(size(result)==0)
411    {
412        return(1);
413    }
414    else
415    {
416        if(size(module(matrix(M)*matrix(result)))!=0)
417        {
418            return(0);
419        }
420        return(1);
421    }
422}
423
424////////////////////////////////////////////////////////////////////////////////
425// ------------------------ test in characteristic p ------------
426static proc pTest_syzmod(string command, list args, def result, int p)
427{
428     /*
429      * This procedure performs the first test in positive characteristic to
430      * verify whether the set generating 'result' also generates the syzygy
431      * module of the submodule args[1]. Note that this test works only
432      * over Z_p
433      */
434     def br = basering;
435     if(size(result)==0)
436     {
437         return(1);
438     }
439     list lbr = ringlist(br);
440     if (typeof(lbr[1]) == "int")
441     {
442         lbr[1] = p;
443     }
444     else
445     {
446         lbr[1][1] = p;
447     }
448     def rp = ring(lbr);
449     setring(rp);
450     module Jp = imap(br, args)[1];
451     module Gp = imap(br, result);
452     module Ip = syz(Jp);
453     // test if Ip is contained in Gp
454     attrib(Gp, "isSB", 1);
455     for (int i = ncols(Ip); i > 0; i--)
456     {
457          if (reduce(Ip[i], Gp, 1) != 0)
458          {
459              setring(br);
460              return(0);
461          }
462     }
463     // test if Gp is contained in syz(Jp)
464     if(size(module(matrix(Jp)*matrix(Gp)))!=0)
465     {
466         setring br;
467         return(0);
468     }
469     setring br;
470     return(1);
471}
472
473////////////////////////////////////////////////////////////////////////////////
474 // ------------------------ test in characteristic p ------------
475static proc pTest_syz(string command, list args, def result, int p)
476{
477     /*
478      * This procedure performs the first test in positive characteristic to
479      * verify whether the set generating 'result' also generates the syzygy
480      * module of args[1]. Note that this test works only over Z_p(t) where
481      * t is an algebraic number which is not in Z_p.
482     */
483
484     def br = basering;
485     if(size(result)==0)
486     {
487         return(1);
488     }
489     list lbr = ringlist(br);
490     if (typeof(lbr[1]) == "int")
491     {
492         lbr[1] = p;
493     }
494     else
495     {
496         lbr[1][1] = p;
497     }
498     def rp = ring(lbr);
499     setring(rp);
500     def Ip = imap(br, args)[1];
501
502     int u,in,j,i;
503     list LL,Lk,T2;
504     module J,II;
505     vector f;
506     u=ncols(Ip);
507     J=Ip[1..u-1];
508     f=Ip[u];
509     poly ff = f[1];
510     ideal K=factorize(ff,1);
511     in=ncols(K);
512     def Ls = basering;
513     list l = ringlist(Ls);
514     if(l[3][1][1]=="c")
515     {
516        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
517        list(list(l[3][size(l[3])]))+list(ideal(0));
518        l[2] = delete(l[2],size(l[2]));
519        l[3] = delete(l[3],size(l[3]));
520     }
521     else
522     {
523        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
524        list(list(l[3][size(l[3])-1]))+list(ideal(0));
525        l[2] = delete(l[2],size(l[2]));
526        l[3] = delete(l[3],size(l[3])-1);
527     }
528
529     def S1 = ring(l);
530     setring S1;
531     number Num= number(imap(Ls,ff));
532     list l = ringlist(S1);
533     l[1][4][1] = Num;
534     S1 = ring(l);
535     setring S1;
536     ideal K = imap(Ls,K);
537     module Jp = imap(Ls,J);
538     def S2;
539     module Ip;
540     number Num;
541     /* ++++++ if the minpoly is irreducible then K = ideal(0) +++ */
542     if(size(K)==0)
543     {
544          module M = syz(Jp);
545          Ip = normalize(M);
546     }
547     else
548     {
549         for(j=1;j<=ncols(K);j++)
550         {
551              LL[j]=K[j];
552              Num = number(K[j]);
553              T2 = ringlist(S1);
554              T2[1][4][1] = Num;
555              S2 = ring(T2);
556              setring S2;
557              module M = syz(imap(Ls,J));
558              setring S1;
559              Lk[j]= imap(S2,M);
560         }
561         if(check_leadmonom_and_size(Lk))
562         {
563              // apply CRT for polynomials
564              setring Ls;
565              II =chinrempoly(imap(S1,Lk),imap(S1,LL));
566              setring S1;
567              Ip = normalize(imap(Ls,II));
568         }
569         else
570         {
571              setring S1;
572              Ip=[0];
573         }
574     }
575     setring S1;
576     module Gp = imap(br, result);
577     // test if Ip is contained in Gp
578     attrib(Gp, "isSB", 1);
579     for (i = ncols(Ip); i > 0; i--)
580     {
581          if (reduce(Ip[i], Gp, 1) != 0)
582          {
583              setring(br);
584              return(0);
585          }
586     }
587     // test if Gp is contained in syz(Jp)
588     if(size(module(matrix(Jp)*matrix(Gp)))!=0)
589     {
590         setring br;
591         return(0);
592     }
593     setring br;
594     return(1);
595}
596
597 ////////////////////////////////////////////////////////////////////////////////
598
599 static proc cleardenomIdeal(def I)
600 {
601     int t=ncols(I);
602     if(size(I)==0)
603     {
604         return(I);
605     }
606     else
607     {
608         for(int i=1;i<=t;i++)
609         {
610             I[i]=cleardenom(I[i]);
611         }
612     }
613     return(I);
614 }
615
616////////////////////////////////////////////////////////////////////////////////
617
618 static proc modStdparallelized_syzSB(module I, list #)
619 {
620     /* save options */
621     intvec opt = option(get);
622     option(redSB);
623     option(returnSB);
624     /*------ if these options are set, the Singular command syz returns the
625      reduced Groebner basis of I ---------------------------------------*/
626
627     // apply modular command from modular.lib
628     if(size(#)>0)
629     {
630         I = modular("syz", list(I), testPrime, Modstd::deleteUnluckyPrimes_std,
631             pTest_syzmod, final_test, 536870909);
632     }
633     else
634     {
635         I = modular("Nfmodsyz::LiftPolyCRT_syz", list(I), PrimeTestTask_syz,
636             Modstd::deleteUnluckyPrimes_std,pTest_syz, final_Test_syz,536870909);
637     }
638     attrib(I, "isSB", 1);
639     option(set,opt);
640     return(I);
641 }
642
643////////////////////////////////////////////////////////////////////////////////
644/* main procedure */
645proc nfmodSyz(def I)
646"USAGE:  nfmodSyz(I); I ideal or module
647RETURN:  syzygy module of I over an algebraic number field
648SEE ALSO: syz
649EXAMPLE: example nfmodSyz; shows an example
650"
651{
652     if(typeof(I)!="ideal" and typeof(I)!="module")
653     {
654        ERROR("type of input must be either ideal or module");
655     }
656     else
657     {
658        module F = I;
659        kill I;
660        module I = F;
661     }
662     def Rbs=basering;
663     poly f;
664     int n=nvars(Rbs);
665     if(size(I)==0)
666     {
667         return(module([0]));
668     }
669     if(npars(Rbs)==0)
670     {
671        module M = modStdparallelized_syzSB(I,1); //if algebraic number is in Q
672        return(M);
673     }
674
675     def S;
676     list rl=ringlist(Rbs);
677     f=rl[1][4][1];
678
679     if(rl[3][1][1]!="c")
680     {
681        rl[2] = rl[2] + rl[1][2];
682        rl[3] = insert(rl[3], rl[1][3][1],1);
683        rl[1] = rl[1][1];
684     }
685     else
686     {
687        rl[2] = rl[2] + rl[1][2];
688        rl[3][size(rl[3])+1] = rl[1][3][1];
689        rl[1] = rl[1][1];
690     }
691
692     S = ring(rl);
693     setring S;
694     poly f=imap(Rbs,f);
695     def I=imap(Rbs,I);
696     I = simplify(I,2); // eraze the zero generatos
697     if(f==0)
698     {
699         ERROR("minpoly must be non-zero");
700     }
701     I=I,f;
702     def J_I = modStdparallelized_syzSB(I);
703     setring Rbs;
704     def J=imap(S,J_I);
705     J=simplify(J,2);
706     return(J);
707}
708example
709{ "EXAMPLE:"; echo = 2;
710    ring r1 =(0,a),(x,y),(c,dp);
711    minpoly = (a^3+2a+7);
712    module M1 = [(a/2+1)*y, 3*x-a*y],
713                [y-x,y2],
714                [x2-xy, ax-y];
715    nfmodSyz(M1);
716    ring r2 = (0,a),(x,y,z),(dp,c);
717    minpoly = (a3+a+1);
718    module M2 = [x2z+x+(-a)*y,z2+(a+2)*x],
719                [y2+(a)*z+(a),(a+3)*z3+(-a)*x2],
720                [-xz+(a2+3)*yz,xy+(a2)*z];
721    nfmodSyz(M2);
722    ring r3=0,(x,y),dp; // ring without parameter
723    module M3 = [x2 + y, xy], [-7y, 2x], [x2-y, 0];
724    nfmodSyz(M3);
725    ring r4=0,(x,y),(c,dp); // ring without parameter
726    module M4 = [xy, x-y],
727                [x2 + y, 5y],
728                [- 7y, 2x],
729                [x2-y, 0];
730    nfmodSyz(M4);
731}
732
Note: See TracBrowser for help on using the repository browser.