generate_equations.singular

LIB "na.lib.singular";


ring r=0,(aa,ab,ac,ad,ba,bb,bc,bd,ca,cb,cc,cd,da,db,dc,dd),dp;	

// The non-associative variables we will use.
list na_vars = list("a","b","c","d");

// Initialize the non-associative polynomial library
// for use with the variables 'na_vars'.
na_init(na_vars);


/************************
*** UTILITY FUNCTIONS ***
************************/


// Does the list have any repeated elements?

proc list_elements_are_unique(list as)
{
  int j;
  for (int i=1; i<=size(as); i++)
  {
    for (j=i+1; j<=size(as); j++)
    {
      if( as[i] == as[j] )
      {
        return( 0 );
      } 
    }
  }
  return( 1 );
}


proc set_equality(list as, list bs)
{
  assert( list_elements_are_unique(as),
          "'set_equality' assumes the arguments don't have repeated elements");
  assert( list_elements_are_unique(bs),
          "'set_equality' assumes the arguments don't have repeated elements");

  if( size(as) != size(bs) )
  {
    return( 0 );
  }
  int count_in_bs;
  int j;
  for (int i=1; i<=size(as); i++)
  {
    count_in_bs = 0;

    for (j=1; j<=size(bs); j++)
    {
      if( as[i] == bs[j] )
      {
        count_in_bs++;
      }
    }  
    
    if( count_in_bs != 1 )
    {
      return( 0 ); 
    }      
  }
  return( 1 );
}


/************************************
*** THE POLYNOMIAL TRANSFORMATION ***
************************************/


// Maps "L" to "R" and "R" to "L".

proc swapLR(string side)
{
  if(side == "L")
  {
    return( "R" );
  }
  else
  {
    if(side == "R")
    {
      return( "L" );
    }
    else
    {
      ERROR("Invalid input to 'swapLR'");
    }
  }
}


// Coefficient creation for 'mksubst'.

proc mkcoef(string suffix, int i)
{
  // We expect a suffix LL, LR, RR or RL
  // The corresponding indices are LL -> 1, LR -> 2, RR -> 3, RL -> 4
  int suffix_index = find( "LLRRL", suffix ); 

  // 1 -> a, 2 -> b, 3 -> d, 4 -> c (note the order of d and c)
  list suffix_translations = list("a","b","d","c");

  list i_translations = list("a","b","c","d");
  
  string coef_name = suffix_translations[suffix_index] + i_translations[i];

  execute("na_poly p = " + coef_name);

  return( p );
}


// This function creates the RHS (Right Hand Side) of the identities
// that are derived from the associativity of the co-smash.
// The LHS of the identity is described by the arguments of this
// function. The LHS is a product of the monomials 'x', 'y' and 'z',
// with the placement of the parentheses given by 'suffix'.
// 'x' is the monomial whose location in the multiplication tree is
//     given by 'suffix'. It is the monomial "we are pulling one step
//     out of the parentheses".
// 'y' is inside the parentheses along with 'x'.
// 'z' is outside of the parentheses.
// For example, the 'suffix' "LR" corresponds to the LHS '(y*x)*z'.

proc mksubst(na_mono x, na_mono y, na_mono z, string suffix)
{
  na_poly subst_rhs = mkcoef(suffix,1) * x*(y*z)
                    + mkcoef(suffix,2) * x*(z*y)
                    + mkcoef(suffix,3) * (y*z)*x 
                    + mkcoef(suffix,4) * (z*y)*x ;

  return( subst_rhs );
}


// This function maps the non-associative monomial
// 'm' to a non-associative polynomial.
// We do this by "pulling the submonomial 'subm'"
// one step out of the parentheses" using the
// identities generated by 'mksubst'.

proc transform_monomial(na_mono subm, na_mono m)
{
  // Depending on the path to the submonomial, we now make a
  // substitution based on what the end of the path looks like. 
 
  // list patterns = list( "LL", "LR", "RL", "RR" );
  // These correspond to (xu)v, (ux)v, u(xv), u(vx),
  // where x denotes the monomial we are looking for
  // and u, v are submonomials of degree >= 1.
  // Find the monomial for patterns LL, LR, RL, RR
  // Card (path(LL) + ... path(RR)) == 1
  // index = 1 * path(LL) + ... + 4* path(RR);

  assert( na_is_product(m) ,
    "'monomial_substitute_rec' first argument not of degree >=2" );
  
  list occurrences = na_find_submonomial(m, subm);

  // We assume that each non-associative variable appears in the
  // monomial 'm' at most once.
  assert(size(occurrences) <= 1, "Too many occurrences of the submonomial.");

  if( size(occurrences) == 0 )
  {
    // If 'subm' does not appear in 'm' as a submonomial,
    // we return 'm' as a 'na_poly'.
    return( na_poly(m) );
  }
  else
  {
    string path = occurrences[1];

    if(size(path) < 2)
    {
      return( na_poly(m) );
    }


    // We take the last two characters of the path.
    string suffix = path[size(path)-1,2];

    string reduced_path = path[1, size(path)-2]; // 'path' minus 'suffix'

    // The substitution we want to make depends on the 'suffix'.

    na_mono x = subm;
    na_mono y = na_get_submonomial(m, reduced_path + suffix[1] + swapLR(suffix[2]));
    na_mono z = na_get_submonomial(m, reduced_path + swapLR(suffix[1]));

    // 'x', 'y' and 'z' are submonomials of 'm', which combine to form
    // a submonomial such as '(x*y)*z' or 'z*(x*y)'. The order of the
    // factors varies (depending on 'suffix'), but we always have that
    // 'x' is the submonomial we want to "pull one step out of the parentheses",
    // 'y' is also inside the parentheses with 'x' and
    // 'z' is outside of the parentheses.
    
    return( na_subst_na_poly_into_submonomial(m,
                                             mksubst(x,y,z,suffix),
                                             reduced_path) );
  }
}


