////////////////////////////////////////////////////////////////////////////////
version="version nfmodstd.lib 4.1.2.0 Feb_2019 ";  // $Id$
category="Commutative Algebra";
info="

LIBRARY:   nfmodstd.lib  Groebner bases of ideals in polynomial rings
                           over algebraic number fields
AUTHORS:   D.K. Boku       boku@mathematik.uni-kl.de
@*         W. Decker       decker@mathematik.uni-kl.de
@*         C. Fieker       fieker@mathematik.uni-kl.de

OVERVIEW:
  A library for computing the Groebner basis of an ideal in the polynomial
  ring over an algebraic number field Q(t) using the modular methods, where t is
  algebraic over the field of rational numbers Q. For the case Q(t) = Q, the
  procedure is inspired by Arnold [1]. This idea is then extended to the case
  t not in Q using factorization as follows:
  Let f be the minimal polynomial of t.
  For I, I' ideals in Q(t)[X], Q[X,t]/<f> respectively, we map I to I' via the
  map sending  t to t + <f>. We first choose a prime p such that f has at least
  two factors in characteristic p and add each factor f_i to I' to obtain the
  ideal J'_i = I' + <f_i>. We then compute a standard basis G'_i of J'_i for each
  i and combine the G'_i to G_p (a standard basis of I'_p) using chinese remaindering
  for polynomials. The procedure is repeated for many primes p, where we compute the
  G_p in parallel until the number of primes is sufficiently large to recover the
  correct standard basis G' of I'. Finally, by mapping G' back to Q(t)[X], a standard
  basis G of I is obtained.
@*The procedure also works if the input is a module. For this, we consider the
  rings A = Q(t)[X] and A' = (Q[t]/<f>)[X]. For submodules I, I' in A^m, A'^m,
  respectively, we map I to I' via the map sending t to t + <f>. As above, we first
  choose a prime p such that f has at least two factors in characteristic p. For each
  factor f_{i,p} of f_p := (f mod p), we set I'_{i,p} := (I'_p mod f_{i,p}). We then
  compute a standard basis G'_i of I'_{i,p} over F_p[t]/<f_{i,p}> for each i and
  combine the G'_i to G_p (a standard basis of I'_p) using chinese remaindering for
  polynomials. The procedure is repeated for many primes p as described above and we
  finally obtain a standard basis of I.

REFERENCES:
  [1] E. A. Arnold: Modular algorithms for computing Groebner bases.
      J. Symb. Comp. 35, 403-419 (2003).

PROCEDURES:
  chinrempoly(l,m);       chinese remaindering for polynomials
  nfmodStd(I);          standard basis of I over algebraic number field using modular
                        methods
";

LIB "modstd.lib";

////////////////////////////////////////////////////////////////////////////////

static proc testPrime(int p, list L)
{
    /*
     * test whether a prime p divides the denominator(s)
     * and leading coefficients of generating set of ideal
     */
    int i,j,k, tmp;
    def f;
    def I = L[1]; // list L = def I
    number num;
    bigint d1,d2,d3;
    if(typeof(I)=="ideal")
    {
        tmp=1;
    }
    for(i = 1; i <= ncols(I); i++)
    {
        f = cleardenom(I[i]);
        if(f == 0)
        {
             return(0);
        }
        num = leadcoef(I[i])/leadcoef(f);
        d1 = bigint(numerator(num));
        d2 = bigint(denominator(num));
        if( (d1 mod p) == 0)
        {
            return(0);
        }
        if((d2 mod p) == 0)
        {
            return(0);
        }
        if(tmp)
        {
            for(j = size(f); j > 0; j--)
            {
                d3 = bigint(leadcoef(f[j]));
                if( (d3 mod p) == 0)
                {
                    return(0);
                }
            }
        }
        else
        {
            for(j = nrows(f); j > 0; j--)
            {
                for(k=1;k<=size(f[j]);k++)
                {
                    d3 = bigint(leadcoef(f[j][k]));
                    if((d3 mod p) == 0)
                    {
                        return(0);
                    }
                }
            }
        }

    }
    return(1);
}

////////////////////////////////////////////////////////////////////////////////
/* return 1 if the number of factors are in the required bound , 0 else */

