singular-data 1:4.1.0-p3+ds-2build1 / usr / share / singular / LIB / ffmodstd.lib

///////////////////////////////////////////////////////////////////////////////
version="version ffmodstd.lib 4.0.3.0 May_2016 "; // $Id: 35bb14a4aa2f7b5d6ad642767a7026bdf84b6777 $
category="Commutative Algebra";
info="
LIBRARY:   ffmodstd.lib  Groebner bases of ideals in polynomial rings
                           over rational function 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 a Groebner basis of an ideal in a polynomial
  ring over an algebraic function field Q(T):=Q(t_1,...,t_m) using modular
  methods and sparse multivariate rational interpolation, where the
  t_i are transcendental over Q. The idea is as follows:
  Given an ideal I in Q(T)[X], we map I to J via the map sending T to
  Tz:=(t_1z+s_1,..., t_mz+s_m) for a suitable point s in Q^m\\{(0,...,0)} and for some
  extra variable z so that J is an ideal in Q(Tz)[X]. For a suitable point b in
  Z^m\\{(0,...,0)}, we map J to K via the map sending (T,z) to (b,z), where
  b:=(b_1,...,b_m) (usually the b_i's are distinct primes), so that K is an ideal in
  Q(z)[X]. For such a rational point b, we compute a Groebner basis G_b of K using
  modular algorithms [1] and univariate rational interpolation [2,7]. The
  procedure is repeated for many rational points b until their number is sufficiently
  large to recover the correct coeffcients in Q(T). Once we have these points, we
  obtain a set of polynomials G by applying the sparse multivariate rational interpolation
  algorithm from [4] coefficient-wise to the list of Groebner bases G_b in Q(z)[X],
  where this algorithm makes use of the following algorithms: univariate polynomial
  interpolation [2], univariate rational function reconstruction [7], and multivariate
  polynomial interpolation [3]. The last algorithm uses the well-known Berlekamp/Massey
  algorithm [5] and its early termination version [6]. The set G is then a Groebner
  basis of I with high probability.

REFERENCES:
    [1] E. A. Arnold: Modular algorithms for computing Groebner bases.
      J. Symb. Comput. 35, 403-419 (2003).
 @* [2] R. L. Burden and J. D. Faires: Numerical analysis. 9th ed. (1993).
 @* [3] M. Ben-Or and P. Tiwari: A deterministic algorithm for sparse multivariate
      polynomial interpolation. Proc. of the 20th Annual ACM Symposium on
      Theory of Computing, 301-309 (1988).
 @* [4] A. Cuyt and W.-s. Lee: Sparse interpolation of multivariate rational functions.
      Theor. Comput. Sci. 412, 1445-1456 (2011).
 @* [5] E. Kaltofen and W.-s. Lee: Early termination in sparse interpolation algorithms.
      J. Symb. Comput. 36, 365-400 (2003).
 @* [6] E. Kaltofen, W.-s. Lee and A. A. Lobo: Early termination in Ben-Or/Tiwari
      sparse interpolation and a hybrid of Zippel's algorithm. Proc. ISSAC
      (ISSAC '00), 192-201 (2000).
 @* [7] K. Sara and M. Monagan: Fast Rational Function Reconstruction. Proc. ISSAC
      (ISSAC '06), 184-190 (2006).

PROCEDURES:
  fareypoly(g,f);         univariate rational function reconstruction
  polyInterpolation(l,m);  univariate polynomial interpolation
  BerlekampMassey(L,i);    Berlekamp/Massey algorithm
  sparseInterpolation(f,L,n); sparse multivariate polynomial interpolation
  ffmodStd(I);          Groebner bases over algebraic function fields
                          using modular methods and sparse multivariate rational
                          interpolation
";
LIB "modstd.lib";
LIB "linalg.lib";


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

static proc collect_coeffs(ideal I)
{
    // return the numerators of the coefficients in I
    list L, J1;
    int i,j;
    poly g;
    number n1,N;
    while(i < ncols(I))
    {
        i++;
        g=I[i];
        for(j=1;j<=size(g);j++)
        {
            N=leadcoef(g[j]);
            n1=numerator(N);
            J1=J1+list(n1);
        }
    }
    return(J1);
}

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

static proc Testlist_all(list L)
{
    // discard all constants from the list L
    // base ring is a polynomial ring over polynomal ring with block ordering
    ideal I=0;
    for(int j=1;j<=size(L);j++)
    {
        if(deg(L[j])>0)
        {
            I=I+L[j];
        }
    }
    return(I);
}

////////////////////////////////////////////////////////////////////////////////
// ++++++++++++++++++++  polynomial Interpolation +++++++++++++++

proc polyInterpolation(list d, list e,list #)
"USAGE:  polyInterpolation(d, e[, n, L]); d list, e list, n int, L list
RETURN:  a list l_p where f:=l_p[1] is a polynomial of degree at most size(d)-1
         which satisfies the conditions f(d[i])=e[i] for all i, l_p[2] is the product
         of all (var(n)-d[i]) for 1 <= i <= size(d) and l_p[3]=d.
NOTE: The procedure applies the Newton interpolation algorithm to the pair (d,e)
      and returns the output w.r.t. the first variable (default) of the ground
      ring. If an optional parameter n, 1<=n<=N (N is the number of variables in the
      current basering), is given, then the procedure returns the list l_p w.r.t. the
      n-th variable. Moreover, if the number of points (d'[i],e'[i]) is not large enough
      to obtain the target polynomial, L = polyInterpolation(d', e', n) can be provided
      as an optional parameter to add more interpolation points.
      The elements in the first list must be distinct.
EXAMPLE: example polyInterpolation; shows an example
"
{
    /*  compute a polynomial from given numerical data
     *  size of d and e must be equal
     *  d is list of distinct elements
     */

    // optional parameters

    int vr,i,dt,j;
    int sz=size(#);
    int s_d=size(d);
    poly f,g;
    list l_p,ltd;
    f=e[1];
    l_p=f,g,d;
    vr =1;
    number s,t;
    if(sz<=1)
    {
        if(sz)
        {
            vr=#[1];
        }
        if(s_d==1)
        {
            return(l_p);
        }
        g = (var(vr)-d[1]);
        for(j=2;j<=s_d;j++)
        {
            s = (d[j]-d[1]);
            t = e[j] - number(subst(f,var(vr),d[j]));
            for(i=2;i< j;i++)
            {
                s= s*(d[j]-d[i]);
            }
            t = t/s;
            f = f + t*g;
            g = g*(var(vr)-d[j]);
        }
        l_p=f,g,d;
        return(l_p);
    }
    else
    {
        // ================ interpolate at additional points ======
        vr = #[1];
        # = #[2];
        ltd = #[3]+d;
        dt = size(#[3]);
        f = #[1];
        g = #[2];
        for(j=1;j<=s_d;j++)
        {
            s = d[j]-ltd[1];
            t = e[j] - number(subst(f,var(vr),d[j]));
            for(i=2;i < dt+j;i++)
            {
                s= s*(ltd[dt+j]-ltd[i]);
            }
            t = t/s;
            f = f + t*g;
            g = g*(var(vr)-ltd[dt+j]);
        }
        l_p=f,g,ltd;
        return(l_p);
    }
}
example
{ "EXAMPLE:"; echo = 2;
    ring rr=23,(x,y),dp;
    list d = 1,2,3,4;
    list e = -1,10,3,8;
    polyInterpolation(d,e);
    polyInterpolation(d,e,2)[1];
    list d1 = 5,6;
    list e1 = -7,6;
    list L = polyInterpolation(d,e);
    L = polyInterpolation(d1,e1,1,L); // add points
    L;
    ring R = (499,a),x,dp;
    list d2 = 2,3a,5;
    list e2 = (a-2), (9a2-8a), (a+10);
    polyInterpolation(d2,e2);
}

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

static proc NewtonInterpolationNormal(list d, list e,int vr)
{
    /* compute a polynomial from given numerical data
     * size of d and e must be equal
     * d is list of distinct elements
     * number of points here are sufficiently large
     */
    int i;
    int s_d=size(d);
    poly f=e[1];
    poly g=(var(vr)-d[1]);
    number s,t;
    if(s_d==1)
    {
        return(f);
    }
    for(int j=2;j<=s_d;j++)
    {
        s = (d[j]-d[1]);
        t = e[j] - number(subst(f,var(vr),d[j]));
        for(i=2;i< j;i++)
        {
            s= s*(d[j]-d[i]);
        }
        t = t/s;
        f = f + t*g;
        g = g*(var(vr)-d[j]);
    }
    return(list(f,g));
}

///////////////////////////////////////////////////////////////////////////////
// +++++++++++++++ choose a shift +++++++++++++++

static proc test_the_shift(ideal I, int n, int pa)
{
    /* generators of I are the coefficients of a given input ideal  w.r.t main
     * variables n
     * pa is the number of parameters */

    list sh = choose_a_shift(pa);
    if(size(I)==0)
    {
        return(sh);
    }
    ideal J=Evaluate_givenI(I,sh,1,n);
    int i;
    while(size(J)!=ncols(I))
    {
        i++;
        sh = choose_a_shift(pa);
        J=Evaluate_givenI(I,sh,1,n);
    }
    return(sh);
}

///////////////////////////////////////////////////////////////////////////////
// evaluate the polynomials in J at given values

static proc Evaluate_givenI(ideal J,list pr, int i,int n)
{
    // n is the number of main variables
    int sz=ncols(J);
    int sr=size(pr);
    int k;
    for(int j=1;j<=sz;j++)
    {
        for(k=n+1;k<=n+sr;k++)
        {
            J=subst(J,var(k),number(pr[k-n])**i);
        }
    }
    return(J);
}

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

static proc choose_a_shift(int pa)
{
    // choose a shift w.r.t pa (the number of parameters)
    list h = random(100,150);
    for(int i=2;i<=pa;i++)
    {
        h[i]=random(h[i-1]+2, h[i-1]+7);
    }
    return(h);
}

///////////////////////////////////////////////////////////////////////////////
// +++++++++++++++ choose a random distinct primes +++++++++++++++

static proc choose_a_prime(int p)
{
    // p must be a prime number and the procedure returns the next prime
    if(p==2){return(3);}
    int i=p;
    while(1)
    {
        i = i+2;
        if(prime(i)==i)
        {
            return(i);
        }
    }
}

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

static proc list_of_primes(int pa, list #)
{
    // find distinct pa prime(s) up to permutations
    int p=3;
    int j,k;
    list L,l,l1;
    if(size(#)>0)
    {
        p = #[1];
    }
    for(j=1;j<=pa;j++)
    {
        L[j]=p;
        p=choose_a_prime(p);
    }
    l1 = L;
    for(j=1;j<=size(L);j++)
    {

        k = random(1,size(l1));
        l = l+list(l1[k]);
        l1 = delete(l1,k);
    }
    return(l);
}

///////////////////////////////////////////////////////////////////////////////
// +++++++++++ BerlekampMassey Algorithm +++++++++++++++++++++++

static proc reverse_coef(poly f, int i)
{
    /* keeping the monomials (in f) fixed returns a polynomial g by reversing
     * the coeffcients in f, example (1,2,3) to (3,2,1) */
    poly g;
    for(int j=size(f);j>=1;j--)
    {
        g = g + var(i)**(j-1)*leadcoef(f[j]);
    }
    return(g);
}

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

static proc rev_coef_new(poly f, int i)
{
    /* return a list of numbers which are the coefficients of f
     * starting from deg(f) */
    matrix M=coeffs(f,var(i));
    int t=nrows(M);
    list  L=number(M[1..t,1]);
    L=L[t..1];
    return(L);
}

///////////////////////////////////////////////////////////////////////////////
// compute the minimal polynomial of L using Berlekamp/Massey algorithm

proc BerlekampMassey(list L, int i,list #)
"USAGE:  BerlekampMassey(L, i[, M]); L list, i int, M list
RETURN:  a list Tr where f:=Tr[1] is the minimal polynomial (w.r.t. the i-th variable)
         generated by the sequence (L[j]), 1<=j<= Tr[2], if the
         length of the sequence is long enough. In this case, the coefficients c_i of
         the polynomial f satisfy the relation -L[j+t] = c_0*L[j] + ... + c_{t-1}*L[j+t-1]
         for all j >=1 where t=deg(f).
NOTE: The procedure applies the Berlekamp/Massey algorithm to the sequence L[j]
      (elements from the field Q) for j>0 and returns a polynomial f. If the polynomial
      f splits into linear factors with no multiplicity greater than one, then we say that
      the length of the sequence L is long enough. If this polynomial does not split into
      linear factors, an optional parameter M = BerlekampMassey(L',i) can be provided to
      add more elements to the sequence.
REFERENCES:
@* [1] E. Kaltofen and W.-s. Lee: Early termination in sparse interpolation
      algorithms. J. Symb. Comput. 36, 365-400 (2003).
@* [2] E. Kaltofen, W.-s. Lee and A. A. Lobo: Early termination in Ben-Or/Tiwari
  sparse interpolation and a hybrid of Zippel's algorithm. Proc. ISSAC
  (ISSAC '00), 192-201 (2000).
EXAMPLE: example BerlekampMassey; shows an example
"
{
    /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     * L is stream (sequence), unbounded,  of elements from any field
     * i is variable position
     * Note that we may not obtain the minimal polynomial because the length of
     *  the sequence may not be long enough so we need to update it this is where
     *  we need the optional parameters  +++++++++++++++ */

    list la,Tr,Z;
    int s,j,k,n, sp;
    number De=1;
    number d;
    if(size(#)==0)
    {
        n=size(L);
        poly g(0..n);
        poly B(0..n);
        int l(0..n);
        number D(1..n);
        B(1)=0;
        l(1)=0;
        g(0)=1;
   }
   else
   {
        /*********************************************************************
         ************  update BerlekampMassey procedure *********************/
        n=size(L);
        list M=#[6];
        sp = #[5];
        int ik=n+sp;
        M[sp+1..ik]=L[1..n];
        L=M;
        poly g(0..ik);
        poly B(0..ik);
        int l(0..ik);
        number D(1..ik);
        g(sp)=#[1];
        B(sp)=#[2];
        l(sp)=#[3];
        De=#[4];
   }
   for(j=sp+1;j<=n+sp;j++)
   {
       s=deg(g(j-1));
       la=rev_coef_new(g(j-1),i);
       d=0;
       for(k=0;k<=s;k++)
       {
           d = d + la[k+1]*number(L[j-s+k]);
       }
       D(j)=d;
       if(D(j)==0)
       {
           if((2*l(j-1)) < j && j>1)
           {
               return(list(reverse_coef (g(j-1),i),j-1));
           }
           g(j)=g(j-1);
           B(j)=var(i)*B(j-1);
           l(j)=l(j-1);
       }
       else
       {
           if(D(j)!=0 && (2*l(j-1)) < j)
           {
               B(j)=g(j-1);
               g(j)=g(j-1)- (D(j)*var(i)*B(j-1))/De;
               l(j)=j-l(j-1);
               De=D(j);
           }
           else
           {
               if( D(j)!=0 && (2*l(j-1)) >= j)
               {
                   g(j)=g(j-1) - (D(j)*var(i)*B(j-1))/De;
                   B(j)=var(i)*B(j-1);
                   l(j)=l(j-1);
               }
           }
       }
    }
    Tr=g(n+sp),B(n+sp),l(n+sp),De,n+sp,L;
    return(Tr);
}
example
{  "EXAMPLE:"; echo = 2;
    ring rr=0,x,dp;
    list L = 150,3204,79272,2245968;
    list Tr = BerlekampMassey(L,1);
    Tr[1];
    factorize(Tr[1]); //not linearly factored
    list L1 = 70411680, 2352815424, 81496927872;
    Tr = BerlekampMassey(L1,1,Tr); // increase the length of L by size(L1)
    Tr[1];
    factorize(Tr[1]); //linearly factored and has distinct roots
    Tr[2]; //the length of the sequence required to generate Tr[1]
}

///////////////////////////////////////////////////////////////////////////////
// ++++++++++++++++++ Sparse Multivariate Interpolation +++++++++++++

static proc find_monomials(bigint B, list L,int n)
{
    /* return monomial(s) represented by B w.r.t. L where L[i] -> var(i)
     * L is list of primes */
    int nr=size(L);
    poly f=1;
    list l;
    for(int j=n+1;j<=n+nr;j++)
    {
        l = p_adic_valuation(B,L[j-n]);
        f = f*var(j)**(l[1]);
        B = l[2];
    }
    return(f);
}

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

static proc p_adic_valuation(bigint B, bigint p)
{
    // return the exponent j of the greatest power of p that divides B
    int j=-1;
    bigint H=1;
    while(1)
    {
        j++;
        H = H*p;
        if((B mod H)!=0)
        {
            return(list(j,(B div (H div p))));
        }
    }
}

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

static proc rootsofpoly(poly f, int n)
{
    // return  roots of f only for linearly factorizable polynomal f
    if(n==0){n=1;}
    ideal J=factorize(f,1);
    list L;
    for(int i=1;i<=ncols(J);i++)
    {
        L[i]=bigint(-1*subst(J[i],var(n),0));
    }
    return(L);
}

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

static proc generate_tVandermondeMatrix(list lroot, list lprime, int n, list #)
{
    /*
     * return list L of monomials w.r.t. list of primes lprime and
     * transposed Vandermonde matrix M
     */
    int nr=size(lroot);
    int i,j,mr;
    int k=1;
    mr = nr;
    if(size(#)>0)
    {
        k = #[1];
        mr = mr+k-1;
    }
    list L;
    matrix M[nr][nr];
    for(j=1;j<=nr;j++)
    {
        L[j]=find_monomials(lroot[j],lprime,n);
        for(i=k;i<=mr;i++)
        {
            M[i-k+1,j]= lroot[j]**i;
        }
    }
    return(list(L,M));
}

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

proc sparseInterpolation(poly Br,list La,list lpr,int n, list #)
"USAGE:  sparseInterpolation(Br, La, lpr, n[, m]); Br poly, La list, lpr list, n int, m int
RETURN:  a polynomial B in the polynomial ring Q[var(n+1),...,var(n+size(lpr))]
         satisfying the relation La[i] = B(lpr[1]^i,...,lpr[size(lpr)]^i).
NOTE: The polynomial Br in Q[var(n)] is the minimal polynomial obtained by
      applying the SINGULAR command @ref{BerlekampMassey} to the sequence (La[j]),
      1<=j<=size(La). By default the exponent i starts from 1. However, if the optional
      parameter m>=0 is provided, then it starts from m.
      The list lpr must be a list of distinct primes.
SEE ALSO: BerlekampMassey
EXAMPLE: example sparseInterpolation; shows an example
"
{
    /* n is the number of variables of the base ring
     * Br the minimal polynomial of the sequence La
     */
    int na=size(La);

    // compute the roots of f using factorization algorithm

    list Lr=rootsofpoly(Br,n);

    // Monomials and transposed Vandermonde Matrix

    list F=generate_tVandermondeMatrix(Lr, lpr, n, #);

    //  compute the coefficients of the monomials in F[1]

    int nr=nrows(F[2]);
    list la=F[1]; // list of monomials
    matrix T[nr][1]=La[1..nr];
    matrix V=inverse(F[2]); // compute inverse of the matrix F[2]
    matrix Z=V*T; // the coefficient ci

    // the procedure ends here

    matrix C[1][nr]=la[1..nr];
    poly g=(C*Z)[1,1];
    return(g);
}
example
{   "EXAMPLE:"; echo = 2;
    ring rr=0,(x,y),dp;
    list lpr = 2,3; // assign 2 for x and 3 for y
    list La = 150,3204,79272,2245968,70411680, 2352815424, 81496927872;
    // La[i] = number(subst(f,y,lpr[1]^i,z,lpr[2]^i)); for f = x2y2+2x2y+5xy2 and i=1,...,7
    poly Br = BerlekampMassey(La,1)[1];
    Br;
    sparseInterpolation(Br,La,lpr,0); // reconstruct f default
    La = 97,275,793,2315,6817;
    // La[i] = number(subst(g,y,lpr[1]^i,z,lpr[2]^i)); for g = x+y and i=4,...,8
    Br = BerlekampMassey(La,1)[1];
    Br;
    sparseInterpolation(Br,La,lpr,0,4);
}

///////////////////////////////////////////////////////////////////////////////
// ++++++++++++++++++++++ univariate rational function reconstruction ++++

proc fareypoly(poly g, poly f,list #)
"USAGE:  fareypoly(f, g[, m]); f poly, g poly, m int
RETURN:  a list l where r/t (r:=l[1], t:=l[2]) is a univariate rational function
         such that r/t = g mod f, gcd(r,t)=gcd(f,t)=1 and deg(r) + deg(t) < deg(f)
NOTE: An optional parameter m can be provided to define the way how t is normalized.
      If m = 0 (default), then the leading coefficient of t is 1. Otherwise,
      assuming the polynomial t has a non-zero constant term, the procedure
      returns the uniquely determined rational function r/t where the constant term
      in t is equal to 1.
      If the ground ring has n variables and f and g are in a polynomial
      ring k[var(i)] (k is a field) for some i<=n, then the function r/t is returned
      as an element in k(var(i)).
      In positive characteristic, the condition r/t = g mod f may not be satisfied.
      The degree deg(f) of f must be higher than the degree deg(g) of g.
SEE ALSO: polyInterpolation, farey
EXAMPLE: example fareypoly; shows an example
"
{
    int normalize_constant_term = 0; // default
    if(size(#) > 0 && typeof(#[1]) == "int")
    {
        normalize_constant_term = #[1];
    }
    poly r1,r2,r3,t1,t2,q_m,r_m,t_m,q1;
    q_m = 1;
    if(g==0)
    {
        return(list(poly(0),poly(1)));
    }
    if(2*deg(g)<deg(f))
    {
        // the degree of f is large enough
        return(list(g,poly(1)));
    }
    number h=number(1)/lu(g);
    r2=g*h;
    r1=f/lu(f);
    t1=0;
    t2=h;
    list ls,l1,l,T;
    int i=0;
    // a modified while loop in the Extended Euclidean algorithm
    while(r2!=0)
    {
        i++;
        ls=division(r1,r2);
        r3=r2;
        q1=ls[1][1,1]; // quotient
        h=number(1)/lu(ls[2][1]);
        r2=ls[2][1]*h; // remainder times h
        r1=r3;
        r3=t2;
        t2=(t1-q1*t2)*h;
        t1=r3;
        /***** find a quotient q_m whose degree is maximal and polynomials r_m & t_m
        correspond to q_m *********************************************************/
        if( deg(q1) > deg(q_m))
        {
            q_m=q1;
            r_m=r1;
            t_m=t1;
        }
    }
    if(deg(q_m)==1)
    {
        return(list(g,poly(1))); //trivial solution
    }
    else
    {
        poly vd = gcd(r_m,t_m);
        if(vd!=1)
        {
            //gcd condition is not satisfied;
            t_m = t_m/vd;
            r_m = r_m/vd;
            if(normalize_constant_term)
            {
                // here we normalize only for internal use
                number ut=number(1)/lu(t_m[size(t_m)]);
                return(list(ut*r_m,ut*t_m));
            }
            number ut=lu(t_m);
            return(list(r_m/ut,t_m/ut));
        }
        else
        {
            if(normalize_constant_term)
            {
                // here we normalize only for internal use
                number ut=number(1)/lu(t_m[size(t_m)]);
                return(list(ut*r_m,ut*t_m));
            }
            number ut=lu(t_m);
            return(list(r_m/ut,t_m/ut));
        }
    }
}
example
{   "EXAMPLE:"; echo = 2;
    ring rr=23,x,dp;
    poly g = 10x5-5x4+3x3+3x2-x-11;
    poly f = x6+2x5-9x4+x3-9x2+7x+7;
    fareypoly(g,f);
    fareypoly(g,f,1);

    ring R = 0, x,dp;
    poly g = (24/1616615)*x6-(732/1616615)*x5+(9558/1616615)*x4-(14187/323323)*x3+
    (1148101/1616615)*x2+(4089347/1616615)*x+547356/230945;
    poly f = x7-28x6+322x5-1960x4+6769x3-13132x2+13068x-5040;
    fareypoly(g,f);
    fareypoly(g,f,1);

    ring r = (499,a),x,dp;
    number N = (-113a4+170a3-29a2+226a+222)/(a7-56a6+114a5+144a4+171a3-64a2+192a);
    poly h1 = x4+(-55a5-18a4-141a3+233a2+66a-40)/(a4-28a3+40a2-2a+210)*x3;
    poly h2 = (107a6-221a5-68a4-93a3+112a2-54a+216)/(a4-28a3+40a2-2a+210)*x2;
    poly h3 = (-53a7+214a6+27a5+12a4+15a3+60a2-167a-83)/(a4-28a3+40a2-2a+210)*x;
    poly h4 = (10a6-75a5+47a4+246a3-20a2-217a+196)/(a4-28a3+40a2-2a+210);
    poly g = N*(h1+h2+h3+h4);
    poly f = x5+(-2a-119)*x4+(a2+237a+3437)*x3+(-118a2-6756a-29401)*x2+
    (3319a2+55483a+26082)*x+(-26082a2-26082a);
    fareypoly(g,f);
}

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

static proc lu(poly f)
{
    // lu leading unit of f
    if(f!=0)
    {
        return(leadcoef(f));
    }
    return(1);
}

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

static proc list_coef_index(list L, int idx1, int idx2,int lmt)
{
    // return list of leadcoef of the list L w.r.t. idx1, idx2,lmt
    ideal K;
    list lv;
    for(int j=1;j<=lmt;j++)
    {
        K = L[j];
        lv[j] = leadcoef(K[idx1][idx2]);
    }
    return(lv);
}

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

static proc scalIdeal(ideal I)
{
    //clear the denominators in the ideal I
    int t=ncols(I);
    if(size(I)==0)
    {
        return(I);
    }
    else
    {
        for(int i=1;i<=ncols(I);i++)
        {
            I[i]=cleardenom(I[i]);
        }
    }
    return(I);
}

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

static proc evaluate_f_at_given_points(poly f,list shft,int n)
{
    // evaluate f at var(n+k) = bigint(shft[k])**j for each j
    if(deg(f)==0 or f==0)
    {
        return(list(poly(0)));
    }
    int sr=size(shft);
    int k;
    poly g;
    for(k=n+1;k<=n+sr;k++)
    {
        f=subst(f,var(k),var(n)*var(k)+shft[k-n]);
    }
    matrix M = coeffs(f,var(n));
    list L = M[(nrows(M)-1)..1,1];
    return(L);
}

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

static proc evaluatef_ataprime_power(poly f, list prm, int n, int in, int st)
{
    // evaluate f at var(n+k) = bigint(prm[k])**j for in <= j <= st
    list L;
    poly v;
    int k;
    for(int j=in;j<=st;j++)
    {
        v=subst(f,var(n+1),bigint(prm[1])**j);
        for(k=2;k<=size(prm);k++)
        {
            v=subst(v,var(k+n),bigint(prm[k])**j);
        }
        L= L + list(number(v));
    }
    return(L);
}

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

static proc Add_poly_list(list lup)
{
    // return a list of the sum of two polynomials in the list lup
    if(size(lup[2])==0)
    {
        return(lup[1]);
    }
    else
    {
        list l1 = lup[1];
        list l2 = lup[2];
        for(int j=1;j<=size(l1);j++)
        {
            l1[j] = l1[j]+l2[j];
        }
    }
    return(l1);
}

///////////////////////////////////////////////////////////////////////////////
// subtract L[1][i]-L[2][i]

static proc SubList(list L)
{
    if(size(L[2])==0)
    {
        return(L[1]);
    }
    else
    {
        list l1 = L[1];
        list l2 = L[2];
        for(int j=1;j<=size(l1);j++)
        {
            l1[j] = l1[j]-l2[j];
        }
    }
    return(l1);
}

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

static proc Add_the_list_farey(list L)
{
    /* the procedure returns list of the denominators and  the numerators of the
    rational functions */
    list lst,lyt,Yt,lm;
    int i,j;
    for(j=1;j<=size(L[1]);j++)
    {
        lst = L[1];
        if(size(lst[j])!=0)
        {
            lyt = lst[j];
            for(i=2;i<=size(L);i++)
            {
                lst = L[i];
                lyt = Add_two_lists(lyt,lst[j]);
            }
        }
        Yt[j]= lyt;
        lyt = list();
    }
    return(Yt);
}
example
{   "EXAMPLE:"; echo = 2;
    ring rr=0,y,dp;
    list l1 = list(list(-2y-5/6,poly(1))); // (-2y-5/6)/(1) is a rational function
    list l2 = list(list(poly(-7/6),7/2y+1));
    list m1 = list(list(-4y-5/6,poly(1)));
    list m2 = list(list(poly(-7/6),21/2y+1));
    list L = list(l1,l2),list(m1,m2);
    Add_the_list_farey(L);
}

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

static proc Add_list_of_list(list l1,list l2, int m)
{
    // the procedure returns a list containing list of polynomials w.r.t. m
    list lst,lyt,Yt,lm;
    if(size(l1)!=size(l2))
    {
        ERROR("wrong size: sizes of lists do not coincide");
    }
    int i,j;
    for(j=m;j<=size(l1);j++)
    {
        if(size(l1[j])!=0)
        {
            for(i=1;i<=size(l1[j]);i++)
            {
                l1[j][i][1] = l1[j][i][1] + l2[j][i][1];
                l1[j][i][2] = l1[j][i][2] + l2[j][i][2];
            }
        }
    }
    return(l1);
}
example
{
    "EXAMPLE:"; echo = 2;
    ring rr=0,y,dp;
    list l2,m2;
    l2[1] = list(list(-2y-5/6,y+1));
    m2[1] = list(list(-4y-5/6,y-5));
    l2[2] = list(list(-2y-5/6,2y+3));
    m2[2] = list(list(-4y-5/6,3y-7));
    Add_list_of_list(l2,m2,1);
}
///////////////////////////////////////////////////////////////////////////////

static proc Add_two_lists(list l1, list l2)
{
    // the procedure returns a list containing list of polynomials
    int im=size(l1);
    int k,i;
    list l,m;
    for(k=1;k<=im;k++)
    {
        for(i=1;i<=2;i++)
        {
            if(typeof(l1[k][i])=="poly")
            {
                l[i] = list(l1[k][i]) + list(l2[k][i]);
            }
            else
            {
                l[i] = l1[k][i] + list(l2[k][i]);
            }
        }
        m[k] = l;
        l = list();
    }
    return(m);
}
example
{
    "EXAMPLE:"; echo = 2;
    ring rr=0,y,dp;
    list l2,m2;
    l2[2] = list(list(-2y-5/6,poly(1)));
    m2[2] = list(list(-4y-5/6,poly(1)));
    Add_two_lists(l2[2],m2[2]);
}

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

static proc arrange_list_first(list L)
{
    // arrange a given list
    list T,TT;
    int j,u,l1,l2;
    l1=size(L[1]);
    l2=size(L);
    for(u=1;u<=l1;u++)
    {
        for(j=1;j<=l2;j++)
        {
            TT[j]=L[j][u];
        }
        T[u]=TT;
    }
    return(T);
}
example
{   "EXAMPLE:"; echo = 2;
    ring rr=0,y,dp;
    list l1 = 1,2,3;
    list l2 = 1,4,9;
    list l3 = 1,8,27;
    list l4 = 1,16,81;
    list L = l1,l2,l3,l4;
    arrange_list_first(L);
}

///////////////////////////////////////////////////////////////////////////////
// return list of leadcoef of I w.r.t. i, I is of type list of list

static proc return_coef_indx(def I, int i)
{
    list l;
    for(int j=1;j<=size(I);j++)
    {
        l[j] = leadcoef(I[j][i]);
    }
    return(l);
}

///////////////////////////////////////////////////////////////////////////////
// return list of leadcoef of I of size k w.r.t. i, I is of type list of list

static proc return_coef_indx_wrtk(def I, int i, int k)
{
    list l;
    for(int j=1;j<=k;j++)
    {
        l[j] = leadcoef(I[j][i]);
    }
    return(l);
}

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

static proc arrangeListofIdeals(list L, list #)
{
    /* the procedure returns list of coeffcients in the given list of ideals. The
     * coeffcients are correspond to polynomials whose leading monomials are the same.
     * Moreover, if optional parameter # is given, then it also returns a list
     * of monomials in one of these ideals
    */
    ideal Tr = L[1];
    L = arrange_list_first(L);
    list l1,l2,ln,ld, J,ln1,ld1,Ld;
    int i,j, k;
    number N;
    poly f;
    for(i=1;i<=size(L);i++)
    {
        J = L[i];
        if(size(J[1])>1)
        {
            for(j=2;j<=size(J[1]);j++)
            {
                for(k=1;k<=size(J);k++)
                {
                    N = number(leadcoef(J[k][j]));
                    ln[k] = numerator(N);
                    ld[k] = denominator(N);
                }
                l1=ln,ld;
                l2[j-1] = l1;
                ln = list();
                ld = list();
            }
        }
        Ld[i] = l2;
        l2 = list();
    }
    if(size(#)>0)
    {
        list Zr = list_all_monom(Tr);
        return(list(Ld,Zr));
    }
    return(Ld);
}
example
{
    "EXAMPLE:"; echo = 2;
    ring rr=(0,a),y,dp;
    ideal I1 = y2+(3a2+6a+5)/(a+2)*y+1, y-6a;
    ideal I2 = y2+(3a2+5a+7)/(a+4)*y+1, 7y-17a;
    list L = I1,I2;
    arrangeListofIdeals(L);
    arrangeListofIdeals(L,1);
}

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

static proc Add_the_shift_and_evaluate_new(ideal J,list pr, list shft, int i)
{
    // evaluate J at a given point
    int k;
    number Nm;
    ideal Jc=J;
    int sc = size(pr);
    for(k=1;k<=sc;k++)
    {
        Nm = (par(sc+1)*number(pr[k])**i)+shft[k];
        Jc =subst(Jc,par(k),Nm);
    }
    return(Jc);
}

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

static proc generate_uniRationalFunctions(ideal I, list pr, list shift, int in,
    int fn, string Command, list JL,list #)
{
    // generate a set of ideals whose coefficients are univariate rational functions
    def Gt = basering;
    int i,i1;
    int tp = 0;
    if(size(#)>0){tp = 1; ideal Jc = #[1];}
    list L;
    for(i=in;i<=fn;i++)
    {
        L = L + list(Add_the_shift_and_evaluate_new(I, pr, shift, i));
    }
    list rl = ringlist(Gt);
    rl[1][2] = list("AXVR");
    def St = ring(rl);
    setring St;
    list L = imap(Gt, L);
    list optL = imap(Gt,JL);
    list T;
    int tmp;
    int c_z = size(L);
    for(i1=1;i1<=c_z;i1++)
    {
        if(Command == "slimgb")
        {
            task tk(i1) = "slimgb", list(L[i1]);
        }
        else
        {
            task tk(i1) = "Ffmodstd::ffmodStdOne",list(L[i1],list(optL,1));
        }
    }
    startTasks(tk(1..c_z));
    waitAllTasks(tk(1..c_z));
    for(i1 = 1;i1 <= c_z; i1++)
    {
        T[i1] = getResult(tk(i1));
        killTask(tk(i1));
        kill tk(i1);
    }
    if(tp)
    {
        ideal Jc = imap(Gt,Jc);
        T = list(Jc)+T;
    }
    setring Gt;
    list M = imap(St,T);
    M = normalize_LiftofIdeal(M);
    return(M);
}

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

static proc normalize_LiftofIdeal(list L)
{
    // normalize each ideal L[j]
    for(int j=1;j<=size(L);j++)
    {
        L[j] = normalize(L[j]);
    }
    return(L);
}

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

static proc normalize_constTerm(poly g, poly f)
{
    /* Assumption: f has a non-zero constant. Then the procedure
    returns a pair (g/N, f/N) where N is the constant term of f */
    number N = leadcoef(f[size(f)]);
    g = g/N;
    f = f/N;
    return(list(g,f));
}

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

static proc normalize_constTermAll(list M2, int kk)
{
    /*
      Assumption: each polynomial in the list has a non-zero constant term. Then
      the procedure normalize_constTermAll normalizes the constatnt term of
      each polynomial in this list starting from M2[kk] where kk is a
      positive integer less than or equal to the size of M2
    */
    int i,j,k;
    list l,l1,l2,l3,l4;
    if(kk!=1)
    {
        matrix M[1][kk-1]=0;
        l4 = M[1,1..kk-1];
    }
    for(i=kk;i<=size(M2);i++)
    {
        for(j=1;j<=size(M2[i]);j++)
        {
            for(k=1;k<=size(M2[i][j][1]);k++)
            {
                l =  normalize_constTerm(M2[i][j][1][k], M2[i][j][2][k]);
                l1 = l1+list(l[1]);
                l2 = l2+list(l[2]);
            }
            l3[j] = list(l1,l2);
            l1 = list();
            l2 = list();
        }
        l4[i] = l3;
        l3 = list();
    }
    return(l4);
}

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

static proc stdoverFF(ideal I, list pr, list shft,string Command,list Zr, list JL)
{
    // return std of I with high probability using sparse rational interpolation
    int fn=13;
    def R_1=basering;
    int n,pa;
    n = nvars(R_1);
    pa = npars(R_1);
    list rl=ringlist(R_1);
    list la=rl[1][2];
    list m=rl[2];
    m[(n+1)..(n+pa)]=la[1..pa];
    rl[2]=m;
    rl[1]=rl[1][1];
    rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
    rl[3][size(rl[3])-1]=list("lp",pa);
    def R_2=ring(rl);
    setring R_2;
    setring R_1;
    int in = 2;
    // with respect to given bounds, generate univariate rational functions
    list M2 = generate_uniRationalFunctions(I, pr, shft, in,fn, Command,JL[1],JL[2]);
    M2 = arrangeListofIdeals(M2);
    poly Gn;
    list uM1;
    ideal J;
    setring R_2;
    list M2 = imap(R_1,M2);
    M2 = normalize_constTermAll(M2,1); // starts from 1
    list Zr = imap(R_1, Zr); // list of monomials in the std of I
    int u_w,i1,i2,i3;
    poly Gn,pn,pl,plm;
    list l1,l2,l3,bp,lup,lk1,lk2,lk3,py,l3n;
    ideal J;
    list #;
    while(u_w < size(Zr))
    {
        u_w++;
        if(size(Zr[u_w])>1)
        {
            setring R_1;
            Gn = Zr[u_w][1];
            setring R_2;
            /******** start lifting elements from Q to Q(var(1),.., var(pa)) *****/
            l1 = M2[u_w];
            for(i1=1;i1<=size(l1);i1++)
            {
                for(i2=1;i2<=2;i2++)
                {
                    // the procedure first lifts the numerator and then the denominator
                    for(i3=1;i3<=size(l1[i1][i2][1]);i3++)
                    {
                        l3 = return_coef_indx(l1[i1][i2],i3);
                        l3 = l3,lk1;
                        l3 = SubList(l3); // adjust the coefficients
                        // early termination for BerlekampMassey algorithm(BMA)
                        bp = BerlekampMassey(l3,n);
                        if(size(bp)==2)
                        {
                            // at this step BerlekampMassey algorithm terminated early
                            if(bp[1]==1)
                            {
                                lup = lk3,list(); // if l3 = 0, ..., 0
                            }
                            else
                            {
                                l3 = l3[1..bp[2]];
                                pn = sparseInterpolation(bp[1],l3,pr,n);
                                pl = pn+pl;
                                lup = evaluate_f_at_given_points(pn, shft, n);
                                lup = lup,lk3;
                            }
                        }
                        else
                        {
                            /* elements in the sequence l3[i] are not large enough.
                             * We thus add more elements to l3 and continue with the
                             * early termination strategy until size(bp)=2 where bp[2]
                             * is the length of the sequence that the early temination
                             * of BMA requires
                             */
                            l3n = bp[6];
                            while(1)
                            {
                                in = fn+1;
                                fn = 2*in;
                                setring R_1;
                                /* add more points to l3 by generating rational
                                functions*/
                                uM1 = generate_uniRationalFunctions(I, pr, shft, in,
                                                                    fn, Command,JL[1]);
                                uM1 = arrangeListofIdeals(uM1);
                                setring R_2;
                                list uM2 = imap(R_1, uM1);
                                uM2 = normalize_constTermAll(uM2,u_w); // u_w optional
                                M2 = Add_list_of_list(M2,uM2,u_w);
                                l3 = return_coef_indx(uM2[u_w][i1][i2],i3);
                                kill uM2;
                                if(i3>1)
                                {
                                    lk1 = evaluatef_ataprime_power(plm, pr,n,in,fn);
                                    l3 = l3,lk1;
                                    l3 = SubList(l3);
                                }
                                # = bp;
                                l3n = l3n+l3;
                                bp = BerlekampMassey(l3,n,#);
                                if(size(bp)==2)
                                {
                                    // at this step BerlekampMassey algorithm terminated early
                                    l3n = l3n[1..bp[2]];
                                    pn = sparseInterpolation(bp[1],l3n,pr,n);
                                    pl = pn+pl;
                                    lup = evaluate_f_at_given_points(pn,shft,n);
                                    lup = lup,lk3;
                                    break;
                                }
                            }
                            l1 = M2[u_w];
                        }
                        if(i3 < size(l1[i1][i2][1]))
                        {
                            // unshift the shifted parameters see also SubList
                            lk3 = Add_poly_list(lup);
                            plm = lk3[1];
                            lk1 = evaluatef_ataprime_power(lk3[1], pr, n,1,fn);
                            lk3 = delete(lk3,1);
                        }
                    }
                    py[i2]=pl;
                    pl=0;
                    lk1 = list();
                    lk2 = list();
                    lk3 = list();
                }
                setring R_1;
                list H = imap(R_2,py);
                // numerator H[1] & denominator H[2] are recovered for each i1=1,2,...
                Gn = Gn + (H[1]/H[2])*Zr[u_w][i1+1];
                kill H;
                setring R_2;
            }
        }
        else
        {
            setring R_1;
            Gn = Zr[u_w][1];
            setring R_2;
        }
        setring R_1;
        Gn = cleardenom(Gn);
        J[u_w] = Gn;
        setring R_2;
    }
    setring R_1;
    return(J);
}


///////////////////////////////////////////////////////////////////////////////
// +++++++++++++++++ std for one parameter begins here +++++++++++++++++++

static proc test_fmodI(poly f,ideal I)
{
    // test whether f in I or not
    ideal If=f;
    attrib(If,"isSB",1);
    if(size(reduce(I,If))!=ncols(I))
    {
        return(0);
    }
    return(1);
}

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

static proc  choose_evaluation_points(ideal I, ideal cI,int n1,int nt, int i_s)
{
    /*
     * I is given ideal
     * cI is an output from Testlist_all
     * n1 is number of variables
     * nt number of prime ideals
     * i_s is an integer
     */
    int j,ss,si;
    list m,n;
    poly f;
    poly g=1;
    while(1)
    {
        ss++;
        if(i_s==0)
        {
            ERROR("no more points");
        }
        f=var(n1)-i_s;
        if(test_fmodI(f,cI)==1)
        {
            j=j+1;
            // specialize var(n1) with i_s
            n=n+list(subst(I,var(n1),i_s));
            m[j]=i_s;
        }
        i_s = i_s-1;
        if(j==nt)
        {
            n=n,m;
            return(n);
        }
    }
}

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

static proc firststdmodp(ideal I,ideal cI,int in_value)
{
    /* return std FJ of I together with the coefficients, as a list of univariate
    rational functions, of FJ */
    def St=basering;
    // define a new ring
    int nx=nvars(St);
    int nr=nx-1;
    list ly=ringlist(St);
    list L,l1,l3,l4,l5,A1,l2;
    l1=ly[2];
    l2=l1[1..nr];
    L[2]=l2;
    l3=l1[(nr+1)..nx];
    l5=ly[3];
    l4=l5[2];
    A1=ly[1],l3,list(l4),ly[4];
    L[1]=A1;
    l5=delete(l5,2);
    L[3]=l5;
    L[4]=ly[4];
    def Gt=ring(L);
    kill L,l1,l3,l4,l5,A1,l2;
    setring Gt;
    setring St;
    int i,k1,k2;
    list T,T1,L1,L2,L3,m_l;
    int r_d = char(basering)-10000000; // choose an integer r_d which is d/t from
    // char(basering)
    list l_1 =  choose_evaluation_points(I,cI, nx,in_value,r_d);
    list lus = l_1[2];
    l_1 = l_1[1];
    m_l[1] = std(l_1[1]);
    if(size(m_l[1])==1 and deg(m_l[1][1])==0)
    {
        return(list(ideal(1)));
    }
    // compute std in parallel w.r.t. distinct r_d
    for(k1 = 2; k1<= in_value; k1++)
    {
        task tk(k1-1) = "std", list(l_1[k1]);
    }
    startTasks(tk(1..in_value-1));
    waitAllTasks(tk(1..in_value-1));
    for(k1 = 1;k1<=in_value-1;k1++)
    {
        m_l = m_l + list(getResult(tk(k1)));
        killTask(tk(k1));
        kill tk(k1);
    }
    r_d = lus[in_value]-1;
    // DeleteUnluckyEvaluationPoints
    list indices = Modstd::deleteUnluckyPrimes_std(m_l);
    for(i = size(indices); i > 0; i--)
    {
        m_l = delete(m_l, indices[i]);
        lus = delete(lus, indices[i]);
    }
    m_l = normalize_LiftofIdeal(m_l);
    in_value = size(lus);
    int ug;
    for(int cZ =1;cZ<=ncols(m_l[1]);cZ++)
    {
        if(size(m_l[1][cZ])>1)
        {
            ug = cZ;
            break;
        }
    }
    if(!ug)
    {
        return(list(m_l[1]));
    }
    // dense univariate rational interpolation
    list lev = list_coef_index(m_l,ug,2,in_value-1);
    list L = polyInterpolation(list(lus[1..in_value-1]),lev,nx);
    poly G1 = L[2];
    poly G2 = G1*(var(nx)-lus[in_value-1]);
    poly DR,NR;
    list fry = fareypoly(L[1],G2);
    DR = fry[2];
    NR = fry[1];
    list # = list(NR)+list(DR)+L;
    list M = m_l;
    ideal J = m_l[1];
    list Zr = list_all_monom(J);
    lus = lus[in_value];
    int n_z,m_x;
    list M1,M2,Fr;
    setring Gt;
    poly g_t;
    list Zr = imap(St,Zr);
    ideal FJ = imap(St,J);
    setring St;
    for(k1=ug;k1<=ncols(J);k1++)
    {
        n_z = size(J[k1]);
        setring Gt;
        g_t = Zr[k1][1];
        setring St;
        if(n_z>1)
        {
            for(k2=2;k2<=size(J[k1]);k2++)
            {
                if(size(#)<2)
                {
                    lev = list_coef_index(M,k1,k2,size(lus));
                }
                else
                {
                    lev = list_coef_index(list(M[size(M)]), k1, k2,1);
                }
                T1 = fareypolyEarlyTermination(I,cI,M,nx,lus,lev, k1, k2, r_d,#);
                // save maximum number of data m_x that are used in each lifting
                if(m_x < T1[3][2]){ m_x = T1[3][2];}
                L1[k2-1] = T1[3];
                Fr = T1[4];
                setring Gt;
                list Fr = imap(St,Fr);
                // add the coeffcients to g_t which are already lifted
                g_t = g_t + (Fr[1]/Fr[2])*Zr[k1][k2];
                kill Fr;
                setring St;
                M = T1[1];
                lus = T1[2];
                r_d = T1[5];
                # = nx;
            }
            L2[k1] = L1;
            L1 = list();
        }
        setring Gt;
        g_t = cleardenom(g_t);
        FJ[k1]= g_t;
        setring St;
    }
    ideal FJ = imap(Gt,FJ);
    FJ = normalize(FJ);
    L2 = FJ,L2,Zr,m_x;
    return(L2);
}

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

static proc list_all_monom(ideal T)
{
    // list all monomials in T
    int nr=ncols(T);
    list L,E;
    list l;
    int i,j;
    poly f;
    for(j=1;j<=nr;j++)
    {
        f=T[j];
        for(i=1;i<=size(f);i++)
        {
            l[i]=leadmonom(f[i]);
        }
        L[j]=l;
        l=E;
    }
    return(L);
}

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

static proc fareypolyEarlyTermination(ideal I, ideal cI, list M, int nx,
                       list lus, list lev, int k1, int k2, int r_d, list #)
{
    /* the early termination version of the farey rational funtion map for
    polynomials */
    ideal Id;
    list Tr,fr,L,l_1,m_l;
    int ksz,us,sz, deg_fr, mx_data,i,k11;
    number ev;
    poly G1,G2,DR,NR;
    int in_value = 10;
    list indices;
    if(size(#)==1)
    {
        // dense univariate rational interpolation
        sz = size(lus);
        L = polyInterpolation(list(lus[1..sz-1]),list(lev[1..sz-1]),#);
        fr = fareypoly(L[1],L[2]);
        NR = fr[1];
        DR = fr[2];
        # = L;
        L = polyInterpolation(lus[sz],lev[sz],nx,#);
        fr = fareypoly(L[1],L[2]);
    }
    else
    {
        // update by adding more evaluation points
        NR = #[1];
        DR = #[2];
        # = #[3..size(#)];
        L = polyInterpolation(lus,lev,nx,#);
        sz=size(lus);
        fr = fareypoly(L[1],L[2]);
    }
    us = r_d;
    while(fr[1]!=NR or fr[2]!=DR)
    {
        l_1 =  choose_evaluation_points(I, cI, nx, in_value, us);
        lus = l_1[2];
        l_1 = l_1[1];
        // compute std in parallel w.r.t. lus
        for(k11 = 1; k11 <= in_value; k11++)
        {
            task tk(k11) = "std", list(l_1[k11]);
        }
        startTasks(tk(1..in_value));
        waitAllTasks(tk(1..in_value));
        for(k11 = 1;k11<=in_value;k11++)
        {
            m_l[k11] = getResult(tk(k11));
            killTask(tk(k11));
            kill tk(k11);
        }
        us = lus[size(lus)]-1;
        // DeleteUnluckyEvaluationPoints
        indices = Modstd::deleteUnluckyPrimes_std(m_l);
        for(i = size(indices); i > 0; i--)
        {
            m_l = delete(m_l, indices[i]);
            lus = delete(lus, indices[i]);
        }
        sz = size(m_l);
        m_l = normalize_LiftofIdeal(m_l);
        lev = list_coef_index(m_l,k1,k2,sz);
        # = L;
        // dense univariate rational interpolation
        L = polyInterpolation(list(lus[1..sz-1]),list(lev[1..sz-1]),nx,#);
        fr = fareypoly(L[1],L[2]);
        NR = fr[1];
        DR = fr[2];
        M = M + m_l;
        # = L;
        L = polyInterpolation(lus[sz],lev[sz],nx,#);
        fr = fareypoly(L[1],L[2]);
        if(fr[1]==NR and fr[2]==DR)
        {
            deg_fr = deg(fr[1]);
            mx_data = deg_fr + deg(fr[2])+2;
            return(list(M,L[3],list(deg_fr, mx_data),fr,us));
            break;
        }
    }
    deg_fr = deg(fr[1]);
    mx_data = deg_fr + deg(fr[2])+2;
    return(list(M,L[3],list(deg_fr, mx_data), fr,us));
}

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

static proc  RecoverCoeffsForAFixedData(list stdResults,list distElmnt,list maxData,
                                       list Zr)
{
    /* here a bound on the number of interpolation points is known and, hence,
     * w.r.t this bound we apply rational interpolation algorithm to obtain
     * a set of polynomials over a function field modulo p. stdResults is a list
     * of std over F_p*/
    def St=basering;
    // define a new ring
    int nx=nvars(St);
    int nr=nx-1;
    list ly=ringlist(St);
    list L,l1,l3,l4,l5,A1,l2;
    l1=ly[2];
    l2=l1[1..nr];
    L[2]=l2;
    l3=l1[(nr+1)..nx];
    l5=ly[3];
    l4=l5[2];
    A1=ly[1],l3,list(l4),ly[4];
    L[1]=A1;
    l5=delete(l5,2);
    L[3]=l5;
    L[4]=ly[4];
    def Gt=ring(L);
    kill L,l1,l3,l4,l5,A1,l2;
    setring Gt;
    setring St;
    int i,k1,k2,n_z;
    list M = stdResults;
    ideal J = M[1];
    if(J[1]==1)
    {
        return(list(ideal(1)));
    }
    int ug;
    for(int cZ =1;cZ<=ncols(J);cZ++)
    {
        if(size(J[cZ])>1)
        {
            ug = cZ;
            break;
        }
    }
    if(!ug)
    {
        return(list(J));
    }
    setring Gt;
    poly g_t;
    list Zr = imap(St,Zr);
    ideal FJ = imap(St,J);
    setring St;
    for(k1=ug;k1<=ncols(J);k1++)
    {
        n_z = size(J[k1]);
        setring Gt;
        g_t = Zr[k1][1];
        setring St;
        if(n_z>1)
        {
            // dense univariate rational interpolation in parallel
            list TL = interpolation_farey_lift_parallel(maxData, M, distElmnt, n_z,
                                                       nx, k1);
            setring Gt;
            list LT = imap(St,TL);
            for(k2=2;k2<=n_z;k2++)
            {
                g_t = g_t + (LT[k2][1]/LT[k2][2])*Zr[k1][k2];
            }
            kill LT;
            setring St;
            kill TL;
        }
        setring Gt;
        g_t = cleardenom(g_t);
        FJ[k1]= g_t;
        setring St;
    }
    ideal FJ = imap(Gt,FJ);
    return(normalize(FJ));
}

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

static proc interpolation_farey_lift_parallel(list maxData, list M, list distElmnt,
                                              int m_sz,int nx, int k1)
{
    // return list of univariate rational functions
    list L;
    int k2;
    for(k2=2;k2<=m_sz;k2++)
    {
       task tk(k2) = "Ffmodstd::lift_interp_farey", list(maxData, M, distElmnt,
       nx, k1, k2);
    }
    startTasks(tk(2..m_sz));
    waitAllTasks(tk(2..m_sz));
    for(k2 = 2;k2 <= m_sz; k2++)
    {
         L[k2] = getResult(tk(k2));
         killTask(tk(k2));
         kill tk(k2);
    }
    return(L);
}

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

static proc lift_interp_farey(list maxData, list M, list distElmnt,
       int nx, int k1, int k2)
{
    // lift points in F_p to a rational function in F_p(t)
    list mD, lev, lus,Fr;
    mD = maxData[k1][k2-1];
    lev = list_coef_index(M,k1,k2,mD[2]);
    lus = distElmnt[1..mD[2]];
    Fr = NewtonInterpolationNormal(lus,lev,nx);
    Fr = fareypoly_wrt_deg_dg(Fr[1],Fr[2], mD[1]);
    return(Fr);
}

///////////////////////////////////////////////////////////////////////////////
// ++++++++++++++++++++++ univariate rational function reconstruction ++++

static proc fareypoly_wrt_deg_dg(poly g, poly f, int dg, list #)
{
    // find a rational function whose degree of numerator is equal to dg
    //system("pid");
    int const_int=0;
    if(size(#)>0 or typeof(#[1])=="int")
    {
        const_int = #[1];
    }
    poly r1,r2,r3,t1,t2,q_m,r_m,t_m,q1;
    q_m = 1;
    if(g==0)
    {
        return(list(poly(0),poly(1)));
    }
    if(2*deg(g)<deg(f))
    {
        // here the degree of f is large enough
        return(list(g,poly(1)));
    }
    number h=number(1)/lu(g);
    r2=g*h;
    r1=f/lu(f);
    t1=0;
    t2=h;
    list ls,l1,l,T;
    int i=0;
    // a modified while loop in the Extended Euclidean algorithm
    while(r2!=0)
    {
        i++;
        ls=division(r1,r2);
        r3=r2;
        q1=ls[1][1,1];
        h=number(1)/lu(ls[2][1]);
        r2=ls[2][1]*h;
        r1=r3;
        r3=t2;
        t2=(t1-q1*t2)*h;
        t1=r3;
        if(deg(r1)==dg)
        {
            break;
        }
    }
    if(const_int)
    {
        number ut=number(1)/lu(t1[size(t1)]);
        return(list(ut*r1,ut*t1));
    }
    number ut=lu(t1);
    return(list(r1/ut,t1/ut));
}

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

static proc passThePreList(ideal I, ideal cI, int nva, int ma_x)
{
    /* return list of normalized std output and distict points used in the
     computations */

    list m_l,l_1,R_d, indices,G,Fd;
    int k11,v_r,i;
    int r_d = char(basering)-10000000; // choose an integer r_d d/f from char(basering)
    while(size(G) < ma_x)
    {
        l_1 =  choose_evaluation_points(I, cI,nva,ma_x,r_d);
        R_d = l_1[2];
        l_1 = l_1[1];
        // compute std in parallel
        for(k11 = 1; k11 <= ma_x; k11++)
        {
            task tk(k11) = "std", list(l_1[k11]);
        }
        startTasks(tk(1..ma_x));
        waitAllTasks(tk(1..ma_x));
        for(k11 = 1;k11 <= ma_x; k11++)
        {
            m_l[k11] = getResult(tk(k11));
            killTask(tk(k11));
            kill tk(k11);
        }
        indices = Modstd::deleteUnluckyPrimes_std(m_l);
        r_d = R_d[size(R_d)]-1;
        for(i = size(indices); i > 0; i--)
        {
            m_l = delete(m_l, indices[i]);
            R_d = delete(R_d, indices[i]);
        }
        m_l = normalize_LiftofIdeal(m_l);
        G = G + m_l;
        Fd = Fd + R_d;
    }
    return(list(G,Fd));
}

///////////////////////////////////////////////////////////////////////////////
/* return 0 if p divides any numerator or any denominator in the coefficients
of I */

static proc prime_test(int p, ideal I)
{

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

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

static proc prime_pass(int p, ideal I)
{
    //return a prime for which prime_test()==true
    // p must be prime
    int i,q;
    q=p;
    if(q<2)
    {
        ERROR("No more Primes");
    }
    else
    {
        while(1)
        {
            i++;
            if(prime_test(p,I)==1)
            {
                return(q);
            }
            q = prime(q-1);
        }
    }
}

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

static proc select_the_command(ideal I)
{
    /* the procedure select_the_command selects a command which finishes the
       computation first but if this command applied in a machine with a single
       core it returns (by default) the command ffmodStdOne
    */
    if(getcores() == 1)
    {
        def F = ffmodStdOne(I);
        return(list("ffmodStdOne", F));
    }
    else
    {
         list commands = list("Ffmodstd::ffmodStdOne", "slimgb");
         list args = list(list(I), list(I));
         list L = parallelWaitFirst(commands, args);
         if(typeof(L[1])!="none")
         {
            return(list("ffmodStdOne", L[1]));
         }
        else
        {
            return(list("slimgb", L[2]));
        }
    }
}

////////////////////////////////////////////////////////////////////////////////
// connect ffmodStdOne to modular

static proc modp_tran(list newL)
{
    /* compute a standard basis of I modulo p using a dense rational
     *   interpolation */
    int nva = nvars(basering);
    ideal I = newL[1];
    ideal cI = newL[2];
    list L = newL[3];
    list tedY = passThePreList(I, cI, nva, L[3]);
    list stdResults = tedY[1];
    list distElmnt = tedY[2];
    ideal F =  RecoverCoeffsForAFixedData(stdResults, distElmnt, L[1], L[2]);
    return(F);
}

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

static proc ffmodStdOne(ideal I,list #)
{
    // compute std of I using dense univariate rationa interpolation
    int tmp2 = 0;
    int tmp = 0;
    int tmp1 = 1;
    def G_t=basering;
    intvec opt = option(get);
    option(redSB);
    list pL;
    // optional parameters
    if(size(#)>0)
    {
        tmp = 1;
        if(size(#)==1)
        {
            pL = #[1];
        }
        else
        {
            pL = #[1];
            tmp2 = #[2];
        }
    }
    int n,pa,kr;
    n=nvars(G_t);
    pa=1;
    I = normalize(I); // make each element of I monic
    I=scalIdeal(I);
    list L=collect_coeffs(I);
    // define a new ring
    list rl=ringlist(G_t);
    list la=rl[1][2];
    list m=rl[2];
    m[(n+1)..(n+pa)]=la[1..pa];
    rl[2]=m;
    rl[1]=rl[1][1];
    rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
    rl[3][size(rl[3])-1]=list("lp",pa);
    def S_t=ring(rl);
    setring S_t;
    if(tmp){ list opL = imap(G_t, pL);}
    list L=imap(G_t,L);
    ideal I=imap(G_t,I);
    ideal cI=Testlist_all(L);
    if(size(cI)==0)
    {
        cI = 1;
    }
    int pr = 536870909;
    int paSS = prime_pass(pr, cI); // test whether pr is a suitable for cI
    pr = paSS;
    if(!tmp)
    {
        int in_value = 13; // initial number of prime ideals (or number of points)
        list lbr = ringlist(S_t);
        lbr[1] = pr;
        def rp = ring(lbr);
        setring(rp);
        // initial std computation using dense rational interpolation
        list rL = firststdmodp(imap(S_t,I),imap(S_t,cI), in_value);
        setring S_t;
        list rL = imap(rp,rL);
        if(size(rL)==1)
        {
            setring G_t;
            list EI = imap(S_t,rL);
            return(EI[1]);
        }
        /*
        * from the initial std we obtain bounds on the number of interpolation points
        * and degree of each numerator and denominator which we use for later
        * interpolation
        */
        rL= rL[2..size(rL)];
    }
    list newL,optionL;
    ideal J;
    if(tmp)
    {
        newL = I, cI, opL;
        // apply modular command from the Singular library parallel.lib
        if(tmp2)
        {
            // with final test
            J = modular("Ffmodstd::modp_tran", list(newL), primeTest_tran,
            Modstd::deleteUnluckyPrimes_std, pTest_tran, finalTest_tran, pr);
        }
        else
        {
            // without final test
            J = modular("Ffmodstd::modp_tran", list(newL), primeTest_tran,
            Modstd::deleteUnluckyPrimes_std, pTest_tran, pr);
        }
    }
    else
    {
        // apply modular command from the Singular library parallel.lib with final test
        newL = I, cI, rL;
        J = modular("Ffmodstd::modp_tran", list(newL),  primeTest_tran,
        Modstd::deleteUnluckyPrimes_std, pTest_tran, finalTest_tran, pr);
    }
    setring G_t;
    ideal J = imap(S_t, J);
    attrib(J, "isSB", 1);
    option(set,opt);
    if(!tmp)
    {
        optionL = imap(S_t,rL);
        return(list(J,optionL));
    }
    return(J);
}

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

static proc primeTest_tran(int p, list L)
{
    /*
     * test whether a prime p divides the denominator(s)
     * and leading coefficients of generating set of ideal
     */
    int i,j;
    ideal I = L[1][2]; // I = L[1][1]
    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);
}

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

static proc ideal_Inclusion(ideal I, ideal fareyresult)
{
    //return 1 if I is in fareyresult otherwise 0
    attrib(fareyresult,"isSB",1);
    int i;
    for(i = ncols(I); i >= 1; i--)
    {
        if(reduce(I[i],fareyresult,1)!= 0)
        {
            return(0);
        }
    }
    return(1);
}

///////////////////////////////////////////////////////////////////////////////
// pTest

static proc pTest_tran(string Command,list args, ideal fareyResult,int Testp)
{
    def St = basering;
    list lpr = ringlist(St);
    lpr[1] = Testp;
    def Stp = ring(lpr);
    setring Stp;
    list args = imap(St, args);
    ideal Jp = modp_tran(args[1]);
    setring St;
    ideal Jp = imap(Stp, Jp);
    list l = ringlist(St);
    list l1 = l[1], list(l[2][size(l[2])]), list(l[3][size(l[3])-1]), l[4];
    list l2 = delete(l[2],size(l[2]));
    list l3 = delete(l[3],2);
    list rl = l1,l2,l3,l[4];
    rl[1][1] = Testp;
    def Gt=ring(rl);
    setring Gt;
    ideal Jp = imap(St, Jp);
    ideal fry = imap(St,fareyResult);
    attrib(fry,"isSB",1);
    attrib(Jp,"isSB",1);
    if(ideal_Inclusion(Jp,fry))
    {
        // test if fry is in Jp
        if(size(reduce(fry,Jp))==0)
        {
            setring St;
            return(1);
        }
    }
    setring St;
    return(0);
}
///////////////////////////////////////////////////////////////////////////////

static proc finalTest(ideal I, ideal fareyresult)
{
    //return 1 if I included in fareyresult otherwise 0
    attrib(fareyresult,"isSB",1);
    int i;
    for(i = ncols(I); i >= 1; i--)
    {
        if(reduce(I[i],fareyresult,1)!= 0)
        {
            return(0);
        }
    }
    //return 1 if std(fareyresult) included in fareyresult otherwise 0
    attrib(fareyresult,"isSB",1);
    ideal J=std(fareyresult);
    if(size(reduce(J,fareyresult,1))!=0)
    {
        return(0);
    }
    return(1);
}

// =============== a procedure for one parameter ends here ==========

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

static proc finalTest_tran(string Command,list args, ideal fareyresult)
{
    //return 1 if I included in fareyresult otherwise 0
    def St = basering;
    list l = ringlist(St);
    list l1 = l[1], list(l[2][size(l[2])]), list(l[3][size(l[3])-1]), l[4];
    list l2 = delete(l[2],size(l[2]));
    list l3 = delete(l[3],2);
    list rl = l1,l2,l3,l[4];
    def Gt=ring(rl);
    setring Gt;
    list T = imap(St,args);
    ideal I = T[1][1];
    ideal fareyresult = imap(St, fareyresult);
    attrib(fareyresult,"isSB",1);
    int i;
    for(i = ncols(I); i >= 1; i--)
    {
        if(reduce(I[i],fareyresult,1)!= 0)
        {
            setring St;
            return(0);
        }
    }
    //return 1 if std(fareyresult) included in fareyresult otherwise 0
    attrib(fareyresult,"isSB",1);
    ideal J=std(fareyresult);
    if(size(reduce(J,fareyresult,1))!=0)
    {
        setring St;
        return(0);
    }
    setring St;
    return(1);
}

///////////////////////////////////////////////////////////////////////////////
// =========== the main procedure for multi parameters begins here =========

proc ffmodStd(ideal I)
"USAGE:  ffmodStd(I); I ideal
RETURN:  Groebner basis of I over an algebraic function field
SEE ALSO: nfmodStd
EXAMPLE: example ffmodStd; shows an example
"
{
    intvec opt = option(get);
    option(redSB); // to obtain the reduced standard basis
    def G_t=basering;
    int n,pa,kr;
    n=nvars(G_t);
    pa=npars(G_t);
    if(pa==0)
    {
        ERROR("the coefficient field is not rational function field");
    }
    if(size(I)==0)
    {
        return(ideal(0));
    }
    I = simplify(I,2);
    for(int hj=1;hj<=ncols(I);hj++)
    {
        if(deg(I[hj])==0)
        {
            return(ideal(1));
        }
    }
    // optional parameters
    if(pa==1)
    {
        // if the current ring has one parameter
        def GF = ffmodStdOne(I);
        if(size(GF)==1)
        {
            return(GF);
        }
        return(GF[1]);
    }
    I = normalize(I);
    I=scalIdeal(I);
    list L=collect_coeffs(I);
    // define a new ring
    list rl=ringlist(G_t);
    list la=rl[1][2];
    list m=rl[2];
    m[(n+1)..(n+pa)]=la[1..pa];
    rl[2]=m;
    rl[1]=rl[1][1];
    rl[3][size(rl[3])+1]=rl[3][size(rl[3])];
    rl[3][size(rl[3])-1]=list("lp",pa);
    def S_t=ring(rl);
    setring S_t;
    list L=imap(G_t,L);
    ideal I=imap(G_t,I);
    ideal cI =Testlist_all(L);
    if( size(cI) == 0)
    {
        cI = 1;
    }
    list shft = test_the_shift(cI,n,pa);
    int j;
    // shift the parameters
    for(j=1;j<=size(shft);j++)
    {
        cI = subst(cI, var(n+j), var(n+j) + shft[j]);
    }
    list pr=list_of_primes(pa); // list distinct primes
    // define a new ring
    setring G_t;
    rl = ringlist(G_t);
    rl[1][2][size(rl[1][2])+1] = "AXVR";
    def St = ring(rl);
    setring St;
    ideal Jc = Add_the_shift_and_evaluate_new(imap(G_t,I), pr, shft, 1);
    setring G_t;
    rl = ringlist(G_t);
    rl[1][2]= list("AXVR");
    def StA = ring(rl);
    setring StA;
    ideal Jc = imap(St,Jc);
    // make selection to use relitively fast command
    list Lcom = select_the_command(Jc);
    string Command = Lcom[1];
    def GF = Lcom[2];
    list Zr;
    if(typeof(GF)!="ideal")
    {
        Jc = GF[1];
        Zr = GF[2][2];
        Lcom = GF[2];
    }
    else
    {
        Jc = GF;
        Zr = list_all_monom(Jc);
        Lcom = poly(0);
    }
    Lcom = Lcom,Jc;
    if(ncols(Jc)==1 and Jc[1]==1)
    {
        setring G_t;
        return(ideal(1));
    }
    int cd;
    for(j=1;j<=ncols(Jc);j++)
    {
        if(size(Jc[j])>1)
        {
            cd = cd+1;
        }
    }
    if(!cd)
    {
        setring G_t;
        return(imap(StA, Jc));
    }
    setring St;
    ideal I = imap(G_t, I);
    Jc = imap(StA,Jc);
    list Zr = imap(StA,Zr);
    list FG = imap(StA, Lcom);
    // compute std using sparse multivariate rational interpolation
    ideal J = stdoverFF(I,pr,shft, Command, Zr, FG);
    setring G_t;
    ideal J = imap(St,J);
    if(finalTest(I,J))
    {
        attrib(J, "isSB", 1);
        option(set,opt);
        return(J);
    }
    else
    {
        print("FAIL");
        return(ffmodStd(I));
    }
}
example
{   "EXAMPLE:"; echo = 2;
    ring Ra=(0,a),(x,y,z),dp;
    ideal I = (a^2+2)*x^2*y+a*y*z^2, x*z^2+(a+1)*x^2-a*y^2;
    ffmodStd(I);
    ideal J = x^2*y+y*z^2, x*z^2+x^2-y^2;
    ffmodStd(J);
    ring R1=(0,a,b),(x,y,z),dp;
    ideal I = x^2*y^3*z+2*a*x*y*z^2+7*y^3,
              x^2*y^4*z+(a-7b)*x^2*y*z^2-x*y^2*z^2+2*x^2*y*z-12*x+by,
              (a2+b-2)*y^5*z+(a+5b)*x^2*y^2*z-b*x*y^3*z-x*y^3+y^4+2*a2*y^2*z,
              a*x^2*y^2*z-x*y^3*z+3a*x*y*z^3+(-a+4)*y^3*z^2+4*z^2-bx;
    ffmodStd(I);
}