#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<setjmp.h>
#include<gmp.h>
#include<omp.h>
#include<math.h>
#define SPEED 64

typedef signed int SI;
typedef signed long long SLL;

#define DMX(a1,a2,d,r); { r=(SI)(((SLL)a1*(SLL)a2) % d); }
#define usermod(a,r,d); { r=(SI)((SLL)a % d); }

#pragma pack(1)

int *restrict IA,*restrict JA,*restrict JAP;

int MAXthreads;

mpz_t chtmp1, chtmp2, chtmp3;

mpz_t *restrict zmatrix;
mpz_t **restrict row = NULL, **restrict colm = NULL;
mpz_t *restrict resultZ = NULL;

/* definition */
#define MAXlett 65535
#define MAXline 655350000

int MAXmat;
int nonzero;

jmp_buf ma;

SI
ivM (SI a, SI CM, SI CMh)
{

  SI d, q, r;
  SI s, t, x;

  a = a+CM;
  d = CM;
  x = 0;
  s = 1;
  while (a != 0)
    {
      q = d / a;
      r = d % a;
      d = a;
      a = r;
      t = x - q * s;
      x = s;
      s = t;
    }

  if (d != 1)
    {
      printf ("Miller-Rabin method is false.\n");
      longjmp (ma, -1);
    }

  if (x>CMh)
    x=x-CM;
  else if (x<-CMh)
    x=x+CM;

  return x;

}

int
getline2 (char * restrict line, FILE * restrict fp)
{

  int c;
  int i;

  for (i = 0; i < MAXlett; i++)
    {
      c = fgetc (fp);
      if (c == EOF || c == '\n')
	{
	  line[i] = '\0';
	  return 1;
	}

      line[i] = (char) c;
    }

  return 0;

}

void
index2mnp (char * restrict line, int * restrict row, int * restrict colm, mpz_t * restrict pointer)
{

  int m;

  m = strcspn (line, ",\0");
  line[m] = '\0';
  *row = atoi (line);
  line = line + (m + 1);
  m = strcspn (line, ",\0");
  line[m] = '\0';
  *colm = atoi (line);
  line = line + (m + 1);

  mpz_set_str (*pointer, line, 10);

  return;
}

void
makemv (int i, SI * restrict vmatrix, SI * restrict xt, SLL * restrict polyx, int *restrict che0, int *restrict che)
{

  int j, k, row0;

  for (k = 0; k < MAXmat; k++)
    {
      polyx[k] = 0;
    }

  for (k = JAP[che0[i]]; k < JAP[che0[i]+1]; k++)
    {
      polyx[che[IA[k]]] = vmatrix[k];
    }

  for (j = i + 1; j < MAXmat; j++)
    {

      if (xt[j] != 0)
	{
	  for (k = JAP[che0[j]]; k < JAP[che0[j]+1]; k++)
	    {
	      row0 = che[IA[k]];
	      polyx[row0] = polyx[row0] 
		+ (SLL) vmatrix[k] * (SLL) xt[j];
	    }
	}

    }

}

void
mult1 (int k, SI t, SI * restrict poly0, SLL * restrict poly1)
{

  int j;

  if (t != 0)
    {

      t = - t;

      poly1[0] = (SLL) t *(SLL) poly0[0];

      for (j = 1; j < k; j++)
	{
	  poly1[j] = (SLL) t *(SLL) poly0[j] + (SLL) poly0[j - 1];
	}

      poly1[k] = (SLL) t + (SLL) poly0[k - 1];

    }
  else
    {

      poly1[0] = 0;

      for (j = 0; j < k; j++)
	{
	  poly1[j + 1] = poly0[j];
	}

    }

  return;
}

void
mult2 (int l, SI CM, SI CMh, SI t, SI u, SI * restrict poly0, SLL * restrict poly1)
{

  int k;

  DMX (t, u, CM, t);
  if (t>CMh)
    t=t-CM;
  else if (t<-CMh)
    t=t+CM;

  t = - t;

  for (k = 0; k < l; k++)
    {

      poly1[k] = (SLL) t *(SLL) poly0[k] + (SLL)poly1[k];

    }

  poly1[l] = (SLL) t + (SLL)poly1[l];

  return;
}

