////////////////////////////////////////////////////////////////////////////////
version="version nfmodstd.lib 4.0.1.0 Sep_2014 ";  // $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.

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, ideal I)
{
    /*
     * test whether a prime p divides the denominator(s)
     * and leading coefficients of generating set of ideal
     */
    int i,j;
    poly f;
    number num;
    bigint d1,d2,d3;
    for(i = 1; i <= size(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);
        }
        for(j = size(f); j > 0; j--)
        {
            d3 = bigint(leadcoef(f[j]));
            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 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,dg;
    ideal J=L[1];
    sz=ncols(J);
    poly f=J[sz];
    dg=deg(f);
    if(!testPrime(p,J) or !minpolyTask(f,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)
{
    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);
}

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

proc chinrempoly(list l,list m)
"USAGE:  chinrempoly(l, m); l list, m list
RETURN:  a polynomial (resp. ideal) 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 same type
      which can be polynomial or ideal. The moduli must be of type polynomial. 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;
    uz = size(l);
    sz = ncols(ideal(l[1]));
    poly f=1;
    for(i=1;i<=uz;i++)
    {
        f=f*m[i];
    }

    ideal I,J;
    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];
    }

    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;
}

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

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
     */
    ideal J=L[1];
    int i=size(L);
    int sc=ncols(J);
    int j,k;
    poly 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(ideal I)
{
    /*
     * compute std for each factor and combine this result
     * to modulo minpoly via CRT for poly over char p>0
     */
    int u,in,j;
    list LL,Lk;
    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);
}

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

static proc PtestStd(string command, list args, ideal result, int p)
{
    /*
     * let G be std of I which is not yet known whether it is the correct
     *  standard basis or not. 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);
    ideal Ip = fetch(br, args)[1];
    ideal Gp = fetch(br, result);
    attrib(Gp, "isSB", 1);
    int i;
    for (i = ncols(Ip); i > 0; i--)
    {
        if (reduce(Ip[i], Gp, 1) != 0)
        {
            setring(br);
            return(0);
        }
    }
    Ip = LiftPolyCRT(Ip);
    attrib(Ip,"isSB",1);
    for (i = ncols(Gp); i > 0; i--)
    {
        if (reduce(Gp[i], Ip, 1) != 0)
        {
            setring(br);
            return(0);
        }
    }
    setring(br);
    return(1);
}

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

static proc cleardenomIdeal(ideal 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(ideal I)
{
    // apply modular.lib
    /* save options */
    intvec opt = option(get);
    option(redSB);
    I = modular("Nfmodstd::LiftPolyCRT", list(I), PrimeTestTask, Modstd::deleteUnluckyPrimes_std,
              PtestStd, Modstd::finalTest_std,536870909);
    attrib(I, "isSB", 1);
    option(set,opt);
    return(I);
}

