// Singular-library
///////////////////////////////////////////////////////////////////////////////
// version string automatically expanded by CVS
version="$Id: ringgb.lib,v 1.16 2001/01/16 13:48:40 Singular Exp $";
category="Miscellaneous";
info="
LIBRARY:  ringgb.lib     Functions for coefficient rings
AUTHOR:  Oliver Wienand, email: wienand@mathematik.uni-kl.de

KEYWORDS: vanishing polynomial; zeroreduce; polynomial functions; library, ringgb.lib; ringgb.lib, functions for coefficient rings

PROCEDURES:
 findZeroPoly(f);        finds a vanishing polynomial for reducing f
 zeroReduce(f);          normal form of f concerning the ideal of vanishing polynomials
 testZero(poly f);       tests f defines the constant zero function
 noElements(def r);      the number of elements of the coefficient ring, if of type (integer, ...)
";

LIB "general.lib";
///////////////////////////////////////////////////////////////////////////////

proc findZeroPoly (poly f)
"USAGE:   findZeroPoly(f); f - a polynomial
RETURN:  zero polynomial with the same leading term as f if exists, otherwise 0
EXAMPLE: example findZeroPoly; shows an example
"
{
  list data = getZeroCoef(f);
  if (data[1] == 0)
  {
    return(0);
  }
  number q = leadcoef(f) / data[1];
  if (q == 0)
  {
    return(0);
  }
  poly g = getZeroPolyRaw(data[2]);
  g = leadmonom(f) / leadmonom(g) * g;
  return(q * data[1] * g);
  //return(system("findZeroPoly", f));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (integer, 65536), (y,x), dp;
  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
  findZeroPoly(f);
}

