\\\\\\\\\\\\\\\\\\\\\\\
\\ Stampa vettore in verticale

printVector(v) =
{
  local(i);
  
  for(i = 1,matsize(v)[2],
    print(v[i]);
  );
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ Vettore di "punti" rappresentanti i "vertici"
\\ (coefficienti dei polinomi fattori)

vertici(n) =
{
  local(i,v,res);

  v   = vector(n,i,vector(2,i,i-1));
  res = vector(2^n, i, vector(n));
  i = 0;
  forvec(x = v, res[i++] = x);
  
  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ Conversione vertice -> coefficiente

coefficiente(vertice, letter = "a") =
{
  local(i,res = Str(letter));
  for(i=1,matsize(vertice)[2], res = Str(res,vertice[i]));

  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ Conversion vertex -> monomial

termine(vertex, letter = "a", var = "x") =
{
  local(i,n,res = Str(letter));
  n = matsize(vertex)[2];
  for(i=1,n, res = Str(res,vertex[i]));
  if ( vertex != 0,
   res = Str(res, "*");
   for(i=1,n, if(vertex[i],res = Str(res,var,"_",n-i)));
  );
  
  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Return a vector with the coefficients of a generic n-variate
\\ square polynomial (with default "a" letter) with maxdegree 1.

coeffSquarePoly(n,a = "a", getIndexes = 0) =
{
  local (tmp = vector(2^n,i,
    concat(vector(n - matsize(binary(i-1))[2],i,0), binary(i-1))) );

  if ( getIndexes == 0,
    for (i=1,matsize(tmp)[2],
      str = "";
      for (j = 1,n,
        str = Str(str,tmp[i][j]);
      );
      tmp[i] = eval(Str(Str(a),str));
    );
  );
  
 return (tmp);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ N-variate terms with degree k;

terminiDiGradoKinNvariabili(n,k,x = "x") = 
{
 local(res);

 if ( k == 0,
   res = [1],
  if ( n == 1,
    res = [eval(Str(x,"_1"))^k];
  , \\ else
    res = concat(terminiDiGradoKinNvariabili(n-1,k,x),
                 eval(Str(x,"_",n))*terminiDiGradoKinNvariabili(n,k-1,x));
  ));

 return(res);
}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Return a vector with the terms of a generic multivariate
\\ square polynomial with maxdegree 1 (with letter "x").

t(n) =
{
  local(tmp = c(n,"x",1), res = 1 );

  for (i=1,matsize(tmp)[2],
    res = 1;
    for (j = 1,n,
      if ( tmp[i][j] == 1, res *= eval(Str("x_",n-j));
    ));
    tmp[i] = eval(res);
  );

  return (tmp);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Return a generic n-variate square polynomial with degree 1
\\ "x" with coefficients indicated with letter given by a.

polySquare(n,a = "a") =
{
  local(co = c(n,Str(a)), x = t(n), tmp = 0);

  for(i=1,matsize(co)[2],
   tmp += x[i]*co[i];
  );
  
  return (tmp);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Return a vector with the coefficients of a generic n-variate
\\ square polynomial (with default "a" letter) with maxdegree 1.

c(n,a = "a", getIndexes = 0) =
{
  local (tmp = vector(2^n,i,
    concat(vector(n - matsize(binary(i-1))[2],i,0), binary(i-1))) );

  if ( getIndexes == 0,
    for (i=1,matsize(tmp)[2],
      str = "";
      for (j = 1,n,
        str = Str(str,tmp[i][j]);
      );
      tmp[i] = eval(Str(Str(a),str));
    );
  );
  
 return (tmp);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Prints in a standard format a generic square polynomial of degree 1
\\ in n "x" variables with coefficients with letter a.

printSquarePoly(n,a="a") =
{
  local(co = c(n,Str(a)), x = t(n),
        str = Str(co[1],if(x[1]!= 1, Str("*",x[1]))));

  for(i=2,matsize(co)[2],
   str = Str(str," + ",co[i],if(x[i]!= 1, Str("*",x[i])));
  );

  print(str);
  return(str);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Conversion from integer to base 3 representation.

daInteroAbaseTre(n, vectorLength) =
{
  local(res = vector(vectorLength), i = n, j = 0, tmp = 0);

  while( i,
   tmp = divrem(i,3);
   res[j++] = tmp[2];
   i = tmp[1];
  );

  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Conversion from base 3 representation to integer.

daBaseTreAintero(v) =
{
  local(vectorLength = matsize(v)[2], i, res = 0);

  forstep(i = vectorLength, 1, -1,
    res *= 3;
    res += v[i];
  );

  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Conversion from integer to base 2 representation.

daInteroAbaseDue(n, vectorLength) =
{
  local(res = vector(vectorLength), j = vectorLength + 1, tmp = 1);

  while( j > 1,
   res[j--] = if(bitand(n,tmp), 1 , 0);
   tmp <<= 1;
  );

  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Conversion from base 2 representation to integer.

daBaseDueAintero(v) =
{
  local(vectorLength = matsize(v)[2], i, res = 0);

  forstep(i = vectorLength, 1, -1,
    res *= 2;
    res += v[i];
  );

  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Vectors with length n with exactly k "1" (we suppose n >= k)
\\ and all other entries equal to 0 or 2.

vettoriLunghiNconKuni(n,k) = 
{
 local(res, tmp1,tmp2, dim1, dim2, i, j = 0, h);

 if ( k == 0,   res = matrix(2^n,n);
                forvec(x=vector(n,i,[0,1]), res[j++,] = vector(n,h,2*x[h]));
 , \\ else
   if ( n == k, res = matrix(1,n,i,j,1);
   , \\ else
    tmp1 = vettoriLunghiNconKuni(n-1,k);   \\ Primo elemento: 0
    tmp2 = vettoriLunghiNconKuni(n-1,k-1); \\ Primo elemento: 1
    dim1 = matsize(tmp1)[1];
    dim2 = matsize(tmp2)[1];
    res = matrix(2*dim1 + dim2, n);
    for(i=1,dim2,  for(j=1,n-1, res[i,j] = tmp2[i,j];
                   res[i,n] = 1;)
    );
    for(i=1,dim1,  for(j=1,n-1, res[i+dim2,j] = tmp1[i,j];
                   res[i+dim2,n] = 0;);
    );
    dim2 += dim1;
    for(i=1,dim1,  for(j=1,n-1, res[i+dim2,j] = tmp1[i,j];
                   res[i+dim2,n] = 2;);
    );
  ));
  return(res);
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ Position of the i-th 1 on the right in the argument vector
\\ (starting from 0).

posizioneDellIesimoUnopiuAdestra(v, i = 1) =
{
  local(res = matsize(v)[2], position = i+1);

  while( position--,
    while( (res >= 1 ) && (v[res] != 1), res--; );
    res--;
  );

  return(res);
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ Posizione dell'1 pił a destra nel vettore argomento (partendo da 0).

posizioneDelPrimoUno(v) =
{
 return (posizioneDellIesimoUnopiuAdestra(v));
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ Count the number of 1 (in base 3) in indice (dimension of the
\\ vertex itself, not to consider the case will all 0, giving 0).

peso(indice) =
{
  local(res=0);

  return(sum(i=1,matsize(indice)[2],if(indice[i]==1,1,0),res));
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ Ordering depending on the number of 1.

ordinaVettorePerPeso(v) =
{
  local(size = matsize(v)[2], v1, v2, halfSize = size >> 1,
        i = 1, j = 1, h = 1, res = vector(size));

  if (size == 1, return (v); ); \\ Basic case
  v1 = ordinaVettorePerPeso(vecextract(v,Str("1..",eval(halfSize))));
  v2 = ordinaVettorePerPeso(vecextract(v,Str(eval(halfSize+1),"..",eval(size))));
  while( i <= halfSize && j <= size - halfSize,
    if( peso(v1[i]) < peso(v2[j]),
      res[h] = v1[i]; i++;
    ,
      if( peso(v1[i]) > peso(v2[j]),
        res[h] = v2[j]; j++;
      ,
        if( posizioneDelPrimoUno(v1[i]) < posizioneDelPrimoUno(v2[j]),
          res[h] = v1[i]; i++;
        ,
          res[h] = v2[j]; j++;
        );
      );
    );
    h++;
  );

  if( i > halfSize,  \\ Visitato tutto il primo vettore: lavora sul secondo.
    while( h <= size,
      res[h] = v2[j]; j++; h++;
  ));
  if( j > size - halfSize,  \\ Visitato il secondo vettore: lavora sul primo.
    while( h <= size,
      res[h] = v1[i]; i++; h++;
  ));

  return (res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ All indexes of product elements represented as 
\\ vectors in base 3. The vertexes are not indicated,
\\ but there are sides (just one "1"), faces (two "1"), etc.
\\ depending on the number of "1" in base 3 expansion.

indiciProdotto(n) =
{
\\  local(res = vector(3^n-2^n, i, daInteroAbaseTre(i-1,n)));
\\
\\  return(ordinaVettorePerPeso(res));

  local(res = matrix(3^n, n), i, j, k = 0, tmp);

  for(i=1,n,
    tmp = vettoriLunghiNconKuni(n,i);
    for(j=1,matsize(tmp)[1],
      res[k++,] = tmp[j,];
    );
  );

  return(res);
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ VALUTATION: currently the argument is n, but one should have the
\\ coefficient vectors of the two factors. In position i there is the
\\ coefficient whose multiindex correpond to the expansion
\\ in base 2 of (i-1) -- in gp indexes start from 1, not from 0.
\\
\\ The total number (optimal) of additions is 2*(3^n - 2^n)
\\ (the 2 derives from the fact that there are two factors, the -2^n
\\ from the fact that on vertexes no addition is needed).

valutazione(n) =
{
  local(resA   = vector(3^n), resB = vector(3^n), res = vector(3^n),
        indici = indiciProdotto(n),i, index, coeffA, coeffB, add = 0);

  coeffA = coeffSquarePoly(n,"a"); \\ Coefficienti dei fattori:
  coeffB = coeffSquarePoly(n,"b"); \\ questi dovrebbero essere l'input.

  \\ First vertexes (no computation, just a copy).

  for(i = 1, 2^n, index = daBaseTreAintero(daInteroAbaseDue(i-1,n))<<1 + 1;
                  resA[index] = coeffA[i]; resB[index] = coeffB[i]; );

  \\ Then all the rest, from sides upwards.
  
  for(i = 1, 3^n-2^n,
    index  = daBaseTreAintero(indici[i,]) + 1;    \\ Index
    offset = 3^posizioneDelPrimoUno(indici[i,]);  \\ Offset
    resA[index] = resA[index-offset] + resA[index+offset]; add++;
    resB[index] = resB[index-offset] + resB[index+offset]; add++;
  );

  print("-------------\nresA ="); printVector(resA); \\ Print the two factors
\\  print("-------------\nresB ="); printVector(resB);
  
  for(i = 1, 3^n, res[i] = resA[i]*resB[i]);  \\ Compute the product

  print("Additions   = ",add,"   :   2*(3^",n," - 2^",n,") = ",(3^n-2^n)<<1);
  return(res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Initialize the matrix with indexes for sections beginnings
\\ and dimension.

inizializzazioneInizioEdimensioneSezioni(n) =
{
 local(res = vector(n,i,vector(2)), k);

 res[1] = [1,n<<(n-1)];
 for(k = 2, n,
  res[k] = [res[k-1][1] + res[k-1][2], binomial(n,k)<<(n-k)]; );

 return(res);
}

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Subtract from all faces the two faces with 0 and 2
\\ corresponding to the "1" in position posizioneUno starting from left.

sottrazioneFacce(v,indici, base, limite, posizioneUno) =
{
  local(i, indice, offset);

  for(i = base, base + limite - 1, 
    indice = daBaseTreAintero(indici[i,]) + 1;
    offset = 3^(posizioneDellIesimoUnopiuAdestra(indici[i,], posizioneUno));

    v[indice] -= (v[indice-offset] + v[indice+offset]);
  );

  return(v);
}

\\\\\\\\\\\\\\\\\\\\\\\
\\ INTERPOLATION: interpolates the values in vector v, having 3^n entries,
\\ obtaining the final result. In total 2n*3^(n-1) subtractions.

interpolazione(v) =
{
  local(n = ceil(log(matsize(v)[2])/log(3)),i,j, sub = 0,
        indiciOrdinati, inizioEdimensioneSezioni);
  print("N. of variables = ",n,"\n"); \\ Print variables number.

  indiciOrdinati = indiciProdotto(n);

  \\ Vettore con gli indici indicanti l'inizio delle sezioni
  \\ corrspondenti agli insiemi delle varie facce.
  inizioEdimensioneSezioni = inizializzazioneInizioEdimensioneSezioni(n);

  forstep(i = n, 1, -1,            \\ Si trattano le varie "facce" con i "1".

   \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
   \\ Da ciascuna faccia con dimensione maggiore si
   \\ sottraggono le due corrispondenti facce aventi 0 e 2 al posto
   \\ del j-esimo 1 pił a destra dove j dipende dalla differenza
   \\ delle dimensioni.

   for( j = i, n,   \\ print("j = ",j);
     v = sottrazioneFacce(v, indiciOrdinati, inizioEdimensioneSezioni[j][1],
         		  inizioEdimensioneSezioni[j][2], j-i+1 );
     sub += inizioEdimensioneSezioni[j][2]<<1;
  ));

  print("Subtractions = ",sub,"   :   2n*3^(n-1) = ",(n*3^(n-1))<<1);
  return(v);
}

\\****************************************************\\
\\****************************************************\\
\\ CODICE (sembra funzioni !)

n = 4;
v = valutazione(n);
v = interpolazione(v);
\\printVector(v);