////////////////////////////////////////////////////////////////////////////////
/* main procedure */
proc nfmodStd(ideal I, list #)
"USAGE:  nfmodStd(I, #); I ideal, # 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;
    ideal J;
    int n=nvars(Rbs);
    if(size(I)==0)
    {
        return(ideal(0));
    }
    if(npars(Rbs)==0)
    {
        J=modStd(I,L);//if algebraic number is in Q
        if(size(#)>0)
        {
            return(cleardenomIdeal(J));
        }
        return(J);
    }
    list rl=ringlist(Rbs);
    f=rl[1][4][1];
    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);
    def S=ring(rl);
    setring S;
    poly f=imap(Rbs,f);
    ideal I=imap(Rbs,I);
    I = simplify(I,2);// eraze the zero generatos
    ideal J;
    if(f==0)
    {
        ERROR("minpoly must be non-zero");
    }
    I=I,f;
    J=modStdparallelized(I);
    setring Rbs;
    J=imap(S,J);
    J=simplify(J,2);
    if(size(#)>0)
    {
        return(cleardenomIdeal(J));
    }
    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 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;
    I=nfmodStd(I);
    I;
}

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

/*
Benchmark Problems from

Boku, Decker, Fieker, Steenpass: Groebner Bases over Algebraic Number
Fields.

// 1
ring R = (0,a), (x,y,z), dp;
minpoly = (a^2+1);
poly f1 = (a+8)*x^2*y^2+5*x*y^3+(-a+3)*x^3*z
          +x^2*y*z;
poly f2 = x^5+2*y^3*z^2+13*y^2*z^3+5*y*z^4;
poly f3 = 8*x^3+(a+12)*y^3+x*z^2+3;
poly f4 = (-a+7)*x^2*y^4+y^3*z^3+18*y^3*z^2;
ideal I1 = f1,f2,f3,f4;

// 2
ring R = (0,a), (x,y,z), dp;
minpoly = (a^5+a^2+2);
poly f1 = 2*x*y^4*z^2+(a-1)*x^2*y^3*z
          +(2*a)*x*y*z^2+7*y^3+(7*a+1);
poly f2 = 2*x^2*y^4*z+(a)*x^2*y*z^2-x*y^2*z^2
          +(2*a+3)*x^2*y*z-12*x+(12*a)*y;
poly f3 = (2*a)*y^5*z+x^2*y^2*z-x*y^3*z
          +(-a)*x*y^3+y^4+2*y^2*z;
poly f4 = (3*a)*x*y^4*z^3+(a+1)*x^2*y^2*z
          -x*y^3*z+4*y^3*z^2+(3*a)*x*y*z^3
          +4*z^2-x+(a)*y;
ideal I2 = f1,f2,f3,f4;

// 3a
ring R = (0,a), (v,w,x,y,z), dp;
minpoly = (a^7-7*a+3);
poly f1 = (a)*v+(a-1)*w+x+(a+2)*y+z;
poly f2 = v*w+(a-1)*w*x+(a+2)*v*y+x*y+(a)*y*z;
poly f3 = (a)*v*w*x+(a+5)*w*x*y+(a)*v*w*z
          +(a+2)*v*y*z+(a)*x*y*z;
poly f4 = (a-11)*v*w*x*y+(a+5)*v*w*x*z
          +(a)*v*w*y*z+(a)*v*x*y*z
          +(a)*w*x*y*z;
poly f5 = (a+3)*v*w*x*y*z+(a+23);
ideal I3a = f1,f2,f3,f4,f5;

// 3b
ring R = (0,a), (u,v,w,x,y,z), dp;
minpoly = (a^7-7*a+3);
poly f1 = (a)*u+(a+2)*v+w+x+y+z;
poly f2 = u*v+v*w+w*x+x*y+(a+3)*u*z+y*z;
poly f3 = u*v*w+v*w*x+(a+1)*w*x*y+u*v*z+u*y*z
          +x*y*z;
poly f4 = (a-1)*u*v*w*x+v*w*x*y+u*v*w*z
          +u*v*y*z+u*x*y*z+w*x*y*z;
poly f5 = u*v*w*x*y+(a+1)*u*v*w*x*z+u*v*w*y*z
          +u*v*x*y*z+u*w*x*y*z+v*w*x*y*z;
poly f6 = u*v*w*x*y*z+(-a+2);
ideal I3b = f1,f2,f3,f4,f5,f6;

// 4
ring R = (0,a), (w,x,y,z), dp;
minpoly = (a^6+a^5+a^4+a^3+a^2+a+1);
poly f1 = (a+5)*w^3*x^2*y+(a-3)*w^2*x^3*y
          +(a+7)*w*x^2*y^2;
poly f2 = (a)*w^5+(a+3)*w*x^2*y^2
          +(a^2+11)*x^2*y^2*z;
poly f3 = (a+7)*w^3+12*x^3+4*w*x*y+(a)*z^3;
poly f4 = 3*w^3+(a-4)*x^3+x*y^2;
ideal I4 = f1,f2,f3,f4;

// 5
ring R = (0,a), (w,x,y,z), dp;
minpoly = (a^12-5*a^11+24*a^10-115*a^9+551*a^8
          -2640*a^7+12649*a^6-2640*a^5+551*a^4
          -115*a^3+24*a^2-5*a+1);
poly f1 = (2*a+3)*w*x^4*y^2+(a+1)*w^2*x^3*y*z
          +2*w*x*y^2*z^3+(7*a-1)*x^3*z^4;
poly f2 = 2*w^2*x^4*y+w^2*x*y^2*z^2
          +(-a)*w*x^2*y^2*z^2
          +(a+11)*w^2*x*y*z^3-12*w*z^6
          +12*x*z^6;
poly f3 = 2*x^5*y+w^2*x^2*y*z-w*x^3*y*z
          -w*x^3*z^2+(a)*x^4*z^2+2*x^2*y*z^3;
poly f4 = 3*w*x^4*y^3+w^2*x^2*y*z^3
          -w*x^3*y*z^3+(a+4)*x^3*y^2*z^3
          +3*w*x*y^3*z^3+(4*a)*y^2*z^6-w*z^7
          +x*z^7;
ideal I5 = f1,f2,f3,f4;

// 6
ring R = (0,a), (u,v,w,x,y,z), dp;
minpoly = (a^2+5*a+1);
poly f1 = u+v+w+x+y+z+(a);
poly f2 = u*v+v*w+w*x+x*y+y*z+(a)*u+(a)*z;
poly f3 = u*v*w+v*w*x+w*x*y+x*y*z+(a)*u*v
          +(a)*u*z+(a)*y*z;
poly f4 = u*v*w*x+v*w*x*y+w*x*y*z+(a)*u*v*w
          +(a)*u*v*z+(a)*u*y*z+(a)*x*y*z;
poly f5 = u*v*w*x*y+v*w*x*y*z+(a)*u*v*w*x
          +(a)*u*v*w*z+(a)*u*v*y*z+(a)*u*x*y*z
          +(a)*w*x*y*z;
poly f6 = u*v*w*x*y*z+(a)*u*v*w*x*y
          +(a)*u*v*w*x*z+(a)*u*v*w*y*z
          +(a)*u*v*x*y*z+(a)*u*w*x*y*z
          +(a)*v*w*x*y*z;
poly f7 = (a)*u*v*w*x*y*z-1;
ideal I6 = f1,f2,f3,f4,f5,f6,f7;

// 7
ring R = (0,a), (w,x,y,z), dp;
minpoly = (a^8-16*a^7+19*a^6-a^5-5*a^4+13*a^3
          -9*a^2+13*a+17);
poly f1 = (-a^2-1)*x^2*y+2*w*x*z-2*w
          +(a^2+1)*y;
poly f2 = (a^3-a-3)*w^3*y+4*w*x^2*y+4*w^2*x*z
          +2*x^3*z+(a)*w^2-10*x^2+4*w*y-10*x*z
          +(2*a^2+a);
poly f3 = (a^2+a+11)*x*y*z+w*z^2-w-2*y;
poly f4 = -w*y^3+4*x*y^2*z+4*w*y*z^2+2*x*z^3
          +(2*a^3+a^2)*w*y+4*y^2-10*x*z-10*z^2
          +(3*a^2+5);
ideal I7 = f1,f2,f3,f4;

// 8
ring R = (0,a), (t,u,v,w,x,y,z), dp;
minpoly = (a^7+10*a^5+5*a^3+10*a+1);
poly f1 = v*x+w*y-x*z-w-y;
poly f2 = v*w-u*x+x*y-w*z+v+x+z;
poly f3 = t*w-w^2+x^2-t;
poly f4 = (-a)*v^2-u*y+y^2-v*z-z^2+u;
poly f5 = t*v+v*w+(-a^2-a-5)*x*y-t*z+w*z+v+x+z
          +(a+1);
poly f6 = t*u+u*w+(-a-11)*v*x-t*y+w*y-x*z-t-u
          +w+y;
poly f7 = w^2*y^3-w*x*y^3+x^2*y^3+w^2*y^2*z
          -w*x*y^2*z+x^2*y^2*z+w^2*y*z^2
          -w*x*y*z^2+x^2*y*z^2+w^2*z^3-w*x*z^3
          +x^2*z^3;
poly f8 = t^2*u^3+t^2*u^2*v+t^2*u*v^2+t^2*v^3
          -t*u^3*x-t*u^2*v*x-t*u*v^2*x-t*v^3*x
          +u^3*x^2+u^2*v*x^2+u*v^2*x^2
          +v^3*x^2;
ideal I8 = f1,f2,f3,f4,f5,f6,f7,f8;
*/