void
hess (SI CM, SI * restrict CM_result, SI (* restrict hmatrix)[MAXmat], SI * restrict xt, SLL * restrict polyx, SI * restrict vmatrix, int *restrict che0, int *restrict che)
{

  SI CMh=(CM-1)/2;
  int i, j, k, z;
  SI s;
  int t1, t2;

  mpz_t chtmpu;

  mpz_init (chtmpu);

  for (i = 0; i < MAXmat; i++)
    {
      che[i] = i;
    }

  for (i = 0; i < MAXmat; i++)
    {
      che0[i] = i;
    }

  for (i = 0; i < nonzero; i++)
    {
      vmatrix[i] = mpz_mod_ui (chtmpu, zmatrix[i], CM);
      if (vmatrix[i]>CMh)
        vmatrix[i]=vmatrix[i]-CM;
    }

  if (JAP[1] > 0 && IA[0] == 0)
    {
      hmatrix[0][0] = vmatrix[0];
      z = 1;
    }
  else
    {
      hmatrix[0][0] = 0;
      z = 0;
    }

  for (j = 1; j < MAXmat; j++)
    {
      xt[j] = 0;
    }
  for (j = z; j < JAP[1]; j++)
    {
      xt[che[IA[j]]] = vmatrix[j];
    }

  for (i = 0; i < MAXmat - 1; i++)
    {

      if (xt[i + 1] == 0)
	{

	  for (j = i + 2; j < MAXmat; j++)
	    {
	      if (xt[j] != 0)
		break;
	    }

	  if(j < MAXmat){

	    for (k = 0; k < i; k++)
	      {
		s = hmatrix[k][i + 1];
		hmatrix[k][i + 1] = hmatrix[k][j];
		hmatrix[k][j] = s;
	      }
	    
	    /* GYO KOUKANN */
	    t1 = che0[i + 1];
	    t2 = che0[j];
	    
	    z = che[t1];
	    che[t1] = che[t2];
	    che[t2] = z;
	    
	    che0[i + 1] = t2;
	    che0[j] = t1;

	    s = xt[i + 1];
	    xt[i + 1] = xt[j];
	    xt[j] = s;

	  }
	  else{

	    for (j = i + 1; j < MAXmat; j++)
	      {
		hmatrix[i][j] = 0;
	      }
	    
	    for (j = 0; j < MAXmat; j++)
	      {
		polyx[j] = 0;
	      }
	    
	    for (j = JAP[che0[i+1]]; j < JAP[che0[i+1]+1]; j++)
	      {
		polyx[che[IA[j]]] = vmatrix[j];
	      }
	    
	    goto NEXT;

	  }

	}

      hmatrix[i][i + 1] = xt[i + 1];
      s = ivM (xt[i + 1], CM, CMh);

      for (j = i + 2; j < MAXmat; j++)
	{
	  DMX (xt[j], s, CM, xt[j]);
	  if (xt[j]>CMh)
            xt[j]=xt[j]-CM;
          else if (xt[j]<-CMh)
            xt[j]=xt[j]+CM;

	  hmatrix[i][j] = xt[j];
	}

      makemv (i + 1, vmatrix, xt, polyx, che0, che);

    NEXT:

      usermod (polyx[0], hmatrix[i + 1][0], CM);
      if (hmatrix[i + 1][0]>CMh)
        hmatrix[i + 1][0]=hmatrix[i + 1][0]-CM;
      else if (hmatrix[i + 1][0]<-CMh)
        hmatrix[i + 1][0]=hmatrix[i + 1][0]+CM;

      for (k = 0; k < i + 1; k++)
	{

	  if (polyx[k + 1] != 0)
	    {
	      usermod (polyx[k + 1], s, CM);
	      if (s>CMh)
		s=s-CM;
	      else if (s<-CMh)
		s=s+CM;

	      hmatrix[i + 1][k + 1] = s;

	      if (s != 0)
		{
		  s = - s;
		}
	      else
		{
		  continue;
		}
	    }
	  else
	    {
	      hmatrix[i + 1][k + 1] = 0;
	      continue;
	    }

	  for (j = k + 2; j < MAXmat; j++)
	    {
	      polyx[j] = polyx[j] + (SLL) hmatrix[k][j] * (SLL) (s);
	    }
	}

      for (j = i + 2; j < MAXmat; j++)
	{
	  usermod (polyx[j], xt[j], CM);
	  if (xt[j]>CMh)
            xt[j]=xt[j]-CM;
          else if (xt[j]<-CMh)
            xt[j]=xt[j]+CM;

	}

    }

  for(k=0;k<MAXmat-1;k++){
    xt[k]=hmatrix[k][k+1];
  }

  for(k=1;k<MAXmat;k++){
    hmatrix[0][k]=hmatrix[k][k];
  }

  hmatrix[MAXmat-1][MAXmat-1] = - hmatrix[0][0];

  for (k = 1; k < MAXmat; k++)
    {
      mult1 (k, hmatrix[0][k], &hmatrix[MAXmat-k][MAXmat-k], polyx);
      s = 1;
      for (i = 1; i <= k; i++)
	{

	  if (xt[k - i] == 0)
	    break;

	  DMX (s, xt[k - i], CM, s);
	  if (s>CMh)
            s=s-CM;
          else if (s<-CMh)
            s=s+CM;

	  if (hmatrix[k][k - i] != 0)
	    {
	      mult2 (k - i, CM, CMh, s, hmatrix[k][k - i], &hmatrix[MAXmat-(k-i)][MAXmat-(k-i)], polyx);
	    }
	}

      for (i = 0; i <= k; i++)
	{
	  usermod (polyx[i], hmatrix[MAXmat-1-k][MAXmat-1-(k-i)], CM);

	  if (hmatrix[MAXmat-1-k][MAXmat-1-(k-i)]>CMh)
            hmatrix[MAXmat-1-k][MAXmat-1-(k-i)]=hmatrix[MAXmat-1-k][MAXmat-1-(k-i)]-CM;
          else if (hmatrix[MAXmat-1-k][MAXmat-1-(k-i)]<-CMh)
            hmatrix[MAXmat-1-k][MAXmat-1-(k-i)]=hmatrix[MAXmat-1-k][MAXmat-1-(k-i)]+CM;
	}
    }

  for (j = 0; j < MAXmat; j++)
    {
      CM_result[j] = hmatrix[0][j];
    }
      
  mpz_clear (chtmpu);

  return;
}

