/*  Polynomial Code.   
The type 'exponent' is an array of int length 'rank', which is the number
of variables.  An object exp of this type defines a monic monomial, in which
the exponent of x_i is exp[i].
The function degree(exp) returns i_0+...+i_{rank - 1}.

Exports to group_elt .c
  rank, prime, total:  integer variables;
  initexponent, initpolynomial:  set to zero;
  advance: increments an exponent;
  polyprod, polypower, polysum: return the appropriate polynomials;
  writepoly.
*/

#include "global.h"

#include <stdio.h>

#include <stdlib.h>

/* Import from utilities .c */
int choose(int n, int k);

/* Prototypes */

void initialise(void);
int degree(exponent exp);
void printexp(exponent exp);
int readexp(exponent exp);
int extra(exponent exp, int k);
int card_monos(int deg, int k);
int card_monos_atmost(int deg, int k);
int position(exponent exp);
int readfield(void);
void initpolynomial(polynomial f);  /* Export to group_elt.c */
void initexponent(exponent exp);  /* Export to group_elt.c */
int readpoly(polynomial f);  /* Export to group_elt.c  */
void advance(exponent exp);  /* Export to group_elt.c */
void writepoly(polynomial f);  /* Export to group_elt.c */
void expsum(exponent exp1, exponent exp2, exponent exp3);
void polyprod(polynomial f, polynomial g, polynomial h);  /* Export to group_elt.c */
void polymult(polynomial f, polynomial g);  /* Export to group_elt.c */
void polypower(polynomial f, int n, polynomial h);  /* Export to group_elt.c */
void polysum(polynomial f, polynomial g);  /* Export to group_elt.c */
int polyequal(polynomial f, polynomial g);  /* Export to group_elt.c */

int rank;    /* The number of variables */
int prime;   /* Work modulo prime */
int mdegree; /* The maximum degree of a polynomial */
int total;   /* The total number of monomials in a poly */
int help;    /* Does the user need prompting? */

void initialise(void)  /* Sets the global variables  */
{ 
  printf("Enter the number of variables \n");
  scanf("%d",&rank);
  printf("Enter the prime \n");
  scanf("%d",&prime);
  printf("Enter the maximum degree of a polynomial \n");
  scanf("%d",&mdegree);
  total = card_monos_atmost(mdegree,rank);
  if (rank > MAXR || total > MAXL)
    { printf("Too big, reset constants in global.h\n"); 
      printf(" rank is %d, MAXR is %d, total is %d, MAXL is %d.\n", rank, MAXR, total, MAXL);
      exit(0);
    }
  help = 2;  /* Maximum help  */
}

int degree(exponent exp)
{
  int ans, i;
  for (ans = i = 0; i<rank; i++)
    ans += exp[i];
  return ans;
}

void printexp(exponent exp)
{ 
  int i;
  for (i=0; i < rank; i++)
    if (exp[i]) { 
      printf("x");
      if (rank > 1) printf("_%d", i);
      if (exp[i] != 1) printf("^%d",exp[i]);
    }
}

int readexp(exponent exp)
{ 
  int i, h;
  h = help;
  if (help)
    printf("Enter an exponent\n");
  if ((help > 1)  && (rank > 1))
     { 
       printf("eg 0 1 4 for x_0^0 x_1^1 x_2^4\n");
       printf("you must enter %d non-negative integers\n", rank);
       help = 1;
     }
  for (i=0; i < rank; i++)
   {
     scanf("%d", &exp[i]); 
     if (exp[i] < 0) 
     {
       printf("Negative exponents not allowed. Check syntax.\n");
       help = h;
       return 0;
     }
   }
  if (degree(exp) > mdegree) {
    printf("Degree is too big, or syntax error.  \n"); 
    help = h;
    return 0;
  }
  help = h;
  return 1;
}

int extra(exponent exp, int k)
{
  int deg, j, e;
  for(j = k , deg = 0; j < rank; j++)  deg += exp[j];
  if (!deg) return 0;
  for (j = deg - exp[k] + 1, e = extra(exp, k+1); j <= deg; j++)
    e += card_monos(j, rank - k - 1);
  return e;
}

