/*  Group Operations   */

#include "global.h"
#include <stdio.h>
#include <stdlib.h>

/* Imports from polynomial.c.  */

void initexponent(exponent exp);
void initpolynomial(polynomial f);
void advance(exponent exp);
void polyprod(polynomial f, polynomial g, polynomial h);
void polypower(polynomial f, int n, polynomial g);
void polysum(polynomial f, polynomial g);
void writepoly(polynomial f);
void polymult(polynomial f, polynomial g); 
int readpoly(polynomial f);
int polyequal(polynomial f, polynomial g);

 /* Prototypes  */

void initgroupelt(group_elt g);
void gprod(group_elt f, group_elt g, group_elt h);
void gcopy(group_elt f, group_elt g);
void gmult(group_elt f, group_elt g);
void gpower(group_elt g, int n);
void writegpelt(group_elt g);
int readgpelt(group_elt g);
int checkgroupelt(group_elt g);
int groupequal(group_elt f, group_elt g);
int groupnotid(group_elt f);
void groupinv(group_elt f, group_elt g);
void commutator(group_elt f, group_elt g, group_elt h);
void conjugate(group_elt f, group_elt g, group_elt h);

extern int rank, total, prime, help;

void initgroupelt(group_elt g)
{
  int i;
  for (i=0; i < rank; i++)
    initpolynomial(g[i]);
}

void gprod(group_elt f,group_elt g, group_elt h)
{
  int i, j, k, a;
  exponent expf;
  polynomial term, t;
  initgroupelt(h);
  for(j = 0; j < rank; j++)
    for(i = 0, initexponent(expf); i < total; i++, advance(expf))
    if (a = f[j][i])
    {
       initpolynomial(term);
       term[0] = a;
       for(k = 0; k < rank; k++)
          if (expf[k])
          {
            polypower(g[k], expf[k], t);
            polymult(term, t);
          }
       polysum(h[j], term);
    }
}    

void gcopy(group_elt g, group_elt h)  /* g = h */
{
  int i, j;
  for (i = 0; i < rank; i++)
    for (j = 0; j < total; j++)
      g[i][j] = h[i][j];
}
 
void gmult(group_elt f,group_elt g)
{
  group_elt t;
  gprod(f, g, t);
  gcopy(f, t);
}

void gpower(group_elt g, int n)
{
  group_elt h;
  if (n <= 0) {printf("Bad value of n in gpower %d\n",n); exit(0);}
  if (n == 1) return;
  if (n%2)
  {
    gcopy(h, g);
    gpower(g, n-1);
    gmult(g, h);
  }
  else 
  {
    gmult(g, g);
    gpower(g, n/2);
  }
}

void writegpelt(group_elt g)
{
  int i;
  for(i=0; i < rank; i++)
  {
    writepoly(g[i]);
    printf("\n");
  }
}

int readgpelt(group_elt g)
{
  int i;
  if (help) 
    printf("Enter group element as a sequence of %d polynomials.\n", rank);
  for (i = 0; i < rank; i++)
    if ( !readpoly(g[i]) ) return 0;
  if (!checkgroupelt(g))
  {
    printf("Not an element of the Nottingham group.\n");
    return 0;
  }
  return 1;
}

int checkgroupelt(group_elt g)
{
  int i, j;
  for (i = 0; i < rank; i++)
  for (j = 0; j <= rank; j++)
    if ( g[i][j] != (i + j == rank? 1:0) ) return 0;
  return 1;
}

int groupequal(group_elt f,group_elt g)
{
  int i;
  for(i = 0; i < rank; i++)
    if ( !polyequal(f[i], g[i]) ) return 0;
  return 1;
}

int groupnotid(group_elt f)
{
  int i, j;
  for(i = 0; i < total; i++)
    for(j = 0; j < rank; j++)
      if ( f[j][i] != (i + j == rank?  1:0) ) return i + 1;
  return 0;
}

void groupinv(group_elt f,group_elt g)
{
  group_elt s, t;
  int i, j, k;
  for(i = 0; i < rank; i++)
    for(j = 0; j < total; j++)
      g[i][j] = ( j <= rank? f[i][j]:(prime - f[i][j])%prime );
  if ( k = groupnotid(f))
  {
    gprod(f, g, t);
    groupinv(t, s);
    gmult(g, s);
  }
}

void commutator(group_elt f, group_elt g, group_elt h)
{
  group_elt k;
  groupinv(f, h);
  groupinv(g, k);
  gmult(h, k);
  gmult(h, f);
  gmult(h, g);
}
  
 void conjugate(group_elt f, group_elt g, group_elt h)
{
  groupinv(f, h);
  gmult(h, g);
  gmult(h, f);
}