Changeset 204cf0 in git


Ignore:
Timestamp:
Apr 13, 2016, 8:26:20 PM (7 years ago)
Author:
dere <dkifle16@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
5a4188b414ab493b3e3e697845d6e6e6f3e783e8
Parents:
b1cef7673ad345288c0f013557cd800eaed8d1aa
Message:
add: new library nfmodsyz.lib
Files:
7 added
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/nfmodstd.lib

    rb1cef7 r204cf0  
    1313  A library for computing the Groebner basis of an ideal in the polynomial
    1414  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:
     15  algebraic over the field of rational numbers Q. For the case Q(t) = Q, the
     16  procedure is inspired by Arnold [1]. This idea is then extended to the case
     17  t not in Q using factorization as follows:
    1818  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.
     19  For I, I' ideals in Q(t)[X], Q[X,t]/<f> respectively, we map I to I' via the
     20  map sending  t to t + <f>. We first choose a prime p such that f has at least
     21  two factors in characteristic p and add each factor f_i to I' to obtain the
     22  ideal J'_i = I' + <f_i>. We then compute a standard basis G'_i of J'_i for each
     23  i and combine the G'_i to G_p (a standard basis of I'_p) using chinese remaindering
     24  for polynomials. The procedure is repeated for many primes p, where we compute the
     25  G_p in parallel until the number of primes is sufficiently large to recover the
     26  correct standard basis G' of I'. Finally, by mapping G' back to Q(t)[X], a standard
     27  basis G of I is obtained.
     28@*The procedure also works if the input is a module. For this, we consider the
     29  rings A = Q(t)[X] and A' = (Q[t]/<f>)[X]. For submodules I, I' in A^m, A'^m,
     30  respectively, we map I to I' via the map sending t to t + <f>. As above, we first
     31  choose a prime p such that f has at least two factors in characteristic p. For each
     32  factor f_{i,p} of f_p := (f mod p), we set I'_{i,p} := (I'_p mod f_{i,p}). We then
     33  compute a standard basis G'_i of I'_{i,p} over F_p[t]/<f_{i,p}> for each i and
     34  combine the G'_i to G_p (a standard basis of I'_p) using chinese remaindering for
     35  polynomials. The procedure is repeated for many primes p as described above and we
     36  finally obtain a standard basis of I.
    2837
    2938REFERENCES:
    … …  
    3342PROCEDURES:
    3443  chinrempoly(l,m);       chinese remaindering for polynomials
    35   nfmodStd(I);          standard basis of I over algebraic number field using modular methods
     44  nfmodStd(I);          standard basis of I over algebraic number field using modular
     45                        methods
    3646";
    3747
    … …  
    4050////////////////////////////////////////////////////////////////////////////////
    4151
    42 static proc testPrime(int p, ideal I)
     52static proc testPrime(int p, list L)
    4353{
    4454    /*
    … …  
    4656     * and leading coefficients of generating set of ideal
    4757     */
    48     int i,j;
    49     poly f;
     58    int i,j,k, tmp;
     59    def f;
     60    def I = L[1]; // list L = def I
    5061    number num;
    5162    bigint d1,d2,d3;
    52     for(i = 1; i <= size(I); i++)
     63    if(typeof(I)=="ideal")
     64    {
     65        tmp=1;
     66    }
     67    for(i = 1; i <= ncols(I); i++)
    5368    {
    5469        f = cleardenom(I[i]);
    5570        if(f == 0)
    5671        {
    57             return(0);
     72             return(0);
    5873        }
    5974        num = leadcoef(I[i])/leadcoef(f);
    … …  
    6883            return(0);
    6984        }
    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         }
     85        if(tmp)
     86        {
     87            for(j = size(f); j > 0; j--)
     88            {
     89                d3 = bigint(leadcoef(f[j]));
     90                if( (d3 mod p) == 0)
     91                {
     92                    return(0);
     93                }
     94            }
     95        }
     96        else
     97        {
     98            for(j = nrows(f); j > 0; j--)
     99            {
     100                for(k=1;k<=size(f[j]);k++)
     101                {
     102                    d3 = bigint(leadcoef(f[j][k]));
     103                    if((d3 mod p) == 0)
     104                    {
     105                        return(0);
     106                    }
     107                }
     108            }
     109        }
     110
    78111    }
    79112    return(1);
    … …  
    103136            nr = testfact(ur);
    104137            k=ncols(L[1]);
    105             if(nr < k && k < (ur-nr))// set bound for k
     138            if(nr < k && k < (ur-nr))// set a bound for k
    106139            {
    107140                return(1);
    … …  
    118151{
    119152    /* L=list(I), I=J,f; J ideal , f minpoly */
    120     int sz,nr,dg;
    121     ideal J=L[1];
     153    int sz,nr;
     154    def J=L[1];
    122155    sz=ncols(J);
    123     poly f=J[sz];
    124     dg=deg(f);
    125     if(!testPrime(p,J) or !minpolyTask(f,p))
     156    def f=J[sz];
     157    poly g;
     158    if(typeof(f)=="vector")
     159    {
     160        g = f[1];
     161    }
     162    else
     163    {
     164        g = f;
     165    }
     166    if(!testPrime(p,list(J)) or !minpolyTask(g,p))
    126167    {
    127168        return(0);
    … …  
    206247static proc chinRm(list m, list ll, list lk,list l1,int uz)
    207248{
    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);
     249    if(typeof(l1[1])=="ideal" or typeof(l1[1])=="poly")
     250    {
     251        poly ff,c;
     252        for(int i=1;i<=uz;i++)
     253        {
     254            c = division(l1[i]*ll[i],m[i])[2][1];
     255            ff = ff + c*lk[i];
     256        }
     257        return(ff);
     258    }
     259    else
     260    {
     261        vector ff,c;
     262        for(int i=1;i<=uz;i++)
     263        {
     264            c = vector(m[i]);
     265            attrib(c,"isSB",1);
     266            ff = ff + (reduce(l1[i]*ll[i],c))*lk[i];
     267        }
     268        return(ff);
     269    }
    216270}
    217271
    … …  
    220274proc chinrempoly(list l,list m)
    221275"USAGE:  chinrempoly(l, m); l list, m list
    222 RETURN:  a polynomial (resp. ideal) which is congruent to l[i] modulo m[i] for all i
     276RETURN:  a polynomial (resp. ideal/module) which is congruent to l[i] modulo m[i]
     277         for all i
    223278NOTE: 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.
     279      moduli given in the second. The elements in the first list must be of the same
     280      type which can be polynomial, ideal, or module. The moduli must be of type
     281      polynomial. The elements in the second list must be distinct and co-prime.
    227282SEE ALSO: chinrem
    228283EXAMPLE: example chinrempoly; shows an example
    229284"
    230285{
    231     int i,j,sz,uz;
     286    int i,j,sz,uz, tmp;
    232287    uz = size(l);
    233     sz = ncols(ideal(l[1]));
     288    if(typeof(l[1])=="ideal" or typeof(l[1])=="poly")
     289    {
     290        sz = ncols(ideal(l[1]));
     291        tmp = 1;
     292    }
     293    else
     294    {
     295        sz = ncols(module(l[1]));
     296    }
    234297    poly f=1;
    235298    for(i=1;i<=uz;i++)
    … …  
    238301    }
    239302
    240     ideal I,J;
    241303    list l1,ll,lk,l2;
    242304    poly c,ff;
    … …  
    247309    }
    248310
    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);
     311    if(tmp)
     312    {
     313        ideal I,J;
     314        for(i=1;i<=sz;i++)
     315        {
     316            for(j=1;j<=uz;j++)
     317            {
     318                I = l[j];
     319                l1[j] = I[i];
     320            }
     321            J[i] = chinRm(m,ll,lk,l1,uz);
     322        }
     323        return(J);
     324    }
     325    else
     326    {
     327        module I,J;
     328        for(i=1;i<=sz;i++)
     329        {
     330            for(j=1;j<=uz;j++)
     331            {
     332                I = l[j];
     333                l1[j] = I[i];
     334            }
     335            J[i] = chinRm(m,ll,lk,l1,uz);
     336        }
     337        return(J);
     338    }
    259339}
    260340example
    … …  
    276356    ring r = p, (x,y,a), (dp(2),dp(1));
    277357    poly minpolynomial = a^2+1;
    278     ideal kf=factorize(minpolynomial,1);//return factors without multiplicity
     358    ideal kf=factorize(minpolynomial,1); //return factors without multiplicity
    279359    kf;
    280360    ideal k=(a+1)*x2+y, 3x-ay+ a+2;
    … …  
    289369    I=simplify(I,2);
    290370    I;
    291 }
    292 
     371    l = module(k1[2..ncols(k1)]), module(k2[2..ncols(k2)]);
     372    module M = chinrempoly(l,m);
     373    M;
     374}
    293375////////////////////////////////////////////////////////////////////////////////
    294376
    … …  
    300382     * size(L)>=2
    301383     */
    302     ideal J=L[1];
     384    def J=L[1];
    303385    int i=size(L);
    304386    int sc=ncols(J);
    305387    int j,k;
    306     poly g=leadmonom(J[1]);
     388    def g=leadmonom(J[1]);
    307389    for(j=1;j<=i;j++)
    308390    {
    … …  
    327409////////////////////////////////////////////////////////////////////////////////
    328410
    329 static proc LiftPolyCRT(ideal I)
     411static proc LiftPolyCRT(def I)
    330412{
    331413    /*
    … …  
    333415     * to modulo minpoly via CRT for poly over char p>0
    334416     */
     417    def sl;
    335418    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);
     419    list LL,Lk,T2;
     420    if(typeof(I)=="ideal")
     421    {
     422
     423        ideal J,K,II;
     424        poly f;
     425        u=ncols(I);
     426        J=I[1..u-1];
     427        f=I[u];
     428        K=factorize(f,1);
     429        in=ncols(K);
     430        for(j=1;j<=in;j++)
     431        {
     432            LL[j]=K[j];
     433            ideal I(j)=J,K[j];
     434            I(j)=std(I(j));
     435            if(size(I(j))==1)
     436            {
     437                Lk[j]=I(j);
     438            }
     439            else
     440            {
     441                I(j)[1]=0;
     442                I(j)=simplify(I(j), 2);
     443                Lk[j]=I(j);
     444            }
     445        }
     446        if(check_leadmonom_and_size(Lk))
     447        {
     448            // apply CRT for polynomials
     449            II =chinrempoly(Lk,LL),f;
    352450        }
    353451        else
    354452        {
    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 
    374 static proc PtestStd(string command, list args, ideal result, int p)
     453            II=0;
     454        }
     455        return(II);
     456    }
     457    else
     458    {
     459        module J,II;
     460        vector f;
     461        u=ncols(I);
     462        J=I[1..u-1];
     463        f=I[u];
     464        poly ff = f[1];
     465        ideal K=factorize(ff,1);
     466        in=ncols(K);
     467        def Ls = basering;
     468        list l = ringlist(Ls);
     469        if(l[3][1][1]=="c")
     470        {
     471             l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
     472             list(list(l[3][size(l[3])]))+list(ideal(0));
     473             l[2] = delete(l[2],size(l[2]));
     474             l[3] = delete(l[3],size(l[3]));
     475        }
     476        else
     477        {
     478            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
     479            list(list(l[3][size(l[3])-1]))+list(ideal(0));
     480            l[2] = delete(l[2],size(l[2]));
     481            l[3] = delete(l[3],size(l[3])-1);
     482        }
     483
     484        def S1 = ring(l);
     485        setring S1;
     486        number Num= number(imap(Ls,ff));
     487        list l = ringlist(S1);
     488        l[1][4][1] = Num;
     489        S1 = ring(l);
     490        setring S1;
     491        ideal K = imap(Ls,K);
     492        def S2;
     493        module II;
     494        number Num;
     495        /* ++++++ if minpoly is irreducible then K will be the zero ideal +++ */
     496        if(size(K)==0)
     497        {
     498            module M = std(imap(Ls,J));
     499            if(size(M)==1 && M[1]==[1])
     500            {
     501                setring Ls;
     502                return(module([1]));
     503            }
     504            II = normalize(M);
     505        }
     506        else
     507        {
     508            for(j=1;j<=in;j++)
     509            {
     510                LL[j]=K[j];
     511                Num = number(K[j]);
     512                T2 = ringlist(S1);
     513                T2[1][4][1] = Num;
     514                S2 = ring(T2);
     515                setring S2;
     516                module M = std(imap(Ls,J));
     517                if(size(M)== 1 && M[1]==[1])
     518                {
     519                    setring Ls;
     520                    return(module([1]));
     521                }
     522                setring S1;
     523                Lk[j]= imap(S2,M);
     524            }
     525
     526            if(check_leadmonom_and_size(Lk))
     527            {
     528                // apply CRT for polynomials
     529                setring Ls;
     530                II =chinrempoly(imap(S1,Lk),imap(S1,LL));
     531                setring S1;
     532                II = normalize(imap(Ls,II));
     533            }
     534            else
     535            {
     536                setring S1;
     537                II=[0];
     538            }
     539        }
     540        setring Ls;
     541        return(imap(S1,II));
     542    }
     543}
     544
     545////////////////////////////////////////////////////////////////////////////////
     546
     547/* test if 'result' is a GB of the input ideal */
     548static proc final_Test_minpolyzero(string command, alias list args, module result)
     549{
     550    int i;
     551    list arg = args;
     552    attrib(result, "isSB", 1);
     553    for (i = ncols(args[1]); i > 0; i--)
     554    {
     555         if (reduce(args[1][i], result, 1) != 0)
     556         {
     557              return(0);
     558         }
     559    }
     560    /* test if result is a GB */
     561    module G = std(result);
     562    if (size(reduce(G, result,1))!=0) //Modstd::reduce_parallel(G, result)
     563    {
     564         return(0);
     565    }
     566    return(1);
     567}
     568
     569////////////////////////////////////////////////////////////////////////////////
     570
     571/* test if 'result' is a GB of the input ideal */
     572static proc final_Test(string command, alias list args, def result)
     573{
     574    int i;
     575    list arg = args;
     576    // modified for module case
     577    if(typeof(args[1])=="ideal")
     578    {
     579        /* test if args[1] is in result */
     580        attrib(result, "isSB", 1);
     581        for (i = ncols(args[1]); i > 0; i--)
     582        {
     583            if (reduce(args[1][i], result, 1) != 0)
     584            {
     585                return(0);
     586            }
     587        }
     588
     589        /* test if result is a GB */
     590        ideal G = std(result);
     591        if (size(reduce(G, result,1))!=0) //Modstd::reduce_parallel(G, result)
     592        {
     593            return(0);
     594        }
     595        return(1);
     596    }
     597    else
     598    {
     599        /* test if args[1] is in result */
     600        def Ls = basering;
     601        list l = ringlist(Ls);
     602        if(l[3][1][1]=="c")
     603        {
     604             l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
     605             list(list(l[3][size(l[3])]))+list(ideal(0));
     606             l[2] = delete(l[2],size(l[2]));
     607             l[3] = delete(l[3],size(l[3]));
     608        }
     609        else
     610        {
     611            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
     612            list(list(l[3][size(l[3])-1]))+list(ideal(0));
     613            l[2] = delete(l[2],size(l[2]));
     614            l[3] = delete(l[3],size(l[3])-1);
     615        }
     616        def sL = ring(l);
     617        kill l;
     618        setring sL;
     619        list arg = imap(Ls,arg);
     620        module arg_s =arg[1];
     621        list l = ringlist(sL);
     622        l[1][4][1] = arg_s[ncols(arg_s)][1];
     623        arg_s = arg_s[1..ncols(arg_s)-1];
     624        def cL = ring(l);
     625        setring cL;
     626        module ar_gs = imap(sL,arg_s);
     627        def Result = imap(Ls,result);
     628        attrib(Result, "isSB", 1);
     629        for (i = ncols(ar_gs); i > 0; i--)
     630        {
     631            if (reduce(ar_gs[i], Result, 1) != 0)
     632            {
     633                setring Ls;
     634                return(0);
     635            }
     636        }
     637        // test if result is a GB
     638        module G = std(Result);
     639        if (size(reduce(G,Result,1))!=0) //Modstd::reduce_parallel(G, Result)
     640        {
     641            setring Ls;
     642            return(0);
     643        }
     644        setring Ls;
     645        return(1);
     646    }
     647}
     648
     649////////////////////////////////////////////////////////////////////////////////
     650
     651static proc PtestStd_minpolyzero(string command, list args, module result, int p)
    375652{
    376653    /*
    377654     * 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
     655     *  standard basis. So this procedure does the first test
    379656     */
    380657    def br = basering;
    … …  
    390667    def rp = ring(lbr);
    391668    setring(rp);
    392     ideal Ip = fetch(br, args)[1];
    393     ideal Gp = fetch(br, result);
     669    module Ip = imap(br, args)[1];
     670    int i;
     671    module Gp = imap(br, result);
    394672    attrib(Gp, "isSB", 1);
    395     int i;
    396673    for (i = ncols(Ip); i > 0; i--)
    397674    {
    … …  
    402679        }
    403680    }
    404     Ip = LiftPolyCRT(Ip);
     681    Ip = std(Ip);
    405682    attrib(Ip,"isSB",1);
    406683    for (i = ncols(Gp); i > 0; i--)
    … …  
    418695////////////////////////////////////////////////////////////////////////////////
    419696
    420 static proc cleardenomIdeal(ideal I)
     697static proc PtestStd(string command, list args, def result, int p)
     698{
     699    /*
     700     * let G be std of I which is not yet known whether it is the correct
     701     *  standard basis. So this procedure does the first test
     702     */
     703    def br = basering;
     704
     705    list lbr = ringlist(br);
     706    if (typeof(lbr[1]) == "int")
     707    {
     708        lbr[1] = p;
     709    }
     710    else
     711    {
     712        lbr[1][1] = p;
     713    }
     714    def rp = ring(lbr);
     715    setring(rp);
     716    def Ip = imap(br, args)[1];
     717
     718    int u,in,j,i;
     719    list LL,Lk,T2;
     720    if(typeof(Ip)!="ideal")
     721    {
     722        module J,II;
     723        vector f;
     724        u=ncols(Ip);
     725        J=Ip[1..u-1];
     726        f=Ip[u];
     727        poly ff = f[1];
     728        ideal K=factorize(ff,1);
     729        in=ncols(K);
     730        def Ls = basering;
     731        list l = ringlist(Ls);
     732        if(l[3][1][1]=="c")
     733        {
     734            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
     735            list(list(l[3][size(l[3])]))+list(ideal(0));
     736            l[2] = delete(l[2],size(l[2]));
     737            l[3] = delete(l[3],size(l[3]));
     738        }
     739        else
     740        {
     741            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
     742            list(list(l[3][size(l[3])-1]))+list(ideal(0));
     743            l[2] = delete(l[2],size(l[2]));
     744            l[3] = delete(l[3],size(l[3])-1);
     745        }
     746        def S1 = ring(l);
     747        setring S1;
     748        number Num= number(imap(Ls,ff));
     749        list l = ringlist(S1);
     750        l[1][4][1] = Num;
     751        S1 = ring(l);
     752        setring S1;
     753        ideal K = imap(Ls,K);
     754        module Jp = imap(Ls,J);
     755        def S2;
     756        module Ip;
     757        number Num;
     758        /* ++++++ if the minpoly is irreducible then K = ideal(0) +++ */
     759        if(size(K)==0)
     760        {
     761            module M = std(Jp);
     762            Ip = normalize(M);
     763        }
     764        else
     765        {
     766            for(j=1;j<=ncols(K);j++)
     767            {
     768                LL[j]=K[j];
     769                Num = number(K[j]);
     770                T2 = ringlist(S1);
     771                T2[1][4][1] = Num;
     772                S2 = ring(T2);
     773                setring S2;
     774                module M = std(imap(Ls,J));
     775                setring S1;
     776                Lk[j]= imap(S2,M);
     777            }
     778            if(check_leadmonom_and_size(Lk))
     779            {
     780                // apply CRT for polynomials
     781                setring Ls;
     782                II =chinrempoly(imap(S1,Lk),imap(S1,LL));
     783                setring S1;
     784                Ip = normalize(imap(Ls,II));
     785            }
     786            else
     787            {
     788                setring S1;
     789                Ip=[0];
     790            }
     791        }
     792        setring S1;
     793        module Gp = imap(br, result);
     794        attrib(Gp, "isSB", 1);
     795        for (i = ncols(Jp); i > 0; i--)
     796        {
     797            if (reduce(Jp[i], Gp, 1) != 0)
     798            {
     799                setring(br);
     800                return(0);
     801            }
     802        }
     803
     804        attrib(Ip,"isSB",1);
     805        for (i = ncols(Gp); i > 0; i--)
     806        {
     807            if (reduce(Gp[i], Ip, 1) != 0)
     808            {
     809                setring(br);
     810                return(0);
     811            }
     812        }
     813        setring(br);
     814        return(1);
     815    }
     816    else
     817    {
     818        ideal Gp = imap(br, result);
     819        attrib(Gp, "isSB", 1);
     820        for (i = ncols(Ip); i > 0; i--)
     821        {
     822            if (reduce(Ip[i], Gp, 1) != 0)
     823            {
     824                setring(br);
     825                return(0);
     826            }
     827        }
     828        Ip = LiftPolyCRT(Ip);
     829        attrib(Ip,"isSB",1);
     830        for (i = ncols(Gp); i > 0; i--)
     831        {
     832            if (reduce(Gp[i], Ip, 1) != 0)
     833            {
     834                setring(br);
     835                return(0);
     836            }
     837        }
     838        setring(br);
     839        return(1);
     840    }
     841}
     842
     843////////////////////////////////////////////////////////////////////////////////
     844
     845static proc cleardenomIdeal(def I)
    421846{
    422847    int t=ncols(I);
    … …  
    437862////////////////////////////////////////////////////////////////////////////////
    438863
    439 static proc modStdparallelized(ideal I)
     864static proc modStdparallelized(def I, list #)
    440865{
    441866    // apply modular.lib
    … …  
    443868    intvec opt = option(get);
    444869    option(redSB);
    445     I = modular("Nfmodstd::LiftPolyCRT", list(I), PrimeTestTask, Modstd::deleteUnluckyPrimes_std,
    446               PtestStd, Modstd::finalTest_std,536870909);
     870    if(size(#)>0)
     871    {
     872        I = modular("std", list(I), testPrime, Modstd::deleteUnluckyPrimes_std,
     873                PtestStd_minpolyzero, final_Test_minpolyzero,536870909);
     874    }
     875    else
     876    {
     877        I = modular("Nfmodstd::LiftPolyCRT", list(I), PrimeTestTask,
     878                    Modstd::deleteUnluckyPrimes_std,PtestStd, final_Test,536870909);
     879    }
    447880    attrib(I, "isSB", 1);
    448881    option(set,opt);
    … …  
    452885////////////////////////////////////////////////////////////////////////////////
    453886/* main procedure */
    454 proc nfmodStd(ideal I, list #)
    455 "USAGE:  nfmodStd(I, #); I ideal, # optional parameters
     887proc nfmodStd(def I, list #)
     888"USAGE:  nfmodStd(I, #); I ideal or module, # optional parameters
    456889RETURN:  standard basis of I over algebraic number field
    457890NOTE: The procedure passes to @ref{modStd} if the ground field has no
    … …  
    462895"
    463896{
    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);
     897     list L=#;
     898     def Rbs=basering;
     899     poly f;
     900     int tmp=1;
     901     if(typeof(I)!="ideal")
     902     {
     903         tmp =0;
     904     }
     905     int n=nvars(Rbs);
     906     if(size(I)==0)
     907     {
     908         if(!tmp)
     909         {
     910             return(module([0]));
     911         }
     912         return(ideal(0));
     913     }
     914     if(npars(Rbs)==0)
     915     {
     916         //if algebraic number is in Q
     917         if(typeof(I)=="module")
     918         {
     919             def J = modStdparallelized(I,1);
     920             return(J);
     921         }
     922         else
     923         {
     924             def J=modStd(I,L);
     925             return(J);
     926         }
     927     }
     928     def S;
     929     list rl=ringlist(Rbs);
     930     f=rl[1][4][1];
     931     if(tmp)
     932     {
     933         rl[2][n+1]=rl[1][2][1];
     934         rl[1]=rl[1][1];
     935         rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
     936         rl[3][size(rl[3])-1]=list("dp",1);
     937     }
     938     else
     939     {
     940         if(rl[3][1][1]!="c")
     941         {
     942             rl[2] = rl[2] + rl[1][2];
     943             rl[3] = insert(rl[3], rl[1][3][1],1);
     944             rl[1] = rl[1][1];
     945         }
     946         else
     947         {
     948             rl[2] = rl[2] + rl[1][2];
     949             rl[3][size(rl[3])+1] = rl[1][3][1];
     950             rl[1] = rl[1][1];
     951         }
     952     }
     953     S = ring(rl);
     954     setring S;
     955     poly f=imap(Rbs,f);
     956     def I=imap(Rbs,I);
     957     I = simplify(I,2); // eraze the zero generatos
     958     if(f==0)
     959     {
     960         ERROR("minpoly must be non-zero");
     961     }
     962     I=I,f;
     963     def J_I=modStdparallelized(I);
     964     setring Rbs;
     965     def J=imap(S,J_I);
     966     J=simplify(J,2);
     967     return(J);
    508968}
    509969example
    510970{ "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);
     971    ring r1 = (0,a),(x,y),dp;
     972    minpoly = a^2+1;
     973    ideal k = (a/2+1)*x^2+2/3y, 3*x-a*y+ a/7+2;
     974    ideal I = nfmodStd(k);
    515975    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);
     976    ring rm = (0,a),(x,y),(c,dp);
     977    minpoly = a^3+2a+7;
     978    module M = [(a/2+1)*x^2+2/3y, 3*x-a*y+ a/7+2], [ax, y];
     979    M = nfmodStd(M);
     980    M;
     981    ring r2 = (0,a),(x,y,z),dp;
     982    minpoly = a^3 +2;
     983    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;
     984    ideal IJ = nfmodStd(k);
    520985    IJ;
    521     ring r3=0,(x,y),dp;// ring without parameter
     986    ring r3 = 0, (x,y), dp; // ring without parameter
    522987    ideal I = x2 + y, xy - 7y + 2x;
     988    ideal J1 = nfmodStd(I);
     989    J1;
     990    module J2 = nfmodStd(module(I));
     991    J2;
     992    ring r4 = 0, (x,y), (c,dp);
     993    module I = [x2, x-y], [xy,0], [0,-7y + 2x];
    523994    I=nfmodStd(I);
    524995    I;
    525996}
    526 
    527 //////////////////////////////////////////////////////////////////////////////
    528 
    529 /*
    530 Benchmark Problems from
    531 
    532 Boku, Decker, Fieker, Steenpass: Groebner Bases over Algebraic Number
    533 Fields.
    534 
    535 // 1
    536 ring R = (0,a), (x,y,z), dp;
    537 minpoly = (a^2+1);
    538 poly f1 = (a+8)*x^2*y^2+5*x*y^3+(-a+3)*x^3*z
    539           +x^2*y*z;
    540 poly f2 = x^5+2*y^3*z^2+13*y^2*z^3+5*y*z^4;
    541 poly f3 = 8*x^3+(a+12)*y^3+x*z^2+3;
    542 poly f4 = (-a+7)*x^2*y^4+y^3*z^3+18*y^3*z^2;
    543 ideal I1 = f1,f2,f3,f4;
    544 
    545 // 2
    546 ring R = (0,a), (x,y,z), dp;
    547 minpoly = (a^5+a^2+2);
    548 poly f1 = 2*x*y^4*z^2+(a-1)*x^2*y^3*z
    549           +(2*a)*x*y*z^2+7*y^3+(7*a+1);
    550 poly f2 = 2*x^2*y^4*z+(a)*x^2*y*z^2-x*y^2*z^2
    551           +(2*a+3)*x^2*y*z-12*x+(12*a)*y;
    552 poly f3 = (2*a)*y^5*z+x^2*y^2*z-x*y^3*z
    553           +(-a)*x*y^3+y^4+2*y^2*z;
    554 poly f4 = (3*a)*x*y^4*z^3+(a+1)*x^2*y^2*z
    555           -x*y^3*z+4*y^3*z^2+(3*a)*x*y*z^3
    556           +4*z^2-x+(a)*y;
    557 ideal I2 = f1,f2,f3,f4;
    558 
    559 // 3a
    560 ring R = (0,a), (v,w,x,y,z), dp;
    561 minpoly = (a^7-7*a+3);
    562 poly f1 = (a)*v+(a-1)*w+x+(a+2)*y+z;
    563 poly f2 = v*w+(a-1)*w*x+(a+2)*v*y+x*y+(a)*y*z;
    564 poly f3 = (a)*v*w*x+(a+5)*w*x*y+(a)*v*w*z
    565           +(a+2)*v*y*z+(a)*x*y*z;
    566 poly f4 = (a-11)*v*w*x*y+(a+5)*v*w*x*z
    567           +(a)*v*w*y*z+(a)*v*x*y*z
    568           +(a)*w*x*y*z;
    569 poly f5 = (a+3)*v*w*x*y*z+(a+23);
    570 ideal I3a = f1,f2,f3,f4,f5;
    571 
    572 // 3b
    573 ring R = (0,a), (u,v,w,x,y,z), dp;
    574 minpoly = (a^7-7*a+3);
    575 poly f1 = (a)*u+(a+2)*v+w+x+y+z;
    576 poly f2 = u*v+v*w+w*x+x*y+(a+3)*u*z+y*z;
    577 poly f3 = u*v*w+v*w*x+(a+1)*w*x*y+u*v*z+u*y*z
    578           +x*y*z;
    579 poly f4 = (a-1)*u*v*w*x+v*w*x*y+u*v*w*z
    580           +u*v*y*z+u*x*y*z+w*x*y*z;
    581 poly f5 = u*v*w*x*y+(a+1)*u*v*w*x*z+u*v*w*y*z
    582           +u*v*x*y*z+u*w*x*y*z+v*w*x*y*z;
    583 poly f6 = u*v*w*x*y*z+(-a+2);
    584 ideal I3b = f1,f2,f3,f4,f5,f6;
    585 
    586 // 4
    587 ring R = (0,a), (w,x,y,z), dp;
    588 minpoly = (a^6+a^5+a^4+a^3+a^2+a+1);
    589 poly f1 = (a+5)*w^3*x^2*y+(a-3)*w^2*x^3*y
    590           +(a+7)*w*x^2*y^2;
    591 poly f2 = (a)*w^5+(a+3)*w*x^2*y^2
    592           +(a^2+11)*x^2*y^2*z;
    593 poly f3 = (a+7)*w^3+12*x^3+4*w*x*y+(a)*z^3;
    594 poly f4 = 3*w^3+(a-4)*x^3+x*y^2;
    595 ideal I4 = f1,f2,f3,f4;
    596 
    597 // 5
    598 ring R = (0,a), (w,x,y,z), dp;
    599 minpoly = (a^12-5*a^11+24*a^10-115*a^9+551*a^8
    600           -2640*a^7+12649*a^6-2640*a^5+551*a^4
    601           -115*a^3+24*a^2-5*a+1);
    602 poly f1 = (2*a+3)*w*x^4*y^2+(a+1)*w^2*x^3*y*z
    603           +2*w*x*y^2*z^3+(7*a-1)*x^3*z^4;
    604 poly f2 = 2*w^2*x^4*y+w^2*x*y^2*z^2
    605           +(-a)*w*x^2*y^2*z^2
    606           +(a+11)*w^2*x*y*z^3-12*w*z^6
    607           +12*x*z^6;
    608 poly f3 = 2*x^5*y+w^2*x^2*y*z-w*x^3*y*z
    609           -w*x^3*z^2+(a)*x^4*z^2+2*x^2*y*z^3;
    610 poly f4 = 3*w*x^4*y^3+w^2*x^2*y*z^3
    611           -w*x^3*y*z^3+(a+4)*x^3*y^2*z^3
    612           +3*w*x*y^3*z^3+(4*a)*y^2*z^6-w*z^7
    613           +x*z^7;
    614 ideal I5 = f1,f2,f3,f4;
    615 
    616 // 6
    617 ring R = (0,a), (u,v,w,x,y,z), dp;
    618 minpoly = (a^2+5*a+1);
    619 poly f1 = u+v+w+x+y+z+(a);
    620 poly f2 = u*v+v*w+w*x+x*y+y*z+(a)*u+(a)*z;
    621 poly f3 = u*v*w+v*w*x+w*x*y+x*y*z+(a)*u*v
    622           +(a)*u*z+(a)*y*z;
    623 poly f4 = u*v*w*x+v*w*x*y+w*x*y*z+(a)*u*v*w
    624           +(a)*u*v*z+(a)*u*y*z+(a)*x*y*z;
    625 poly f5 = u*v*w*x*y+v*w*x*y*z+(a)*u*v*w*x
    626           +(a)*u*v*w*z+(a)*u*v*y*z+(a)*u*x*y*z
    627           +(a)*w*x*y*z;
    628 poly f6 = u*v*w*x*y*z+(a)*u*v*w*x*y
    629           +(a)*u*v*w*x*z+(a)*u*v*w*y*z
    630           +(a)*u*v*x*y*z+(a)*u*w*x*y*z
    631           +(a)*v*w*x*y*z;
    632 poly f7 = (a)*u*v*w*x*y*z-1;
    633 ideal I6 = f1,f2,f3,f4,f5,f6,f7;
    634 
    635 // 7
    636 ring R = (0,a), (w,x,y,z), dp;
    637 minpoly = (a^8-16*a^7+19*a^6-a^5-5*a^4+13*a^3
    638           -9*a^2+13*a+17);
    639 poly f1 = (-a^2-1)*x^2*y+2*w*x*z-2*w
    640           +(a^2+1)*y;
    641 poly f2 = (a^3-a-3)*w^3*y+4*w*x^2*y+4*w^2*x*z
    642           +2*x^3*z+(a)*w^2-10*x^2+4*w*y-10*x*z
    643           +(2*a^2+a);
    644 poly f3 = (a^2+a+11)*x*y*z+w*z^2-w-2*y;
    645 poly f4 = -w*y^3+4*x*y^2*z+4*w*y*z^2+2*x*z^3
    646           +(2*a^3+a^2)*w*y+4*y^2-10*x*z-10*z^2
    647           +(3*a^2+5);
    648 ideal I7 = f1,f2,f3,f4;
    649 
    650 // 8
    651 ring R = (0,a), (t,u,v,w,x,y,z), dp;
    652 minpoly = (a^7+10*a^5+5*a^3+10*a+1);
    653 poly f1 = v*x+w*y-x*z-w-y;
    654 poly f2 = v*w-u*x+x*y-w*z+v+x+z;
    655 poly f3 = t*w-w^2+x^2-t;
    656 poly f4 = (-a)*v^2-u*y+y^2-v*z-z^2+u;
    657 poly f5 = t*v+v*w+(-a^2-a-5)*x*y-t*z+w*z+v+x+z
    658           +(a+1);
    659 poly f6 = t*u+u*w+(-a-11)*v*x-t*y+w*y-x*z-t-u
    660           +w+y;
    661 poly f7 = w^2*y^3-w*x*y^3+x^2*y^3+w^2*y^2*z
    662           -w*x*y^2*z+x^2*y^2*z+w^2*y*z^2
    663           -w*x*y*z^2+x^2*y*z^2+w^2*z^3-w*x*z^3
    664           +x^2*z^3;
    665 poly f8 = t^2*u^3+t^2*u^2*v+t^2*u*v^2+t^2*v^3
    666           -t*u^3*x-t*u^2*v*x-t*u*v^2*x-t*v^3*x
    667           +u^3*x^2+u^2*v*x^2+u*v^2*x^2
    668           +v^3*x^2;
    669 ideal I8 = f1,f2,f3,f4,f5,f6,f7,f8;
    670 */

  • Singular/singular-libs

    rb1cef7 r204cf0  
    4747        gitfan.lib gradedModules.lib \
    4848        locnormal.lib modnormal.lib modular.lib \
    49         JMBTest.lib JMSConst.lib multigrading.lib\
     49        JMBTest.lib JMSConst.lib multigrading.lib nfmodsyz.lib \
    5050        numerAlg.lib numerDecom.lib \
    5151        orbitparam.lib parallel.lib polymake.lib\
Note: See TracChangeset for help on using the changeset viewer.