// Apply 'transform_monomial' to all monomials of the
// non-associative polynomial 'p'.

proc H(na_mono subm, na_poly p)
{
  na_poly new_p = 0;

  for(int i=1; i<=size(p.terms); i++)
  { 
    new_p = new_p + na_get_coef(p.terms[i]) * transform_monomial(subm, na_get_mon(p.terms[i]));
  }
  return( new_p );
}




/*******************************
*** GENERATING THE EQUATIONS ***
*******************************/


// The pattern for generating our system of equations:

// 1) Define a polynomial.
// na_poly p = (a*b)*c;

// 2) Transform the polynomial in two different ways which would
//    both be the identity transformation in the variety V.
// na_poly lhs = H(b,p);
// na_poly rhs = H(b,H(a,p));

// 3) Take the difference and add the coefficients to the list of equations.
// eqs = eqs + na_get_coefficients( rhs - lhs );




proc generate_deg2_equations()
{
  list eqs_deg2 = list();

  na_poly p = (a*b)*c;
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(c,H(a,p)) - p );
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(c,H(b,p)) - p );
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(b,H(a,p)) - H(b,p) );
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(a,H(b,p)) - H(a,p) );

  p = a*(b*c);
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(a,H(b,p)) - p );
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(a,H(c,p)) - p );
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(b,H(c,p)) - H(b,p) );
  eqs_deg2 = eqs_deg2 + na_get_coefficients ( H(c,H(b,p)) - H(c,p) );

  // Each of the assignments above adds 4 equations to the list
  // 'eqs_deg2', so in total 'eqs_deg2' contains 8*4 = 32 equations.
  // assert( size(eqs_deg2) == 32 );

  return( eqs_deg2 );
}


proc generate_deg3_equations()
{
  list eqs_deg3 = list();

  na_poly p = a*(b*(c*d));
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(c*a,H(a*c,H(d,H(b,p))))
                                            - H(d*b,H(b*d,H(c,p))) );
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(d*a,H(a*d,H(c,H(b,p))))
                                            - H(c*b,H(b*c,H(d,p))) );

  p = a*((b*c)*d);
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(c*a,H(a*c,H(b,H(d,p))))
                                            - H(d*b,H(b*d,H(c,p))) );
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(b*a,H(a*b,H(c,H(d,p))))
                                            - H(d*c,H(c*d,H(b,p))) );

  p = (a*(b*c))*d;
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(c*d,H(d*c,H(b,H(a,p))))
                                            - H(a*b,H(b*a,H(c,p))) );
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(b*d,H(d*b,H(c,H(a,p))))
                                            - H(a*c,H(c*a,H(b,p))) );

  p = ((a*b)*c)*d;
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(b*d,H(d*b,H(a,H(c,p))))
                                            - H(a*c,H(c*a,H(b,p))) );
  eqs_deg3 = eqs_deg3 + na_get_coefficients ( H(a*d,H(d*a,H(b,H(c,p))))
                                            - H(b*c,H(c*b,H(a,p))) );

  // Each of the assignments above adds 8 equations to the list
  // 'eqs_deg3', so in total 'eqs_deg3' contains 8*8 = 64 equations.
  // assert( size(eqs_deg3) == 64 );

  return( eqs_deg3 );
}



proc write_equations(list eqs)
{
  string eqs_string = "";

  for(int i=1; i<=size(eqs); i++)
  {
    eqs_string = eqs_string + string(eqs[i]);
    if( i != size(eqs) )
    {
      eqs_string = eqs_string + "," + newline + newline; 
    }
  } 
  write(":w equations.txt",eqs_string);
}


proc read_equations()
{ 
  string eqs_string = read("equations.txt");

  list eqs;

  execute("eqs = list(" + eqs_string + ")");

  return( eqs );
}


/****************
*** CALCULATE ***
****************/

//////////////////////////////
int start_time = rtimer;
print("GENERATING EQUATIONS");
//////////////////////////////

list eqs = generate_deg2_equations() + generate_deg3_equations();

/////////////
print("DONE (time taken: " + string((rtimer - start_time) div system("--ticks-per-sec")) + " seconds)");
/////////////

print("WRITING EQUATIONS TO DISK");

write_equations(eqs);

// We try to read back what we wrote and see if what
// we read is what we attempted to write.
if( (status(":r equations.txt","exists") == "yes") and
     set_equality(eqs,read_equations() ) )
{
  print("SUCCESS!"); 
}
else
{
  print("FAILED");
}

quit;