/* Returns the number of monic monomials in k variables of degree exactly
deg  */
int card_monos(int deg, int k)
{
  return choose(deg + k - 1, deg);
}

int card_monos_atmost(int deg, int k)
{
  return choose(deg + k, deg);
}

int position(exponent exp)
{
  int  k;
  k = extra(exp,0) + card_monos_atmost(degree(exp) - 1, rank);
  if (k<0 || k >= total) {
      printf("Error in position. "); 
      printexp(exp);
      printf("\n");
      exit(0);
   }
  return k;
}

int readfield(void)
{
  int n;
  char c[2];
  if (help) printf("Enter a field element\n");
  scanf("%1s", c);
  if (! isdigit(c[0]))  
  {
    printf("You should have entered a field element as an integer.\n");
    return -1;
  }
  else
  {
    ungetc(c[0], stdin);
    scanf("%d",&n);
    return ((n % prime) + prime) % prime;
  }
}

void initpolynomial(polynomial f)
{
  int i;
  for(i=0; i < total; i++) f[i] = 0;
}

void initexponent(exponent exp)
{
  int i;
  for(i=0; i < rank; i++) exp[i] = 0;
}

int readpoly(polynomial f)
{
  exponent exp;  int n, k;
  initpolynomial(f); 
  if (help) printf("Enter the polynomial, terminate with zero coeft.\n");
  n = readfield();
  while (n)
  {
    if (n == -1) return 0;
    if (!readexp(exp)) return 0;
    f[position(exp)] = n;
    n = readfield();
  }
  return 1;
}

void advance(exponent exp)
{
  int i;
  for(i = rank-1; !exp[i] && i; i--);
  if (i) {
    exp[rank-1]=exp[i]-1;
    exp[i-1]++;
    if(i<rank-1) exp[i]=0;
  }
  else {
    exp[rank-1] = degree(exp) + 1;
    if (rank > 1) exp[0] = 0;
  }
}

void writepoly(polynomial f)
{
  int i; exponent exp; int first = 1;
  initexponent(exp);
  for(i = 0; i < total; i++)
  {  
    if (f[i]) 
    {
      if (first) first = 0;
      else printf("+"); 
      if (i == 0 || f[i] != 1) printf("%d",f[i]);
      printexp(exp);
    }
    advance(exp);
  }
  if (first) printf("%d",0);
}

void expsum(exponent exp1, exponent exp2, exponent exp3)
{
  int i;
  for (i = 0; i < rank; i++)
    exp3[i] = exp1[i] + exp2[i];
}
  
void polyprod(polynomial f, polynomial g, polynomial h)
{
  exponent exp1, exp2, exp3;
  int i=0, j, k, t;
  initexponent(exp1);
  initpolynomial(h);
  while((k = degree(exp1)) <= mdegree)
  { 
    if (f[i])
    {
      initexponent(exp2);
      j = 0;
      while (k + degree(exp2) <= mdegree)
      {
        if (g[j])
        {
          expsum(exp1, exp2, exp3);
          t = position(exp3);
          h[t] += f[i] * g[j];
          h[t] %= prime;
        }
        j++;
        advance(exp2);
      }
    }
    i++;
    advance(exp1);
  }
}

void polymult(polynomial f, polynomial g)
{
  int i;
  polynomial h;
  polyprod(f, g, h);
  for (i = 0; i < total; i++)
    f[i] = h[i];
}

void polypower(polynomial f, int n, polynomial h) 
{
  polynomial g;
  if (n < 0) {
    printf("Error in polypower %d\n", n); 
    exit(0);
  }
  if (!n)
  {
    initpolynomial(h);
    h[0] = 1;
    return;
  }
  if (n%2) 
  {
    polypower(f, n-1, g);
    polyprod(f, g, h);
    return;
  }
  polyprod(f, f, g);
  polypower(g, n/2, h);
}
    
void polysum(polynomial f, polynomial g)  /* f = f + g  */
{
  int i;
  for(i = 0; i < total; i++) f[i] = (f[i] + g[i]) % prime;
}
    
int polyequal(polynomial f, polynomial g)
{
  int i;
  for (i = 0; i < total; i++)
    if (f[i] != g[i]) return 0;
  return 1;
}