int
HE_main ()
{

  SI tmp;
  int i, j, CM_flag;
  SLL s;

  SI CM,*restrict CMD, *restrict CMhD, *restrict tD;
  int num_CM;
  mpz_t *restrict chtmps;
  mpz_t chtmpu;
  SI **restrict CM_result;
  SI (* restrict hmatrix)[MAXmat];
  SI * restrict xt;
  SLL * restrict polyx;
  SI * restrict vmatrix;
  int *restrict che0;
  int *restrict che;

  mpz_ui_pow_ui (chtmp2, 2, 65);
  mpz_sub_ui (chtmp2, chtmp2, 4);
  mpz_tdiv_q_ui (chtmp2, chtmp2, MAXmat);
  mpz_sqrt (chtmp2, chtmp2);
  s = mpz_get_ui (chtmp2)+1;

  if (s > (SLL) 1 << 30)
    s = (SLL) 1 << 30;

  num_CM = 0;
  CMD = NULL;

  if ((chtmps = (mpz_t * restrict ) malloc (sizeof (mpz_t))) == NULL)
    {
      printf ("Out of memory!\n");
      return 0;
    }

  mpz_init_set_ui (chtmps[0], 1);

  for (CM = s; CM >= 3; CM -= 2)
    {

      for (mpz_set_ui (chtmp1, CM);
	   mpz_probab_prime_p (chtmp1, 5) == 0;
	   CM--, mpz_set_ui (chtmp1, CM));

      num_CM = num_CM + 1;

      if ((CMD = (SI * restrict ) realloc (CMD, sizeof (SI)*num_CM)) == NULL)
	{
	  printf ("Out of memory!\n");
	  return 0;
	}

      CMD[num_CM - 1] = CM;

      if ((chtmps =
	   (mpz_t * restrict ) realloc (chtmps, sizeof (mpz_t)*(num_CM + 1))) == NULL)
	{
	  printf ("Out of memory!\n");
	  return 0;
	}

      mpz_init (chtmps[num_CM]);

      mpz_mul_ui (chtmps[num_CM], chtmps[num_CM - 1], CM);

      if (mpz_cmp (chtmps[num_CM], chtmp3) >= 0)
	break;

    }

  if(CM<2){
    printf ("p ga tarinai\n");
    return 0;
  }

  if ((CM_result = (SI ** restrict ) malloc (sizeof (SI * restrict )*num_CM)) == NULL)
    {
      printf ("Out of memory!\n");
      return 0;
    }

  for (i = 0; i < num_CM; i++)
    {
      if ((CM_result[i] = (SI * restrict ) malloc (sizeof (SI)*MAXmat)) == NULL)
	{
	  printf ("Out of memory!\n");
	  return 0;
	}
    }

  CM_flag = 0;

#pragma omp parallel private(hmatrix,xt,polyx,vmatrix,che0,che)
  {

    if ((vmatrix =
	 (SI * restrict ) malloc (sizeof (SI)*nonzero)) == NULL)
      {
	printf ("out of memory.\n");
	CM_flag = 1;
      }

    if ((hmatrix = (SI (* restrict )[MAXmat]) malloc (sizeof (SI)*MAXmat * MAXmat+SPEED)) == NULL)
      {
	printf ("out of memory.\n");
	CM_flag = 1;
      }
    
    if ((xt =
	 (SI * restrict ) malloc (sizeof (SI)*MAXmat + SPEED)) == NULL)
      {
	printf ("out of memory.\n");
	CM_flag = 1;
      }
    if ((polyx =
	 (SLL * restrict ) malloc (sizeof (SLL)*MAXmat + SPEED)) == NULL)
      {
	printf ("out of memory.\n");
	CM_flag = 1;
      }
    
    if ((che = (int * restrict ) malloc (sizeof (int)*MAXmat)) == NULL)
      {
	printf ("out of memory.\n");
	CM_flag = 1;
      }
    
    if ((che0 = (int * restrict ) malloc (sizeof (int)*MAXmat)) == NULL)
      {
	printf ("out of memory.\n");
	CM_flag = 1;
      }
    
#pragma omp for
    for (i = 0; i < num_CM; i++)
      {
	
	if (CM_flag == 1)
	  continue;
	
	hess (CMD[i], CM_result[i],hmatrix,xt,polyx,vmatrix,che0,che);
      }

    if (CM_flag==0){
      free(hmatrix);
      free(xt);
      free(polyx);
      free (vmatrix);
      free (che);
      free (che0);
    }

  }

  if (CM_flag == 1)
    {
      printf ("HE falt!\n");
      return 0;
    }

  if ((CMhD = (SI * restrict ) malloc (sizeof (SI)*num_CM)) == NULL)
    {
      printf ("Out of memory!\n");
      return 0;
    }

  for (i = 0; i < num_CM; i++)
    {
      CMhD[i] = (CMD[i] - 1) / 2;
    }

  if ((tD = (SI * restrict ) malloc (sizeof (SI)*num_CM)) == NULL)
    {
      printf ("Out of memory!\n");
      return 0;
    }

  for (i = 0; i < num_CM; i++)
    {
      tD[i] = ivM (mpz_mod_ui (chtmp1, chtmps[i], CMD[i]), CMD[i], CMhD[i]);
    }

#pragma omp parallel private(chtmpu,tmp,i)
  {

    mpz_init (chtmpu);

#pragma omp for
    for (j = 0; j < MAXmat; j++)
      {

	for (i = 0; i < num_CM; i++)
	  {

	    tmp=(CM_result[i][j] - (SI) mpz_mod_ui (chtmpu, resultZ[j], CMD[i])) % CMD[i];

	    if (tmp!=0){
	      tmp=((SLL)tmp*(SLL)tD[i]) % CMD[i];
	      if (tmp > CMhD[i])
                tmp = tmp - CMD[i];
              else if (tmp < -CMhD[i])
                tmp = tmp + CMD[i];
	      mpz_mul_si (chtmpu, chtmps[i], (long int) (int) tmp);
	      mpz_add (resultZ[j], resultZ[j], chtmpu);
	    }

	  }

      }

    mpz_clear (chtmpu);
  }

  return MAXmat;

}