static proc minpolyTask(poly f,int p)
{
    /*
     * bound for irreducible factor(s) of (f mod p)
     * see testfact()
     */
    int nr,k,ur;
    ur=deg(f);
    list L=factmodp(f,p);
    if(degtest(L[2])==1)
    {
        // now each factor is squarefree
        if(ur<=3)
        {
            return(1);
        }
        else
        {
            nr = testfact(ur);
            k=ncols(L[1]);
            if(nr < k && k < (ur-nr))// set a bound for k
            {
                return(1);
            }
        }
    }
    return(0);
}

////////////////////////////////////////////////////////////////////////////////
/* return 1 if both testPrime(p,J) and minpolyTask(f,p) is true, 0 else */

static proc PrimeTestTask(int p, list L)
{
    /* L=list(I), I=J,f; J ideal , f minpoly */
    int sz,nr;
    def J=L[1];
    sz=ncols(J);
    def f=J[sz];
    poly g;
    if(typeof(f)=="vector")
    {
        g = f[1];
    }
    else
    {
        g = f;
    }
    if(!testPrime(p,list(J)) or !minpolyTask(g,p))
    {
        return(0);
    }
    return(1);
}

////////////////////////////////////////////////////////////////////////////////
/* compute factors of f mod p with multiplicity */

static proc factmodp(poly f, int p)
{
    def R=basering;
    list l=ringlist(R);
    l[1]=p;
    def S=ring(l);
    setring S;
    list L=factorize(imap(R,f),2);
    ideal J=L[1];
    intvec v=L[2];
    list scx=J,v;
    setring R;
    return(imap(S,scx));
    kill S;
}

////////////////////////////////////////////////////////////////////////////////
/* set a bound for number of factors w.r.t degree nr*/

static proc testfact(int nr)
{
    // nr must be greater than 3
    int i;
    if(nr>3 and nr<=5)
    {
        i=1;
    }
    if(nr>5 and nr<=10)
    {
        i=2;
    }
    if(nr>10 and nr<=15)
    {
        i=3;
    }
    if(nr>15 and nr<=20)
    {
        i=4;
    }
    if(nr>20 and nr<=25)
    {
        i=5;
    }
    if(nr>25 and nr<=30)
    {
        i=6;
    }
    if(nr>30)
    {
        i=10;
    }
    return(i);
}

///////////////////////////////////////////////////////////////////////////////
// return 1 if v[i]>1 , 0 else

static proc degtest(intvec v)
{
    for(int j=1;j<=nrows(v);j++)
    {
        if(v[j]>1)
        {
            return(0);
        }
    }
    return(1);
}

////////////////////////////////////////////////////////////////////////////////

static proc chinRm(list m, list ll, list lk,list l1,int uz)
{
    if(typeof(l1[1])=="ideal" or typeof(l1[1])=="poly")
    {
        poly ff,c;
        for(int i=1;i<=uz;i++)
        {
            c = division(l1[i]*ll[i],m[i])[2][1];
            ff = ff + c*lk[i];
        }
        return(ff);
    }
    else
    {
        vector ff,c;
        for(int i=1;i<=uz;i++)
        {
            c = vector(m[i]);
            attrib(c,"isSB",1);
            ff = ff + (reduce(l1[i]*ll[i],c))*lk[i];
        }
        return(ff);
    }
}

////////////////////////////////////////////////////////////////////////////////

proc chinrempoly(list l,list m)
"USAGE:  chinrempoly(l, m); l list, m list
RETURN:  a polynomial (resp. ideal/module) which is congruent to l[i] modulo m[i]
         for all i
NOTE: The procedure applies chinese remaindering to the first argument w.r.t. the
      moduli given in the second. The elements in the first list must be of the same
      type which can be polynomial, ideal, or module. The moduli must be of type
      polynomial. The elements in the second list must be distinct and co-prime.
