source: git/Singular/LIB/nfmodsyz.lib @ 11fd9a

spielwiese
Last change on this file since 11fd9a was 11fd9a, checked in by dere <dkifle16@…>, 8 years ago
Add:new library nfmodsyz.lib
  • Property mode set to 100644
File size: 20.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 over an
14  algebraic number field Q(t), where t is an algebraic number, using modular methods.
15  For the case Q(t)=Q, that is, t is an element of Q, we compute, following
16  [1], the syzygy module of I as follows: For I submodule of A^m with A = Q[X],
17  as described in [1], we first choose sufficiently large set of primes P and compute
18  the reduced Groebner basis of the syzygy module of I_p, for each p in P, in parallel.
19  We then use the Chinese remainder algorithm and rational reconstruction to obtain the syzygy
20  module of I over Q. For the case t is not in Q, we compute, following [2], the syzygy module
21  of I as follows:
22  Let f be the minimal polynomial of t. For I, I' submodules in A^m with A = Q(t)[X],
23  (Q[t]/<f>)[X] respectively, we map I to I' 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 also a generating
31  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 'result' generates the syzygy module of args[1]
352     * 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 'result' generates the syzygy module of args[1]
407     * 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 makes the first test in positive characteristic to verify whether the
430      * set 'result' generates the submodule args[1]. Note that this test works only over Z_p
431      */
432     def br = basering;
433     if(size(result)==0)
434     {
435         return(1);
436     }
437     list lbr = ringlist(br);
438     if (typeof(lbr[1]) == "int")
439     {
440         lbr[1] = p;
441     }
442     else
443     {
444         lbr[1][1] = p;
445     }
446     def rp = ring(lbr);
447     setring(rp);
448     module Jp = imap(br, args)[1];
449     module Gp = imap(br, result);
450     module Ip = syz(Jp);
451     // test if Ip is contained in Gp
452     attrib(Gp, "isSB", 1);
453     for (int i = ncols(Ip); i > 0; i--)
454     {
455          if (reduce(Ip[i], Gp, 1) != 0)
456          {
457              setring(br);
458              return(0);
459          }
460     }
461     // test if Gp is contained in syz(Jp)
462     if(size(module(matrix(Jp)*matrix(Gp)))!=0)
463     {
464         setring br;
465         return(0);
466     }
467     setring br;
468     return(1);
469}
470
471////////////////////////////////////////////////////////////////////////////////
472 // ------------------------ test in characteristic p ------------
473static proc pTest_syz(string command, list args, def result, int p)
474{
475     /*
476      * This procedure makes the first test in positive characteristic to verify whether the
477      * set 'result' generates the submodule args[1]. Note that this test works only over
478      * Z_p(t) where t is an algebraic number which is not Z_p.
479     */
480
481     def br = basering;
482     if(size(result)==0)
483     {
484         return(1);
485     }
486     list lbr = ringlist(br);
487     if (typeof(lbr[1]) == "int")
488     {
489         lbr[1] = p;
490     }
491     else
492     {
493         lbr[1][1] = p;
494     }
495     def rp = ring(lbr);
496     setring(rp);
497     def Ip = imap(br, args)[1];
498
499     int u,in,j,i;
500     list LL,Lk,T2;
501     module J,II;
502     vector f;
503     u=ncols(Ip);
504     J=Ip[1..u-1];
505     f=Ip[u];
506     poly ff = f[1];
507     ideal K=factorize(ff,1);
508     in=ncols(K);
509     def Ls = basering;
510     list l = ringlist(Ls);
511     if(l[3][1][1]=="c")
512     {
513        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
514        list(list(l[3][size(l[3])]))+list(ideal(0));
515        l[2] = delete(l[2],size(l[2]));
516        l[3] = delete(l[3],size(l[3]));
517     }
518     else
519     {
520        l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
521        list(list(l[3][size(l[3])-1]))+list(ideal(0));
522        l[2] = delete(l[2],size(l[2]));
523        l[3] = delete(l[3],size(l[3])-1);
524     }
525
526     def S1 = ring(l);
527     setring S1;
528     number Num= number(imap(Ls,ff));
529     list l = ringlist(S1);
530     l[1][4][1] = Num;
531     S1 = ring(l);
532     setring S1;
533     ideal K = imap(Ls,K);
534     module Jp = imap(Ls,J);
535     def S2;
536     module Ip;
537     number Num;
538     /* ++++++ if the minpoly is irreducible then K = ideal(0) +++ */
539     if(size(K)==0)
540     {
541          module M = syz(Jp);
542          Ip = normalize(M);
543     }
544     else
545     {
546         for(j=1;j<=ncols(K);j++)
547         {
548              LL[j]=K[j];
549              Num = number(K[j]);
550              T2 = ringlist(S1);
551              T2[1][4][1] = Num;
552              S2 = ring(T2);
553              setring S2;
554              module M = syz(imap(Ls,J));
555              setring S1;
556              Lk[j]= imap(S2,M);
557         }
558         if(check_leadmonom_and_size(Lk))
559         {
560              // apply CRT for polynomials
561              setring Ls;
562              II =chinrempoly(imap(S1,Lk),imap(S1,LL));
563              setring S1;
564              Ip = normalize(imap(Ls,II));
565         }
566         else
567         {
568              setring S1;
569              Ip=[0];
570         }
571     }
572     setring S1;
573     module Gp = imap(br, result);
574     // test if Ip is contained in Gp
575     attrib(Gp, "isSB", 1);
576     for (i = ncols(Ip); i > 0; i--)
577     {
578          if (reduce(Ip[i], Gp, 1) != 0)
579          {
580              setring(br);
581              return(0);
582          }
583     }
584     // test if Gp is contained in syz(Jp)
585     if(size(module(matrix(Jp)*matrix(Gp)))!=0)
586     {
587         setring br;
588         return(0);
589     }
590     setring br;
591     return(1);
592}
593
594 ////////////////////////////////////////////////////////////////////////////////
595
596 static proc cleardenomIdeal(def I)
597 {
598     int t=ncols(I);
599     if(size(I)==0)
600     {
601         return(I);
602     }
603     else
604     {
605         for(int i=1;i<=t;i++)
606         {
607             I[i]=cleardenom(I[i]);
608         }
609     }
610     return(I);
611 }
612
613////////////////////////////////////////////////////////////////////////////////
614
615 static proc modStdparallelized_syzSB(module I, list #)
616 {
617     /* save options */
618     intvec opt = option(get);
619     option(redSB);
620     option(returnSB); // forces the procedure to return std of syz(I)
621
622     // apply modular command from modular.lib
623     if(size(#)>0)
624     {
625         I = modular("syz", list(I), testPrime, Modstd::deleteUnluckyPrimes_std,
626             pTest_syzmod, final_test, 536870909);
627     }
628     else
629     {
630         I = modular("Nfmodsyz::LiftPolyCRT_syz", list(I), PrimeTestTask_syz,
631             Modstd::deleteUnluckyPrimes_std,pTest_syz, final_Test_syz,536870909);
632     }
633     attrib(I, "isSB", 1);
634     option(set,opt);
635     return(I);
636 }
637
638////////////////////////////////////////////////////////////////////////////////
639/* main procedure */
640proc nfmodSyz(def I)
641"USAGE:  nfmodSyz(I); I ideal or module
642RETURN:  syzygy module of I over an algebraic number field
643SEE ALSO: syz
644EXAMPLE: example nfmodSyz; shows an example
645"
646{
647     if(typeof(I)!="ideal" and typeof(I)!="module")
648     {
649        ERROR("type of input must be either ideal or module");
650     }
651     else
652     {
653        module F = I;
654        kill I;
655        module I = F;
656     }
657     def Rbs=basering;
658     poly f;
659     int n=nvars(Rbs);
660     if(size(I)==0)
661     {
662         return(module([0]));
663     }
664     if(npars(Rbs)==0)
665     {
666        module M = modStdparallelized_syzSB(I,1); //if algebraic number is in Q
667        return(M);
668     }
669
670     def S;
671     list rl=ringlist(Rbs);
672     f=rl[1][4][1];
673
674     if(rl[3][1][1]!="c")
675     {
676        rl[2] = rl[2] + rl[1][2];
677        rl[3] = insert(rl[3], rl[1][3][1],1);
678        rl[1] = rl[1][1];
679     }
680     else
681     {
682        rl[2] = rl[2] + rl[1][2];
683        rl[3][size(rl[3])+1] = rl[1][3][1];
684        rl[1] = rl[1][1];
685     }
686
687     S = ring(rl);
688     setring S;
689     poly f=imap(Rbs,f);
690     def I=imap(Rbs,I);
691     I = simplify(I,2); // eraze the zero generatos
692     if(f==0)
693     {
694         ERROR("minpoly must be non-zero");
695     }
696     I=I,f;
697     def J_I = modStdparallelized_syzSB(I);
698     setring Rbs;
699     def J=imap(S,J_I);
700     J=simplify(J,2);
701     return(J);
702}
703example
704{ "EXAMPLE:"; echo = 2;
705    ring r1 =(0,a),(x,y),(c,dp);
706    minpoly = (a^3+2a+7);
707    module M1 = [(a/2+1)*y, 3*x-a*y],
708                [y-x,y2],
709                [x2-xy, ax-y];
710    M1 = nfmodSyz(M1);
711    M1;
712    ring r2 = (0,a),(x,y,z),(dp,c);
713    minpoly = (a3+a+1);
714    module M2 = [x2z+x+(-a)*y,z2+(a+2)*x],
715                [y2+(a)*z+(a),(a+3)*z3+(-a)*x2],
716                [-xz+(a2+3)*yz,xy+(a2)*z];
717    M2 = nfmodSyz(M2);
718    M2;
719    ring r3=0,(x,y),dp; // ring without parameter
720    module M3 = [x2 + y, xy], [-7y, 2x], [x2-y, 0];
721    M3=nfmodSyz(M3);
722    M3;
723    ring r4=0,(x,y),(c,dp); // ring without parameter
724    module M4 = [xy, x-y],
725                [x2 + y, 5y],
726                [- 7y, 2x],
727                [x2-y, 0];
728    M4 = nfmodSyz(M4);
729    M4;
730}
731
732/*
733////////////////////////////////////////////////////////////////////////////////
734
735//-------------nfmodSyz--old--version-----------------------------------
736
737static proc ret_indx(matrix M, int rp)
738{
739     int i,j,sc;
740     for(j=1;j<=ncols(M);j++)
741     {
742         sc=re_indx(M, rp, j,sc);
743         if(j>sc)
744         {
745             return(sc);
746             break;
747         }
748     }
749}
750
751////////////////////////////////////////////////////////////////////////////////
752
753static proc re_indx(matrix M, int rp, int j, int sc)
754{
755     int i;
756     if(M[1,j]!=0)
757     {
758         return(sc);
759     }
760     else
761     {
762          for(i=2;i<=rp; i++)
763          {
764             if(M[i,j]!=0)
765             {
766                 return(sc);
767                 break;
768             }
769          }
770          return(j);
771     }
772}
773
774////////////////////////////////////////////////////////////////////////////////
775
776proc syzMod(def I, list #)
777{
778     // compute the syzygy module of a given submodule
779
780     int i,j, s1,s2,s,k1,k2, i1,i2,i_I,j_I,i_3, tmp;
781     if(size(#)>0)
782     {
783         tmp = #[1]; // in this case, the procedure returns the reduced Groebner basis of
784                     // I which may partitioned to 2x2 block matrix
785     }
786     module sy;
787     module N = transpose(I),freemodule(ncols(I));
788     N= transpose(N);
789     matrix M = nfmodStd(N); // nfmodStd
790     i1 = ncols(M);
791     i2 = nrows(M);
792     i_I = ncols(I);
793     j_I = nrows(I);
794     i_3 = i2-i_I; // rows of zeros
795     j = ret_indx(M,j_I);
796     i = i_I; // i is nrows of syz module
797     if(j==0)
798     {
799         i = 0;
800     }
801
802     if(i==0)
803     {
804         if(tmp)
805         {
806             matrix T[i_I][i1] = M[(i_3+1)..i2,1..i1];
807             matrix G[i_3][i1] = M[1..i_3,1..i1];
808             list L = module([0]), module(G), T;
809             return(L);
810         }
811         return(module([0]));
812     }
813
814     if(!tmp)
815     {
816         matrix V[i][j] = M[(i_3+1)..i2,1..j];
817         sy = V;
818         return(sy);
819     }
820     else
821     {
822             matrix V[i][j] = M[(i_3+1)..i2,1..j];
823             matrix T[i_I][i1-j] = M[(i_3+1)..i2,(j+1)..i1];
824             matrix G[i_3][i1-j] = M[1..i_3,(j+1)..i1];
825             list L = module(V), module(G), T;
826             return(L);
827     }
828}
829
Note: See TracBrowser for help on using the repository browser.