int
matrix2Hadamard ()
{

  int i,j;
  int row0, colm0;
  int flag[MAXmat],flag0;
  mpz_t chtmpu;
  int mynumber;

  for(i=0;i<MAXmat;i++)
    flag[i]=0;

  flag0=0;

#pragma omp parallel private(chtmpu,mynumber,j,row0,colm0)
  {

    mpz_init(chtmpu);

    mynumber=omp_get_thread_num();

    if((row[mynumber]=(mpz_t * restrict)malloc(sizeof(mpz_t)*MAXmat))==NULL){
      printf("Out of memory!\n");
      flag0=1;
    }

    if(flag0==0){
      for(j=0;j<MAXmat;j++){
	mpz_init_set_si(row[mynumber][j],0);
      }
    }

    if((colm[mynumber]=(mpz_t * restrict)malloc(sizeof(mpz_t)*MAXmat))==NULL){
      printf("Out of memory!\n");
      flag0=1;
    }

    if(flag0==0){
      for(j=0;j<MAXmat;j++){
	mpz_init_set_si(colm[mynumber][j],0);
      }
    }

#pragma omp for
    for (i = 0; i < nonzero; i++)
      {

        if (flag0==1)
          continue;

        row0 = IA[i];
        colm0 = JA[i];

        if (row0 != colm0)
          {
            mpz_pow_ui (chtmpu, zmatrix[i], 2);
          }
        else
          {
            flag[row0]=1;

            mpz_abs (chtmpu, zmatrix[i]);
            mpz_add_ui (chtmpu, chtmpu, 1);
            mpz_pow_ui (chtmpu, chtmpu, 2);
          }

        mpz_add (row[mynumber][row0], row[mynumber][row0], chtmpu);
        mpz_add (colm[mynumber][colm0], colm[mynumber][colm0], chtmpu);
      }

    mpz_clear(chtmpu);
  }

  if(flag0==1)
    return 0;

#pragma omp for private(i)
  for(j=0;j<MAXmat;j++){
    for(i=1;i<MAXthreads;i++){
      mpz_add (row[0][j], row[0][j], row[i][j]);
      mpz_clear(row[i][j]);
    }
    for(i=1;i<MAXthreads;i++){
      mpz_add (colm[0][j], colm[0][j], colm[i][j]);
      mpz_clear(colm[i][j]);
    }
  }

  for(i=0;i<MAXmat;i++){
    if (flag[i]==0){
      mpz_add_ui (row[0][i], row[0][i], 1);
      mpz_add_ui (colm[0][i], colm[0][i], 1);
    }
  }

  mpz_set_si (chtmp1, 1);
  for (row0 = 0; row0 < MAXmat; row0++)
    {
      mpz_mul (chtmp1, chtmp1, row[0][row0]);
      mpz_clear(row[0][row0]);
    }

  mpz_set_si (chtmp2, 1);
  for (colm0 = 0; colm0 < MAXmat; colm0++)
    {
      mpz_mul (chtmp2, chtmp2, colm[0][colm0]);
      mpz_clear(colm[0][colm0]);
    }

  for(i=0;i<MAXthreads;i++){
    free(row[i]);
  }

  free(row);

  for(i=0;i<MAXthreads;i++){
    free(colm[i]);
  }

  free(colm);

  if (mpz_cmp (chtmp1, chtmp2) > 0)
    mpz_set (chtmp1, chtmp2);

  mpz_sqrtrem (chtmp3, chtmp2, chtmp1);

  if (mpz_cmp_ui (chtmp2, 0) != 0)
    mpz_add_ui (chtmp3, chtmp3, 1);

  mpz_mul_ui (chtmp3, chtmp3, 2);
  mpz_add_ui (chtmp3, chtmp3, 1);

  return 1;
}