proc zeroReduce(poly f, list #)
"USAGE:   zeroReduce(f, [i = 0]); f - a polynomial, i - noise level (if != 0 prints all steps)
RETURN:  reduced normal form of f modulo zero polynomials
EXAMPLE: example zeroReduce; shows an example
"
{
   int i = 0;
   if (size(#) > 0)
   {
     i = #[1];
   }
   poly h = f;
   poly n = 0;
   poly g = findZeroPoly(h);
   if (i <> 0) {
     printf("reducing polyfct  : %s", h);
   }
   while ( h <> 0 ) {
      while ( g <> 0 ) {
         h = h - g;
         if (i <> 0) {
            printf(" reduce with: %s", g);
            printf(" to: %s", h);
         }
         g = findZeroPoly(h);
      }
      n = lead(h) + n;
      if (i <> 0) {
         printf("head irreducible  : %s", lead(h));
         printf("irreducible start : %s", n);
         printf("remains to check  : %s", h - lead(h));
      }
      h = h - lead(h);
      g = findZeroPoly(h);
   }
   return(n);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (integer, 65536), (y,x), dp;
  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
  zeroReduce(f);
  kill r;
  ring r = (integer, 2, 32), (x,y,z), dp;
  // Polynomial 1:
  poly p1 = 3795162112*x^3+587202566*x^2*y+2936012853*x*y*z+2281701376*x+548767119*y^3+16777216*y^2+268435456*y*z+1107296256*y+4244635648*z^3+4244635648*z^2+16777216*z;
  // Polynomial 2:
  poly p2 = 1647678464*x^3+587202566*x^2*y+2936012853*x*y*z+134217728*x+548767119*y^3+16777216*y^2+268435456*y*z+1107296256*y+2097152000*z^3+2097152000*z^2+16777216*z;
  zeroReduce(p1-p2);
}

proc testZero(poly f)
"USAGE:   testZero(f); f - a polynomial
RETURN:  returns 1 if f is zero as a function and otherwise a counterexample as a list [f(x_1, ..., x_n), x_1, ..., x_n]
EXAMPLE: example testZero; shows an example
"
{
  poly g;
  int j;
  bigint i = 0;
  bigint modul = noElements(basering);
  printf("Teste %s Belegungen ...", modul^nvars(basering));
  for (; i < modul^nvars(basering); i = i + 1)
  {
    if ((i + 1) % modul^(nvars(basering)/2) == 0)
    {
      printf("bisher: %s", i);
    }
    g = f;
    for (j = 1; j <= nvars(basering); j++)
    {
      g = subst(g, var(j), number((i / modul^(j-1)) % modul));
    }
    if (g != 0)
    {
      list counter = g;
      for (j = 1; j <= nvars(basering); j++)
      {
        counter = insert(counter, (i / modul^(j-1)) % modul);
      }
      return(counter);
    }
  }
  return(1);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (integer, 12), (y,x), dp;
  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
  zeroReduce(f);
  testZero(f);
  poly g = findZeroPoly(x2y3);
  g;
  testZero(g);
}

proc noElements(def r)
"USAGE:   noElements(r); r - a ring with a finite coefficient ring of type integer
RETURN:  returns the number of elements of the coefficient ring of r
EXAMPLE: example noElements; shows an example
"

{
  list l = ringlist(basering);
  return(l[1][2][1]^l[1][2][2]);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (integer, 233,6), (y,x), dp;
  noElements(r);
}

static proc getZeroCoef(poly f)
{
  if (f == 0)
  {
    return(0, leadexp(f))
  }
  list data = sort(leadexp(f));
  intvec exp = data[1];
  intvec index = data[2];
  intvec nec = 0:size(exp);
  int i = 1;
  int j = 2;
  bigint g;
  bigint G = 1;
  bigint modul = noElements(basering);
  bigint B = modul;

  for (; exp[i] < 2; i++) {if (i == size(exp)) break;}
  for (; i <= size(exp); i++)
  {
    g = gcd(B, G);
    G = G * g;
    B = B / g;
    if (g != 1)
    {
      nec[index[i]] = j - 1;
    }
    if (B == 1)
    {
      return(B, nec);
    }
    for (; j <= exp[i]; j++)
    {
      g = gcd(B, bigint(j));
      G = G * g;
      B = B / g;
      if (g != 1)
      {
        nec[index[i]] = j;
      }
      if (B == 1)
      {
        return(B, nec);
      }
    }
  }
  if (B == modul)
  {
    nec = 0;
    return(0, nec);
  }
  return(B, nec);
}

static proc getZeroPolyRaw(intvec fexp)
{
  list data = sort(fexp);
  intvec exp = data[1];
  intvec index = data[2];
  int j = 0;
  poly res = 1;
  poly tillnow = 1;
  int i = 1;
  for (; exp[i] < 2; i++) {if (i == size(exp)) break;}
  for (; i <= size(exp); i++)
  {
    for (; j < exp[i]; j++)
    {
      tillnow = tillnow * (var(1) - j);
    }
    res = res * subst(tillnow, var(1), var(index[i]));
  }
  return(res);
}

static proc getZeroPoly(poly f)
{
  list data = getZeroCoef(f);
  poly g = getZeroPolyRaw(data[2]);
  g = leadmonom(f) / leadmonom(g) * g;
  return(data[1] * g);
}

static proc findZeroPolyWrap (poly f)
"USAGE:   findZeroPolyWrap(f); f - a polynomial
RETURN:  zero polynomial with the same leading term as f if exists, otherwise 0
NOTE:    just a wrapper, work only in Z/2^n with n < int_machine_size - 1
EXAMPLE: example findZeroPoly; shows an example
"
{
  return(system("findZeroPoly", f));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (integer, 2, 16), (y,x), dp;
  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
  findZeroPoly(f);
}

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

/*
                           Examples:


// POLYNOMIAL EXAMPLES (Singular ready)
// ===================
//
// For each of the examples below, there are three equivalent polynomials. 'm' indicates the bit-widths of the
// input/output variables. For some of the polynomials, I have attached the RTL as well.
//
//
// 1) VOLTERRA MODELS:
//
//        A) CUBIC FILTER: (m = 32, 3 Vars)

LIB "ringgb.lib";
ring r = (integer, 2, 32), (x,y,z), dp;
poly p1 = 3795162112*x^3+587202566*x^2*y+2936012853*x*y*z+2281701376*x+548767119*y^3+16777216*y^2+268435456*y*z \
         +1107296256*y+4244635648*z^3+4244635648*z^2+16777216*z;
poly p2 = 1647678464*x^3+587202566*x^2*y+2936012853*x*y*z+134217728*x+548767119*y^3+16777216*y^2+268435456*y*z \
         +1107296256*y+2097152000*z^3+2097152000*z^2+16777216*z;
poly p3 = 1647678464*x^3+587202566*x^2*y+2936012853*x*y*z+134217728*x+548767119*y^3+16777216*y^2+268435456*y*z \
         +1107296256*y+2097152000*z^3+2097152000*z^2+16777216*z;
zeroReduce(p1-p2);
zeroReduce(p1-p3);
zeroReduce(p2-p3);

//        B) DEGREE-4 FILTER: (m=16 , 3 Vars)

LIB "ringgb.lib";
ring r = (integer, 2, 16), (x,y,z), dp;
poly p1 = 16384*x^4+y^4+57344*z^4+64767*x*y^3+16127*y^2*z^2+8965*x^3*z+19275*x^2*y*z+51903*x*y*z+32768*x^2*y  \
         +40960*z^2+32768*x*y^2+49152*x^2+4869*y;
poly p2 = 8965*x^3*z+19275*x^2*y*z+31999*x*y^3+51903*x*y*z+32768*x*y+y^4+32768*y^3+16127*y^2*z^2+32768*y^2 \
         +4869*y+57344*z^4+40960*z^2;
poly p3 = 8965*x^3*z+19275*x^2*y*z+31999*x*y^3+51903*x*y*z+32768*x*y+y^4+16127*y^2*z^2+4869*y+16384*z^3+16384*z;
zeroReduce(p1-p2);
zeroReduce(p1-p3);
zeroReduce(p2-p3);


// 2) Savitzsky Golay filter(m=16,5 Vars)

LIB "ringgb.lib";
ring r = (integer, 2, 16), (v,w,x,y,z), dp;
poly p1 = 25000*v^2*y+37322*v^2+22142*v*w*z+50356*w^3+58627*w^2+17797*w+17797*x^3+62500*x^2*z+41667*x \
         +22142*y^3+23870*y^2+59464*y+41667*z+58627;
poly p2 = 25000*v^2*y+4554*v^2+22142*v*w*z+32768*v+17588*w^3+25859*w^2+17797*w+17797*x^3+29732*x^2*z+32768*x^2 \
         +32768*x*z+8899*x+22142*y^3+23870*y^2+59464*y+41667*z+58627;
poly p3 = 25000*v^2*y+4554*v^2+22142*v*w*z+32768*v+17588*w^3+25859*w^2+17797*w+17797*x^3+29732*x^2*z+32768*x*z \
         +41667*x+22142*y^3+23870*y^2+59464*y+41667*z+58627;
zeroReduce(p1-p2);
zeroReduce(p1-p3);
zeroReduce(p2-p3);


// 3) Anti-alias filter:(m=16, 1 Var)

LIB "ringgb.lib";
ring r = (integer, 2, 16), c, dp;
poly p1 = 156*c^6+62724*c^5+17968*c^4+18661*c^3+43593*c^2+40224*c+13281;
poly p2 = 156*c^6+5380*c^5+1584*c^4+43237*c^3+27209*c^2+40224*c+13281;
poly p3 = 156*c^6+5380*c^5+1584*c^4+10469*c^3+27209*c^2+7456*c+13281;
zeroReduce(p1-p2);
zeroReduce(p1-p3);
zeroReduce(p2-p3);


// 4) PSK:(m=16, 2 Var)

LIB "ringgb.lib";
ring r = (integer, 2, 16), (x,y), dp;
poly p1 = 4166*x^4+16666*x^3*y+25000*x^2*y^2+15536*x^2+16666*x*y^4+31072*x*y+4166*y^4+15536*y^2+34464;
poly p2 = 4166*x^4+16666*x^3*y+8616*x^2*y^2+16384*x^2*y+15536*x^2+282*x*y^4+47456*x*y+53318*y^4+31920*y^2+34464;
poly p3 = 4166*x^4+16666*x^3*y+8616*x^2*y^2+16384*x^2*y+15536*x^2+282*x*y^4+47456*x*y+4166*y^4+15536*y^2+34464;
zeroReduce(p1-p2);
zeroReduce(p1-p3);
zeroReduce(p2-p3);

// Ref: A. Peymandoust G. De Micheli, “Application of Symbolic Computer Algebra in High-Level Data-Flow
// Synthesis,” IEEE Transactions on CAD/ICAS, Vol. 22, No. 9, September 2003, pp.1154-1165.

*/