//--- Alcune procedure che necessitano general.lib.
LIB "general.lib";
// ** loaded /usr/bin/../share/singular/LIB/general.lib (4.0.0.1,Jan_2014)
// ** loaded /usr/bin/../share/singular/LIB/matrix.lib (4.0.0.0,Jun_2013)
// ** loaded /usr/bin/../share/singular/LIB/nctools.lib (4.0.0.0,Jun_2013)
// ** loaded /usr/bin/../share/singular/LIB/random.lib (4.0.0.0,Jun_2013)
// ** loaded /usr/bin/../share/singular/LIB/ring.lib (4.0.2.2,Jan_2016)
// ** loaded /usr/bin/../share/singular/LIB/primdec.lib (4.0.2.0,Apr_2015)
// ** loaded /usr/bin/../share/singular/LIB/absfact.lib (4.0.0.0,Jun_2013)
// ** loaded /usr/bin/../share/singular/LIB/triang.lib (4.0.0.0,Jun_2013)
// ** loaded /usr/bin/../share/singular/LIB/elim.lib (4.0.0.1,Jan_2014)
// ** loaded /usr/bin/../share/singular/LIB/poly.lib (4.0.0.0,Jun_2013)
// ** loaded /usr/bin/../share/singular/LIB/inout.lib (4.0.0.0,Jun_2013)

proc generatore(int d)
//--- polinomio generatore
{
  list L;
  int i;
  for (i=1;i<=d-1;i++)
    {
     L[i]=var(1)-par(1)^i;
    }; 
  poly g = product(L);
  return(g);
};

proc sindrome(int d, poly y)
//--- polinomio sindrome
{
  list L;
  int i;
  for (i=1;i<=d-1;i++)
    {
     L[i]=subst(y,var(1),par(1)^i)*x^(i-1);
    }; 
  poly s = sum(L);
  return(s);
};

proc codifica (int d, poly w)
{
  return(w - reduce(w,generatore(d)));
};

proc decodifica (int d, poly y)
{
  poly x = var(1);
  number a = par(1);
  int n = 1; while (a^n != 1) {n++;};
  vector v1 = [x^(d-1),0];
  vector v2 = [sindrome(d,y),1];
  module K = v1, v2;
//--- modulo delle soluzioni dell'equazione chiave  
  option(redSB);
  module G = std(K); 
//--- G = base ridotta di Groebner per il modulo K
  if (G[1]<G[2]) 
    {
     poly omega = G[1][1];
     poly lambda = G[1][2]; 
//--- lambda = polinomio locatore degli errori
//    a meno di un fattore costante
    } 
  else
    {
     poly omega = G[2][1];
     poly lambda = G[2][2];
    };
  int i;
  list I;
//---I = insieme delle posizioni degli errori  
  for (i=1; i<=n; i++) 
    {
     if (subst(lambda,x,a^i) == 0) {I=insert(I,n-i);}
    };
  print(I);
  poly e; 
//--- calcolo del polinomio errore e con la formula di Forney
  for (i=1; i<=size(I); i++)
    {
     poly chi(i) = subst(lambda/(1-a^I[i]*x),x,a^(-I[i]));
     e = e + subst(omega,x,a^(-I[i]))/(a^I[i]*chi(i))*x^I[i];
    };
  return(y - e);
};


//--- Esempio:
//---[CLO, Using ..., p. 445], Reed-Solomon [8,4,8-4+1]
//--- si possono correggere 2 errori
//--- c = x^7 + 2*x^5 + x^2 + 2*x +1
//--- y = x^7 + a*x^5 + (a-1)*x^2 + 2*x + 1 
//--- posizioni di informazione: x^7,x^6,x^5,x^4
//--- posizioni di controllo: x^3,x^2,x^1,1
ring r = (3,a),x,(lp,c); minpoly = a^2+a+2;
r;
//   characteristic : 3
//   1 parameter    : a 
//   minpoly        : (a2+a-1)
//   number of vars : 1
//        block   1 : ordering lp
//                  : names    x
//        block   2 : ordering c
> poly w = x^7 + 2*x^5;
> poly c = codifica(5,w);
> c;
x7-x5+x2-x+1
> //---decodifica di y (i coefficienti di x^5 e x^2 sono errati)
. poly y = x^7 + a*x^5 + (a+2)*x^2 + 2*x + 1;
> decodifica(5,y);
[1]:
   2
[2]:
   5
x7-x5+x2-x+1
> poly y = x^7 + a*x^5 + (a+2)*x^2 + 2*x + a;
// ** redefining y **
> decodifica(5,y);
[1]:
   1
[2]:
   6
x7+(-a)*x6+(a)*x5+(a-1)*x2+x+(a)
> poly y = x^7 + a*x^5 + (a+2)*x^2 + a*x + 1;
// ** redefining y **
> decodifica(5,y);
empty list
x7+(a)*x5+(a-1)*x2+(a)*x+1
> sindrome(5,y);
(-a-1)*x3+(-a)*x2+(a-1)*x+(a)
monitor("");