source: git/Singular/LIB/nfmodstd.lib @ 320644

spielwiese
Last change on this file since 320644 was 320644, checked in by Hans Schoenemann <hannes@…>, 9 years ago
Merge branch 'algemodstd' of https://github.com/tigistdereje/Sources into tigistdereje-algemodstd Conflicts: doc/singular.doc
  • Property mode set to 100644
File size: 12.2 KB
Line 
1////////////////////////////////////////////////////////////////////////////////
2version="version nfmodstd.lib 4.0.1.0 Sep_2014 ";  // $Id$
3category="Commutative Algebra";
4info="
5
6LIBRARY:   nfmodstd.lib  Groebner bases of ideals in polynomial rings
7                           over 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 Groebner basis of an ideal in the polynomial
14  ring over an algebraic number field Q(t) using the modular methods, where t is
15  algebraic over the field of rational numbers Q. For the case Q(t) = Q, the procedure
16  is inspired by Arnold [1]. This idea is then extended
17  to the case t not in Q using factorization as follows:
18  Let f be the minimal polynomial of t.
19  For I, I' ideals in Q(t)[X], Q[X,t]/<f> respectively, we map I to I' via the map sending
20  t to t + <f>.
21  We first choose a prime p such that f has at least two factors in characteristic p and
22  add each factor f_i to I' to obtain the ideal J'_i = I' + <f_i>.
23  We then compute a standard basis G'_i of J'_i for each i and combine the G'_i to G_p
24  (a standard basis of I'_p) using chinese remaindering for polynomials. The procedure is
25  repeated for many primes p, where we compute the G_p in parallel until the
26  number of primes is sufficiently large to recover the correct standard basis G' of I'.
27  Finally, by mapping G' back to Q(t)[X], a standard basis G of I is obtained.
28
29REFERENCES:
30  [1] E. A. Arnold: Modular algorithms for computing Groebner bases.
31      J. Symb. Comp. 35, 403-419 (2003).
32
33PROCEDURES:
34  chinrempoly(l,m);       chinese remaindering for polynomials
35  nfmodStd(I);          standard basis of I over algebraic number field using modular methods
36";
37
38LIB "modstd.lib";
39
40////////////////////////////////////////////////////////////////////////////////
41
42static proc testPrime(int p, ideal I)
43{
44    /*
45     * test whether a prime p divides the denominator(s)
46     * and leading coefficients of generating set of ideal
47     */
48    int i,j;
49    poly f;
50    number num;
51    bigint d1,d2,d3;
52    for(i = 1; i <= size(I); i++)
53    {
54        f = cleardenom(I[i]);
55        if(f == 0)
56        {
57            return(0);
58        }
59        num = leadcoef(I[i])/leadcoef(f);
60        d1 = bigint(numerator(num));
61        d2 = bigint(denominator(num));
62        if( (d1 mod p) == 0)
63        {
64            return(0);
65        }
66        if((d2 mod p) == 0)
67        {
68            return(0);
69        }
70        for(j = size(f); j > 0; j--)
71        {
72            d3 = bigint(leadcoef(f[j]));
73            if( (d3 mod p) == 0)
74            {
75                return(0);
76            }
77        }
78    }
79    return(1);
80}
81
82////////////////////////////////////////////////////////////////////////////////
83/* return 1 if the number of factors are in the required bound , 0 else */
84
85static proc minpolyTask(poly f,int p)
86{
87    /*
88     * bound for irreducible factor(s) of (f mod p)
89     * see testfact()
90     */
91    int nr,k,ur;
92    ur=deg(f);
93    list L=factmodp(f,p);
94    if(degtest(L[2])==1)
95    {
96        // now each factor is squarefree
97        if(ur<=3)
98        {
99            return(1);
100        }
101        else
102        {
103            nr = testfact(ur);
104            k=ncols(L[1]);
105            if(nr < k && k < (ur-nr))// set bound for k
106            {
107                return(1);
108            }
109        }
110    }
111    return(0);
112}
113
114////////////////////////////////////////////////////////////////////////////////
115/* return 1 if both testPrime(p,J) and minpolyTask(f,p) is true, 0 else */
116
117static proc PrimeTestTask(int p, list L)
118{
119    /* L=list(I), I=J,f; J ideal , f minpoly */
120    int sz,nr,dg;
121    ideal J=L[1];
122    sz=ncols(J);
123    poly f=J[sz];
124    dg=deg(f);
125    if(!testPrime(p,J) or !minpolyTask(f,p))
126    {
127        return(0);
128    }
129    return(1);
130}
131
132////////////////////////////////////////////////////////////////////////////////
133/* compute factors of f mod p with multiplicity */
134
135static proc factmodp(poly f, int p)
136{
137    def R=basering;
138    list l=ringlist(R);
139    l[1]=p;
140    def S=ring(l);
141    setring S;
142    list L=factorize(imap(R,f),2);
143    ideal J=L[1];
144    intvec v=L[2];
145    list scx=J,v;
146    setring R;
147    return(imap(S,scx));
148    kill S;
149}
150
151////////////////////////////////////////////////////////////////////////////////
152/* set a bound for number of factors w.r.t degree nr*/
153
154static proc testfact(int nr)
155{
156    // nr must be greater than 3
157    int i;
158    if(nr>3 and nr<=5)
159    {
160        i=1;
161    }
162    if(nr>5 and nr<=10)
163    {
164        i=2;
165    }
166    if(nr>10 and nr<=15)
167    {
168        i=3;
169    }
170    if(nr>15 and nr<=20)
171    {
172        i=4;
173    }
174    if(nr>20 and nr<=25)
175    {
176        i=5;
177    }
178    if(nr>25 and nr<=30)
179    {
180        i=6;
181    }
182    if(nr>30)
183    {
184        i=10;
185    }
186    return(i);
187}
188
189///////////////////////////////////////////////////////////////////////////////
190// return 1 if v[i]>1 , 0 else
191
192static proc degtest(intvec v)
193{
194    for(int j=1;j<=nrows(v);j++)
195    {
196        if(v[j]>1)
197        {
198            return(0);
199        }
200    }
201    return(1);
202}
203
204////////////////////////////////////////////////////////////////////////////////
205
206static proc chinRm(list m, list ll, list lk,list l1,int uz)
207{
208    poly ff,c;
209
210    for(int i=1;i<=uz;i++)
211    {
212        c = division(l1[i]*ll[i],m[i])[2][1];
213        ff = ff + c*lk[i];
214    }
215    return(ff);
216}
217
218////////////////////////////////////////////////////////////////////////////////
219
220proc chinrempoly(list l,list m)
221"USAGE:  chinrempoly(l, m); l list, m list
222RETURN:  a polynomial (resp. ideal) which is congruent to l[i] modulo m[i] for all i
223NOTE: The procedure applies chinese remaindering to the first argument w.r.t. the
224      moduli given in the second. The elements in the first list must be of same type
225      which can be polynomial or ideal. The moduli must be of type polynomial. Elements
226      in the second list must be distinct and co-prime.
227SEE ALSO: chinrem
228EXAMPLE: example chinrempoly; shows an example
229"
230{
231    int i,j,sz,uz;
232    uz = size(l);
233    sz = ncols(ideal(l[1]));
234    poly f=1;
235    for(i=1;i<=uz;i++)
236    {
237        f=f*m[i];
238    }
239
240    ideal I,J;
241    list l1,ll,lk,l2;
242    poly c,ff;
243    for(j=1;j<=uz;j++)
244    {
245        lk[j]=f/m[j];
246        ll[j]=extgcd(lk[j],m[j])[2];
247    }
248
249    for(i=1;i<=sz;i++)
250    {
251        for(j=1;j<=uz;j++)
252        {
253            I = l[j];
254            l1[j] = I[i];
255        }
256        J[i] = chinRm(m,ll,lk,l1,uz);
257    }
258    return(J);
259}
260example
261{ "EXAMPLE:"; echo = 2;
262    ring rr=97,x,dp;
263    poly f=x^7-7*x + 3;
264    ideal J=factorize(f,1);
265    J;
266    list m=J[1..ncols(J)];
267    list l= x^2+2*x+3, x^2+5, x^2+7;
268    ideal I=chinrempoly(l,m);
269    I;
270    ring s=0,x,dp;
271    list m= x^2+2*x+3, x^3+5, x^4+x^3+7;
272    list l=x^3 + 2, x^4 + 7, x^5 + 11;
273    ideal I=chinrempoly(l,m);
274    I;
275    int p=prime(536546513);
276    ring r = p, (x,y,a), (dp(2),dp(1));
277    poly minpolynomial = a^2+1;
278    ideal kf=factorize(minpolynomial,1);//return factors without multiplicity
279    kf;
280    ideal k=(a+1)*x2+y, 3x-ay+ a+2;
281    option(redSB);
282    ideal k1=k,kf[1];
283    ideal k2 =k,kf[2];
284    k1=std(k1);
285    k2=std(k2);
286    list l=k1,k2;
287    list m=kf[1..ncols(kf)];
288    ideal I=chinrempoly(l,m);
289    I=simplify(I,2);
290    I;
291}
292
293////////////////////////////////////////////////////////////////////////////////
294
295static proc check_leadmonom_and_size(list L)
296{
297    /*
298     * compare the size of ideals in the list and
299     * check the corresponding leading monomials
300     * size(L)>=2
301     */
302    ideal J=L[1];
303    int i=size(L);
304    int sc=ncols(J);
305    int j,k;
306    poly g=leadmonom(J[1]);
307    for(j=1;j<=i;j++)
308    {
309        if(ncols(L[j])!=sc)
310        {
311            return(0);
312        }
313    }
314    for(k=2;k<=i;k++)
315    {
316        for(j=1;j<=sc;j++)
317        {
318            if(leadmonom(J[j])!=leadmonom(L[k][j]))
319            {
320                return(0);
321            }
322        }
323    }
324    return(1);
325}
326
327////////////////////////////////////////////////////////////////////////////////
328
329static proc LiftPolyCRT(ideal I)
330{
331    /*
332     * compute std for each factor and combine this result
333     * to modulo minpoly via CRT for poly over char p>0
334     */
335    int u,in,j;
336    list LL,Lk;
337    ideal J,K,II;
338    poly f;
339    u=ncols(I);
340    J=I[1..u-1];
341    f=I[u];
342    K=factorize(f,1);
343    in=ncols(K);
344    for(j=1;j<=in;j++)
345    {
346        LL[j]=K[j];
347        ideal I(j)=J,K[j];
348        I(j)=std(I(j));
349        if(size(I(j))==1)
350        {
351            Lk[j]=I(j);
352        }
353        else
354        {
355            I(j)[1]=0;
356            I(j)=simplify(I(j), 2);
357            Lk[j]=I(j);
358        }
359    }
360    if(check_leadmonom_and_size(Lk))
361    {
362        // apply CRT for polynomials
363        II =chinrempoly(Lk,LL),f;
364    }
365    else
366    {
367        II=0;
368    }
369    return(II);
370}
371
372////////////////////////////////////////////////////////////////////////////////
373
374static proc PtestStd(string command, list args, ideal result, int p)
375{
376    /*
377     * let G be std of I which is not yet known whether it is the correct
378     *  standard basis or not. So this procedure does the first test
379     */
380    def br = basering;
381    list lbr = ringlist(br);
382    if (typeof(lbr[1]) == "int")
383    {
384        lbr[1] = p;
385    }
386    else
387    {
388        lbr[1][1] = p;
389    }
390    def rp = ring(lbr);
391    setring(rp);
392    ideal Ip = fetch(br, args)[1];
393    ideal Gp = fetch(br, result);
394    attrib(Gp, "isSB", 1);
395    int i;
396    for (i = ncols(Ip); i > 0; i--)
397    {
398        if (reduce(Ip[i], Gp, 1) != 0)
399        {
400            setring(br);
401            return(0);
402        }
403    }
404    Ip = LiftPolyCRT(Ip);
405    attrib(Ip,"isSB",1);
406    for (i = ncols(Gp); i > 0; i--)
407    {
408        if (reduce(Gp[i], Ip, 1) != 0)
409        {
410            setring(br);
411            return(0);
412        }
413    }
414    setring(br);
415    return(1);
416}
417
418////////////////////////////////////////////////////////////////////////////////
419
420static proc cleardenomIdeal(ideal I)
421{
422    int t=ncols(I);
423    if(size(I)==0)
424    {
425        return(I);
426    }
427    else
428    {
429        for(int i=1;i<=t;i++)
430        {
431            I[i]=cleardenom(I[i]);
432        }
433    }
434    return(I);
435}
436
437////////////////////////////////////////////////////////////////////////////////
438
439static proc modStdparallelized(ideal I)
440{
441    // apply modular.lib
442    /* save options */
443    intvec opt = option(get);
444    option(redSB);
445    I = modular("Nfmodstd::LiftPolyCRT", list(I), PrimeTestTask, Modstd::deleteUnluckyPrimes_std,
446              PtestStd, Modstd::finalTest_std,536870909);
447    attrib(I, "isSB", 1);
448    option(set,opt);
449    return(I);
450}
451
452////////////////////////////////////////////////////////////////////////////////
453/* main procedure */
454proc nfmodStd(ideal I, list #)
455"USAGE:  nfmodStd(I, #); I ideal, # optional parameters
456RETURN:  standard basis of I over algebraic number field
457NOTE: The procedure passes to @ref{modStd} if the ground field has no
458      parameter. In this case, the optional parameters # (if given)
459      are directly passed to @ref{modStd}.
460SEE ALSO: modStd
461EXAMPLE: example nfmodStd; shows an example
462"
463{
464    list L=#;
465    def Rbs=basering;
466    poly f;
467    ideal J;
468    int n=nvars(Rbs);
469    if(size(I)==0)
470    {
471        return(ideal(0));
472    }
473    if(npars(Rbs)==0)
474    {
475        J=modStd(I,L);//if algebraic number is in Q
476        if(size(#)>0)
477        {
478            return(cleardenomIdeal(J));
479        }
480        return(J);
481    }
482    list rl=ringlist(Rbs);
483    f=rl[1][4][1];
484    rl[2][n+1]=rl[1][2][1];
485    rl[1]=rl[1][1];
486    rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
487    rl[3][size(rl[3])-1]=list("dp",1);
488    def S=ring(rl);
489    setring S;
490    poly f=imap(Rbs,f);
491    ideal I=imap(Rbs,I);
492    I = simplify(I,2);// eraze the zero generatos
493    ideal J;
494    if(f==0)
495    {
496        ERROR("minpoly must be non-zero");
497    }
498    I=I,f;
499    J=modStdparallelized(I);
500    setring Rbs;
501    J=imap(S,J);
502    J=simplify(J,2);
503    if(size(#)>0)
504    {
505        return(cleardenomIdeal(J));
506    }
507    return(J);
508}
509example
510{ "EXAMPLE:"; echo = 2;
511    ring r1 =(0,a),(x,y),dp;
512    minpoly =a^2+1;
513    ideal k=(a/2+1)*x^2+2/3y, 3*x-a*y+ a/7+2;
514    ideal I=nfmodStd(k);
515    I;
516    ring r2 =(0,a),(x,y,z),dp;
517    minpoly =a^3 +2;
518    ideal k=(a^2+a/2)*x^2+(a^2 -2/3*a)*yz, (3*a^2+1)*zx-(a+4/7)*y+ a+2/5;
519    ideal IJ=nfmodStd(k);
520    IJ;
521    ring r3=0,(x,y),dp;// ring without parameter
522    ideal I = x2 + y, xy - 7y + 2x;
523    I=nfmodStd(I);
524    I;
525}
Note: See TracBrowser for help on using the repository browser.