SEE ALSO: chinrem
EXAMPLE: example chinrempoly; shows an example
"
{
    int i,j,sz,uz, tmp;
    uz = size(l);
    if(typeof(l[1])=="ideal" or typeof(l[1])=="poly")
    {
        sz = ncols(ideal(l[1]));
        tmp = 1;
    }
    else
    {
        sz = ncols(module(l[1]));
    }
    poly f=1;
    for(i=1;i<=uz;i++)
    {
        f=f*m[i];
    }

    list l1,ll,lk,l2;
    poly c,ff;
    for(j=1;j<=uz;j++)
    {
        lk[j]=f/m[j];
        ll[j]=extgcd(lk[j],m[j])[2];
    }

    if(tmp)
    {
        ideal I,J;
        for(i=1;i<=sz;i++)
        {
            for(j=1;j<=uz;j++)
            {
                I = l[j];
                l1[j] = I[i];
            }
            J[i] = chinRm(m,ll,lk,l1,uz);
        }
        return(J);
    }
    else
    {
        module I,J;
        for(i=1;i<=sz;i++)
        {
            for(j=1;j<=uz;j++)
            {
                I = l[j];
                l1[j] = I[i];
            }
            J[i] = chinRm(m,ll,lk,l1,uz);
        }
        return(J);
    }
}
example
{ "EXAMPLE:"; echo = 2;
    ring rr=97,x,dp;
    poly f=x^7-7*x + 3;
    ideal J=factorize(f,1);
    J;
    list m=J[1..ncols(J)];
    list l= x^2+2*x+3, x^2+5, x^2+7;
    ideal I=chinrempoly(l,m);
    I;
    ring s=0,x,dp;
    list m= x^2+2*x+3, x^3+5, x^4+x^3+7;
    list l=x^3 + 2, x^4 + 7, x^5 + 11;
    ideal I=chinrempoly(l,m);
    I;
    int p=prime(536546513);
    ring r = p, (x,y,a), (dp(2),dp(1));
    poly minpolynomial = a^2+1;
    ideal kf=factorize(minpolynomial,1); //return factors without multiplicity
    kf;
    ideal k=(a+1)*x2+y, 3x-ay+ a+2;
    option(redSB);
    ideal k1=k,kf[1];
    ideal k2 =k,kf[2];
    k1=std(k1);
    k2=std(k2);
    list l=k1,k2;
    list m=kf[1..ncols(kf)];
    ideal I=chinrempoly(l,m);
    I=simplify(I,2);
    I;
    l = module(k1[2..ncols(k1)]), module(k2[2..ncols(k2)]);
    module M = chinrempoly(l,m);
    M;
}
////////////////////////////////////////////////////////////////////////////////

static proc check_leadmonom_and_size(list L)
{
    /*
     * compare the size of ideals in the list and
     * check the corresponding leading monomials
     * size(L)>=2
     */
    def J=L[1];
    int i=size(L);
    int sc=ncols(J);
    int j,k;
    def g=leadmonom(J[1]);
    for(j=1;j<=i;j++)
    {
        if(ncols(L[j])!=sc)
        {
            return(0);
        }
    }
    for(k=2;k<=i;k++)
    {
        for(j=1;j<=sc;j++)
        {
            if(leadmonom(J[j])!=leadmonom(L[k][j]))
            {
                return(0);
            }
        }
    }
    return(1);
}

////////////////////////////////////////////////////////////////////////////////

static proc LiftPolyCRT(def I)
{
    /*
     * compute std for each factor and combine this result
     * to modulo minpoly via CRT for poly over char p>0
     */
    def sl;
    int u,in,j;
    list LL,Lk,T2;
    if(typeof(I)=="ideal")
    {

        ideal J,K,II;
        poly f;
        u=ncols(I);
        J=I[1..u-1];
        f=I[u];
        K=factorize(f,1);
        in=ncols(K);
        for(j=1;j<=in;j++)
        {
            LL[j]=K[j];
            ideal I(j)=J,K[j];
            I(j)=std(I(j));
            if(size(I(j))==1)
            {
                Lk[j]=I(j);
            }
            else
            {
                I(j)[1]=0;
                I(j)=simplify(I(j), 2);
                Lk[j]=I(j);
            }
        }
        if(check_leadmonom_and_size(Lk))
        {
            // apply CRT for polynomials
            II =chinrempoly(Lk,LL),f;
        }
        else
        {
            II=0;
        }
        return(II);
    }
    else
    {
        module J,II;
        vector f;
        u=ncols(I);
        J=I[1..u-1];
        f=I[u];
        poly ff = f[1];
        ideal K=factorize(ff,1);
        in=ncols(K);
        def Ls = basering;
        list l = ringlist(Ls);
        if(l[3][1][1]=="c")
        {
             l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
             list(list(l[3][size(l[3])]))+list(ideal(0));
             l[2] = delete(l[2],size(l[2]));
             l[3] = delete(l[3],size(l[3]));
        }
        else
        {
            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
            list(list(l[3][size(l[3])-1]))+list(ideal(0));
            l[2] = delete(l[2],size(l[2]));
            l[3] = delete(l[3],size(l[3])-1);
        }

        def S1 = ring(l);
        setring S1;
        number Num= number(imap(Ls,ff));
        list l = ringlist(S1);
        l[1][4][1] = Num;
        S1 = ring(l);
        setring S1;
        ideal K = imap(Ls,K);
        def S2;
        module II;
        number Num;
        /* ++++++ if minpoly is irreducible then K will be the zero ideal +++ */
        if(size(K)==0)
        {
            module M = std(imap(Ls,J));
            if(size(M)==1 && M[1]==[1])
            {
                setring Ls;
                return(module([1]));
            }
            II = normalize(M);
        }
        else
        {
            for(j=1;j<=in;j++)
            {
                LL[j]=K[j];
                Num = number(K[j]);
                T2 = ringlist(S1);
                T2[1][4][1] = Num;
                S2 = ring(T2);
                setring S2;
                module M = std(imap(Ls,J));
                if(size(M)== 1 && M[1]==[1])
                {
                    setring Ls;
                    return(module([1]));
                }
                setring S1;
                Lk[j]= imap(S2,M);
            }

            if(check_leadmonom_and_size(Lk))
            {
                // apply CRT for polynomials
                setring Ls;
                II =chinrempoly(imap(S1,Lk),imap(S1,LL));
                setring S1;
                II = normalize(imap(Ls,II));
            }
            else
            {
                setring S1;
                II=[0];
            }
        }
        setring Ls;
        return(imap(S1,II));
    }
}