int
main (int argc, char ** restrict argv)
{

  int i, k, l, t;
  long int bit_bound;
  int row0, colm0;

  FILE * restrict fp;
  char * restrict line;
  int line0;
  int key;

  mpz_init (chtmp1);
  mpz_init (chtmp2);
  mpz_init (chtmp3);

  MAXthreads = omp_get_max_threads ();

  if (argc < 3)
    {
      printf ("input file name.\n");
      return 0;
    }

  fp = fopen (argv[1], "r");
  if ((line = (char * restrict ) malloc (sizeof (char)*MAXlett)) == NULL)
    {
      printf ("out of momory.\n");
      return 0;
    }

  k = 0;
  l = 0;
  printf ("read file.\n");
  for (line0 = 0; line0 < MAXline; line0++)
    {

      if (getline2 (line, fp) != 1)
	{
	  printf ("too many letters in one line.\n");
	  return 0;
	}

      if (line[0] == '\0')
	{
	  printf ("end read file.\n");
	  break;
	}

      switch (line[0])
	{
	case 'Z':
	  nonzero = atoi (line + 1);

	  if ((zmatrix =
	       (mpz_t * restrict ) malloc (sizeof (mpz_t)*nonzero)) == NULL)
	    {
	      printf ("out of memory.\n");
	      return 0;
	    }

	  for(i=0;i<nonzero;i++){
	    mpz_init(zmatrix[i]);
	  }

	  if ((IA =
	       (int * restrict ) malloc (sizeof (int)*nonzero)) == NULL)
	    {
	      printf ("out of memory.\n");
	      return 0;
	    }

	  if ((JA =
	       (int * restrict ) malloc (sizeof (int)*nonzero)) == NULL)
	    {
	      printf ("out of memory.\n");
	      return 0;
	    }

	  break;

	case 'N':
	  MAXmat = atoi (line + 1);

	  if ((JAP =
	       (int * restrict ) malloc (sizeof (int)*(MAXmat+1))) == NULL)
	    {
	      printf ("out of memory.\n");
	      return 0;
	    }

	  JAP[0] = 0;

	  if ((row = (mpz_t ** restrict ) malloc (sizeof (mpz_t * restrict)*MAXthreads)) == NULL)
	    {
	      printf ("out of memory.\n");
	      return 0;
	    }

	  if ((colm = (mpz_t ** restrict ) malloc (sizeof (mpz_t * restrict)*MAXthreads)) == NULL)
	    {
	      printf ("out of memory.\n");
	      return 0;
	    }

	  if ((resultZ = (mpz_t * restrict ) malloc (sizeof (mpz_t)*MAXmat)) == NULL)
	    {
	      printf ("out of memory.\n");
	      return 0;
	    }

	  for (i = 0; i < MAXmat; i++)
	    {
	      mpz_init_set_si (resultZ[i], 0);
	    }

	  break;

	case 'I':
	  index2mnp (line + 1, &row0, &colm0, &zmatrix[l]);

	  if(row0!=k){
	    if (k+1<row0) printf("not row full rank\n");
	    for(k=k+1;k<=row0;k++){
	      JAP[k] = l;
	    }
	    k=row0;
	  }

	  IA[l] = colm0;
	  JA[l] = row0;
	  l++;

	  break;
	}
    }
  if (k+1<MAXmat) printf("not row full rank\n");
  for(k=k+1;k<=MAXmat;k++){
    JAP[k] = l;
  }
  fclose (fp);

  if (line0 == MAXline)
    {
      printf ("too many line.\n");
      return 0;
    }

  printf ("Calcurate Hadamard bound to C[n]\n");
  fflush (stdout);
  if(matrix2Hadamard ()==0)
    return 0;

  mpz_get_d_2exp (&bit_bound, chtmp3);
  printf ("bit bound=%ld\n", bit_bound);
  fflush (stdout);

  key = setjmp (ma);

  if (key == -1)
    return 0;

  if ((t = HE_main ()) == 0)
    return 0;

  fp = fopen (argv[2], "w");

  if (t == 1)
    gmp_fprintf (fp, "(1)*x");
  else
    gmp_fprintf (fp, "(1)*x^%d", t);

  for (i = t - 1; i >= 0; i--)
    {

      if (mpz_cmp_si (resultZ[i], 0) != 0)
	{

	  switch (i)
	    {
	    case 0:
	      gmp_fprintf (fp, "+(%Zd)", resultZ[i]);
	      break;
	    case 1:
	      if (mpz_cmp_si (resultZ[i], 1) == 0)
		gmp_fprintf (fp, "+x");
	      else
		gmp_fprintf (fp, "+(%Zd)*x", resultZ[i]);
	      break;
	    default:
	      if (mpz_cmp_si (resultZ[i], 1) == 0)
		gmp_fprintf (fp, "+x^%d", i);
	      else
		gmp_fprintf (fp, "+(%Zd)*x^%d", resultZ[i], i);
	    }

	}

    }

  gmp_fprintf (fp, "\n");

  fclose (fp);

  return 0;

}