////////////////////////////////////////////////////////////////////////////////

/* test if 'result' is a GB of the input ideal */
static proc final_Test_minpolyzero(string command, alias list args, module result)
{
    int i;
    list arg = args;
    attrib(result, "isSB", 1);
    for (i = ncols(args[1]); i > 0; i--)
    {
         if (reduce(args[1][i], result, 5) != 0)
         {
              return(0);
         }
    }
    /* test if result is a GB */
    module G = std(result);
    if (size(reduce(G, result,5))!=0) //Modstd::reduce_parallel(G, result)
    {
         return(0);
    }
    return(1);
}

////////////////////////////////////////////////////////////////////////////////

/* test if 'result' is a GB of the input ideal */
static proc final_Test(string command, alias list args, def result)
{
    int i;
    list arg = args;
    // modified for module case
    if(typeof(args[1])=="ideal")
    {
        /* test if args[1] is in result */
        attrib(result, "isSB", 1);
        for (i = ncols(args[1]); i > 0; i--)
        {
            if (reduce(args[1][i], result, 5) != 0)
            {
                return(0);
            }
        }

        /* test if result is a GB */
        ideal G = std(result);
        if (size(reduce(G, result,5))!=0) //Modstd::reduce_parallel(G, result)
        {
            return(0);
        }
        return(1);
    }
    else
    {
        /* test if args[1] is in result */
        def Ls = basering;
        list l = ringlist(Ls);
        if(l[3][1][1]=="c")
        {
             l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
             list(list(l[3][size(l[3])]))+list(ideal(0));
             l[2] = delete(l[2],size(l[2]));
             l[3] = delete(l[3],size(l[3]));
        }
        else
        {
            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
            list(list(l[3][size(l[3])-1]))+list(ideal(0));
            l[2] = delete(l[2],size(l[2]));
            l[3] = delete(l[3],size(l[3])-1);
        }
        def sL = ring(l);
        kill l;
        setring sL;
        list arg = imap(Ls,arg);
        module arg_s =arg[1];
        list l = ringlist(sL);
        l[1][4][1] = arg_s[ncols(arg_s)][1];
        arg_s = arg_s[1..ncols(arg_s)-1];
        def cL = ring(l);
        setring cL;
        module ar_gs = imap(sL,arg_s);
        def Result = imap(Ls,result);
        attrib(Result, "isSB", 1);
        for (i = ncols(ar_gs); i > 0; i--)
        {
            if (reduce(ar_gs[i], Result, 5) != 0)
            {
                setring Ls;
                return(0);
            }
        }
        // test if result is a GB
        module G = std(Result);
        if (size(reduce(G,Result,5))!=0) //Modstd::reduce_parallel(G, Result)
        {
            setring Ls;
            return(0);
        }
        setring Ls;
        return(1);
    }
}

////////////////////////////////////////////////////////////////////////////////

static proc PtestStd_minpolyzero(string command, list args, module result, int p)
{
    /*
     * let G be std of I which is not yet known whether it is the correct
     *  standard basis. So this procedure does the first test
     */
    def br = basering;
    list lbr = ringlist(br);
    if (typeof(lbr[1]) == "int")
    {
        lbr[1] = p;
    }
    else
    {
        lbr[1][1] = p;
    }
    def rp = ring(lbr);
    setring(rp);
    module Ip = imap(br, args)[1];
    int i;
    module Gp = imap(br, result);
    attrib(Gp, "isSB", 1);
    for (i = ncols(Ip); i > 0; i--)
    {
        if (reduce(Ip[i], Gp, 5) != 0)
        {
            setring(br);
            return(0);
        }
    }
    Ip = std(Ip);
    attrib(Ip,"isSB",1);
    for (i = ncols(Gp); i > 0; i--)
    {
        if (reduce(Gp[i], Ip, 5) != 0)
        {
            setring(br);
            return(0);
        }
    }
    setring(br);
    return(1);
}

////////////////////////////////////////////////////////////////////////////////

static proc PtestStd(string command, list args, def result, int p)
{
    /*
     * let G be std of I which is not yet known whether it is the correct
     *  standard basis. So this procedure does the first test
     */
    def br = basering;

    list lbr = ringlist(br);
    if (typeof(lbr[1]) == "int")
    {
        lbr[1] = p;
    }
    else
    {
        lbr[1][1] = p;
    }
    def rp = ring(lbr);
    setring(rp);
    def Ip = imap(br, args)[1];

    int u,in,j,i;
    list LL,Lk,T2;
    if(typeof(Ip)!="ideal")
    {
        module J,II;
        vector f;
        u=ncols(Ip);
        J=Ip[1..u-1];
        f=Ip[u];
        poly ff = f[1];
        ideal K=factorize(ff,1);
        in=ncols(K);
        def Ls = basering;
        list l = ringlist(Ls);
        if(l[3][1][1]=="c")
        {
            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
            list(list(l[3][size(l[3])]))+list(ideal(0));
            l[2] = delete(l[2],size(l[2]));
            l[3] = delete(l[3],size(l[3]));
        }
        else
        {
            l[1] = list(l[1]) + list(list(l[2][size(l[2])])) +
            list(list(l[3][size(l[3])-1]))+list(ideal(0));
            l[2] = delete(l[2],size(l[2]));
            l[3] = delete(l[3],size(l[3])-1);
        }
        def S1 = ring(l);
        setring S1;
        number Num= number(imap(Ls,ff));
        list l = ringlist(S1);
        l[1][4][1] = Num;
        S1 = ring(l);
        setring S1;
        ideal K = imap(Ls,K);
        module Jp = imap(Ls,J);
        def S2;
        module Ip;
        number Num;
        /* ++++++ if the minpoly is irreducible then K = ideal(0) +++ */
        if(size(K)==0)
        {
            module M = std(Jp);
            Ip = normalize(M);
        }
        else
        {
            for(j=1;j<=ncols(K);j++)
            {
                LL[j]=K[j];
                Num = number(K[j]);
                T2 = ringlist(S1);
                T2[1][4][1] = Num;
                S2 = ring(T2);
                setring S2;
                module M = std(imap(Ls,J));
                setring S1;
                Lk[j]= imap(S2,M);
            }
            if(check_leadmonom_and_size(Lk))
            {
                // apply CRT for polynomials
                setring Ls;
                II =chinrempoly(imap(S1,Lk),imap(S1,LL));
                setring S1;
                Ip = normalize(imap(Ls,II));
            }
            else
            {
                setring S1;
                Ip=[0];
            }
        }
        setring S1;
        module Gp = imap(br, result);
        attrib(Gp, "isSB", 1);
        for (i = ncols(Jp); i > 0; i--)
        {
            if (reduce(Jp[i], Gp, 5) != 0)
            {
                setring(br);
                return(0);
            }
        }

        attrib(Ip,"isSB",1);
        for (i = ncols(Gp); i > 0; i--)
        {
            if (reduce(Gp[i], Ip, 5) != 0)
            {
                setring(br);
                return(0);
            }
        }
        setring(br);
        return(1);
    }
    else
    {
        ideal Gp = imap(br, result);
        attrib(Gp, "isSB", 1);
        for (i = ncols(Ip); i > 0; i--)
        {
            if (reduce(Ip[i], Gp, 5) != 0)
            {
                setring(br);
                return(0);
            }
        }
        Ip = LiftPolyCRT(Ip);
        attrib(Ip,"isSB",1);
        for (i = ncols(Gp); i > 0; i--)
        {
            if (reduce(Gp[i], Ip, 5) != 0)
            {
                setring(br);
                return(0);
            }
        }
        setring(br);
        return(1);
    }
}

////////////////////////////////////////////////////////////////////////////////

static proc cleardenomIdeal(def I)
{
    int t=ncols(I);
    if(size(I)==0)
    {
        return(I);
    }
    else
    {
        for(int i=1;i<=t;i++)
        {
            I[i]=cleardenom(I[i]);
        }
    }
    return(I);
}

////////////////////////////////////////////////////////////////////////////////

static proc modStdparallelized(def I, list #)
{
    // apply modular.lib
    /* save options */
    intvec opt = option(get);
    option(redSB);
    if(size(#)>0)
    {
        I = modular("std", list(I), testPrime, Modstd::deleteUnluckyPrimes_std,
                PtestStd_minpolyzero, final_Test_minpolyzero,536870909);
    }
    else
    {
        I = modular("Nfmodstd::LiftPolyCRT", list(I), PrimeTestTask,
                    Modstd::deleteUnluckyPrimes_std,PtestStd, final_Test,536870909);
    }
    attrib(I, "isSB", 1);
    option(set,opt);
    return(I);
}

////////////////////////////////////////////////////////////////////////////////
/* main procedure */
proc nfmodStd(def I, list #)
"USAGE:  nfmodStd(I, #); I ideal or module, # optional parameters
RETURN:  standard basis of I over algebraic number field
NOTE: The procedure passes to @ref{modStd} if the ground field has no
      parameter. In this case, the optional parameters # (if given)
      are directly passed to @ref{modStd}.
SEE ALSO: modStd
EXAMPLE: example nfmodStd; shows an example
"
{
     list L=#;
     def Rbs=basering;
     poly f;
     int tmp=1;
     if(typeof(I)!="ideal")
     {
         tmp =0;
     }
     int n=nvars(Rbs);
     if(size(I)==0)
     {
         if(!tmp)
         {
             return(module([0]));
         }
         return(ideal(0));
     }
     if(npars(Rbs)==0)
     {
         //if algebraic number is in Q
         if(typeof(I)=="module")
         {
             def J = modStdparallelized(I,1);
             return(J);
         }
         else
         {
             def J=modStd(I,L);
             return(J);
         }
     }
     def S;
     list rl=ringlist(Rbs);
     f=rl[1][4][1];
     if(tmp)
     {
         rl[2][n+1]=rl[1][2][1];
         rl[1]=rl[1][1];
         rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
         rl[3][size(rl[3])-1]=list("dp",1);
     }
     else
     {
         if(rl[3][1][1]!="c")
         {
             rl[2] = rl[2] + rl[1][2];
             rl[3] = insert(rl[3], rl[1][3][1],1);
             rl[1] = rl[1][1];
         }
         else
         {
             rl[2] = rl[2] + rl[1][2];
             rl[3][size(rl[3])+1] = rl[1][3][1];
             rl[1] = rl[1][1];
         }
     }
     S = ring(rl);
     setring S;
     poly f=imap(Rbs,f);
     def I=imap(Rbs,I);
     I = simplify(I,2); // eraze the zero generatos
     if(f==0)
     {
         ERROR("minpoly must be non-zero");
     }
     I=I,f;
     def J_I=modStdparallelized(I);
     setring Rbs;
     def J=imap(S,J_I);
     J=simplify(J,2);
     return(J);
}
example
{ "EXAMPLE:"; echo = 2;
    ring r1 = (0,a),(x,y),dp;
    minpoly = a^2+1;
    ideal k = (a/2+1)*x^2+2/3y, 3*x-a*y+ a/7+2;
    ideal I = nfmodStd(k);
    I;
    ring rm = (0,a),(x,y),(c,dp);
    minpoly = a^3+2a+7;
    module M = [(a/2+1)*x^2+2/3y, 3*x-a*y+ a/7+2], [ax, y];
    M = nfmodStd(M);
    M;
    ring r2 = (0,a),(x,y,z),dp;
    minpoly = a^3 +2;
    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;
    ideal IJ = nfmodStd(k);
    IJ;
    ring r3 = 0, (x,y), dp; // ring without parameter
    ideal I = x2 + y, xy - 7y + 2x;
    ideal J1 = nfmodStd(I);
    J1;
    module J2 = nfmodStd(module(I));
    J2;
    ring r4 = 0, (x,y), (c,dp);
    module I = [x2, x-y], [xy,0], [0,-7y + 2x];
    I=nfmodStd(I);
    I;
}