#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<setjmp.h>
#include<gmp.h>
#include<sys/types.h>
#include<unistd.h>
#include<omp.h>

#define SPEED 64

typedef signed int SI;
typedef signed long long SLL;

#define DM(a1,a2,d,r); r=(SI)(((SLL)(a1)*(SLL)(a2)) % d);

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

#pragma pack(1)

int *restrict IAP,*restrict JA;

int MAXthreads;

mpz_t chtmp1, chtmp2;

mpz_t *restrict zmatrix;
mpz_t *restrict resultZ = NULL;

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

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
makesind (SI * restrict sind, int i, SI CM, SI CMh, SI (* restrict vmatrix)[MAXmat], SI  * restrict ut, SI * restrict vt, SI * restrict v2t, SLL * restrict u_tmp)
{

  SLL tmp1, tmp3;
  int j, k;

  for(j=0;j<MAXmat;j++){
    u_tmp[j]=0;
  }

  tmp3 = 0;
  for (j = 0; j < MAXmat; j++)
    {
      tmp1 = 0;
      for (k = 0; k < MAXmat; k++)
	{
	  tmp1 = (SLL) vmatrix[j][k] * (SLL) vt[k] + tmp1;
	  u_tmp[k] =(SLL) vmatrix[j][k] * (SLL) ut[j] + u_tmp[k];
	}
      usermod (tmp1, v2t[j], CM);
      if (v2t[j]>CMh)
        v2t[j]=v2t[j]-CM;
      else if (v2t[j]<-CMh)
        v2t[j]=v2t[j]+CM;

      tmp3 = (SLL) v2t[j] * (SLL) ut[j] + tmp3;
    }
  usermod (tmp3, sind[i], CM);
  if (sind[i]>CMh)
    sind[i]=sind[i]-CM;
  else if (sind[i]<-CMh)
    sind[i]=sind[i]+CM;

  tmp3 = 0;
  for(j=0;j<MAXmat;j++){
    usermod (u_tmp[j], ut[j], CM);
    if (ut[j]>CMh)
      ut[j]=ut[j]-CM;
    else if (ut[j]<-CMh)
      ut[j]=ut[j]+CM;

    tmp3 = (SLL) v2t[j] * (SLL) ut[j] + tmp3;
  }
  usermod (tmp3, sind[i+1], CM);
  if (sind[i+1]>CMh)
    sind[i+1]=sind[i+1]-CM;
  else if (sind[i+1]<-CMh)
    sind[i+1]=sind[i+1]+CM;

  return;
}

void
makesind_last (SI * restrict sind, int i,SI CM, SI CMh,SI (* restrict vmatrix)[MAXmat], SI * restrict ut, SI * restrict vt, SI * restrict v2t)
{

  SLL tmp1, tmp3;
  int j, k;

  tmp3 = 0;
  for (j = 0; j < MAXmat; j++)
    {
      tmp1 = 0;
      for (k = 0; k < MAXmat; k++)
        {
          tmp1 = (SLL) vmatrix[j][k] * (SLL) vt[k] + tmp1;
        }
      usermod (tmp1, v2t[j], CM);
      if (v2t[j]>CMh)
        v2t[j]=v2t[j]-CM;
      else if (v2t[j]<-CMh)
        v2t[j]=v2t[j]+CM;

      tmp3 = (SLL) v2t[j] * (SLL) ut[j] + tmp3;
    }
  usermod (tmp3, sind[i], CM);
  if (sind[i]>CMh)
    sind[i]=sind[i]-CM;
  else if (sind[i]<-CMh)
    sind[i]=sind[i]+CM;

  return;
}

size_t g1, g2, g3, h;

int
ber (SI CM, SI * restrict CM_result)
{

  SI CMh=(CM-1)/2;
  int i, j, k;
  int lc = 0;
  int lb = 0;
  int t;
  SI D, DM;
  SLL tmp1;
  int pp,qq;

  SI (*restrict vmatrix)[MAXmat];
  SI *restrict polya, *restrict polyb, *restrict polyc;
  SI *restrict ut, *restrict vt, *restrict v2t;
  SI *restrict sind;

  SLL *restrict u_tmp;

  mpz_t chtmpu;

  mpz_init (chtmpu);

  if ((vmatrix = (SI (*restrict)[MAXmat]) calloc (MAXmat * MAXmat,sizeof (SI))) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((polya = (SI * restrict ) malloc (g2)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((polyb = (SI * restrict ) malloc (g2)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((polyc = (SI * restrict ) malloc (g1)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((ut = (SI * restrict ) malloc (g2)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((u_tmp = (SLL * restrict ) malloc (g3)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((vt = (SI * restrict ) malloc (g2)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((v2t = (SI * restrict ) malloc (g2)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  if ((sind = (SI * restrict ) malloc (h)) == NULL)
    {
      printf ("out of memory.\n");
      return (0);
    }

  for (i = 0; i < MAXmat; i++){
    for(j=IAP[i];j<IAP[i+1];j++){
      vmatrix[i][JA[j]] = mpz_mod_ui (chtmpu, zmatrix[j], CM);

      if (vmatrix[i][JA[j]]>CMh)
        vmatrix[i][JA[j]]=vmatrix[i][JA[j]]-CM;
    }
  }

  for (i = 0; i < MAXmat; i++)
    {
      tmp1 = (SLL) rand () * (SLL) rand ();
      usermod (tmp1, ut[i], CM);
      if (ut[i]>CMh)
        ut[i]=ut[i]-CM;
      else if (ut[i]<-CMh)
        ut[i]=ut[i]+CM;

      tmp1 = (SLL) rand () * (SLL) rand ();
      usermod (tmp1, vt[i], CM);
      if (vt[i]>CMh)
        vt[i]=vt[i]-CM;
      else if (vt[i]<-CMh)
        vt[i]=vt[i]+CM;
    }

  tmp1 = 0;
  for (i = 0; i < MAXmat; i++)
    {
      tmp1 = tmp1 + (SLL) ut[i] * (SLL) vt[i];
    }
  usermod (tmp1, sind[0], CM);
  if (sind[0]>CMh)
    sind[0]=sind[0]-CM;
  else if (sind[0]<-CMh)
    sind[0]=sind[0]+CM;

  polyc[0] = 1;
  pp = 0;
  qq = 0;

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

      if (qq == 1)
        {

	  if(i!=2*MAXmat-1){
	    if(pp==0)
	      makesind (sind, i, CM, CMh, vmatrix, ut, vt, v2t, u_tmp);
	    else
	      makesind (sind, i, CM, CMh, vmatrix, ut, v2t, vt, u_tmp);
	  }
	  else{
	    if(pp==0)
	      makesind_last (sind, i, CM, CMh, vmatrix, ut, vt, v2t);
	    else
	      makesind_last (sind, i, CM, CMh, vmatrix, ut, v2t, vt);
	  }

          pp=1-pp;
        }

      tmp1 = 0;
      for (k = 0; k < lb; k++)
        {
          tmp1 = (SLL) polyb[k] * (SLL) sind[i - (k+1)] + tmp1;
        }
      tmp1=tmp1+sind[i];
      usermod (tmp1, D, CM);
      if (D>CMh)
        D=D-CM;
      else if (D<-CMh)
        D=D+CM;

      if (D == 0)
        {
          break;
        }

      DM = - D;

      if(i<2*MAXmat-1){
        for (k = 0; k < lb; k++)
          {
            DMARX (DM, polyc[k], polyb[k], CM, polya[k]);
	    if (polya[k]>CMh)
              polya[k]=polya[k]-CM;
            else if (polya[k]<-CMh)
              polya[k]=polya[k]+CM;
          }
      }
      else{
        for (k = 0; k < lb; k++)
          {
            DMARX (DM, polyc[k], polyb[k], CM, polyb[k]);
	    if (polyb[k]>CMh)
              polyb[k]=polyb[k]-CM;
            else if (polyb[k]<-CMh)
              polyb[k]=polyb[k]+CM;
          }
        break;
      }

      if (lb == lc)
        {
          DMX (DM, polyc[lc], CM, polya[lc]);
	  if (polya[lc]>CMh)
            polya[lc]=polya[lc]-CM;
          else if (polya[lc]<-CMh)
            polya[lc]=polya[lc]+CM;
        }

      lc = lc + 1;

      if (lb < lc)
        {
          DM = ivM (D, CM, CMh);

	  polyc[0]=DM;
          for (k = 0; k < lb; k++)
            {
              DMX (DM, polyb[k], CM, polyc[k+1]);
	      if (polyc[k+1]>CMh)
                polyc[k+1]=polyc[k+1]-CM;
              else if (polyc[k+1]<-CMh)
                polyc[k+1]=polyc[k+1]+CM;
            }
          t = lc;
          lc = lb;
          lb = t;
        }
      else
        {
          for (k = lc - 1; k >= 0; k--)
            {
              polyc[k + 1] = polyc[k];

            }
          polyc[0] = 0;
        }

      for (k = 0; k < lb; k++)
        {
          polyb[k] = polya[k];
        }

      qq=1-qq;
    }

  for (j = 0; j < lb; j++)
    {
      CM_result[j] = polyb[lb - j-1];
    }

  free (vmatrix);

  free (polya);
  free (polyb);
  free (polyc);
  free (u_tmp);
  free (ut);
  free (vt);
  free (v2t);
  free (sind);

  mpz_clear (chtmpu);

  return lb;
}

int
BM_main ()
{

  SI tmp;
  int i, j, k = 0, l, m, lb0 = -1, lb1, flag, count, pp;
  SLL s;

  mpz_t *restrict chtmpy, *restrict chtmpz;

  int *restrict lb;
  SI CM,*restrict CMD, *restrict CMhD, *restrict tD;
  int mynumber;
  mpz_t *restrict chtmps;
  mpz_t chtmpu;
  SI **restrict CM_result;

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

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

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

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

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

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

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

  for (i = 0; i < MAXthreads; i++)
    {
      mpz_init (chtmps[i]);
    }

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

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

  CM = s;

  for (;;)
    {

      for (i = 0; CM >= 3 && i < MAXthreads; i++, CM -= 2)
	{

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

	  CMD[i] = CM;
	}

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

      flag = 0;

#pragma omp parallel private(mynumber)
      {

	mynumber = omp_get_thread_num ();

	if ((lb[mynumber * SPEED] =
             ber (CMD[mynumber], CM_result[mynumber])) == 0)
          flag=1;

      }

      if (flag == 1)
	{
	  return 0;
	}

      lb1 = -1;
      for (i = 0; i < MAXthreads; i++)
	{
	  if (lb[i * SPEED] > lb1)
	    {
	      lb1 = lb[i * SPEED];
	      m = i;
	    }
	}

      if (lb1 == lb0)
        {
          mpz_set (chtmps[m], chtmp2);
          CMhD[m] = (CMD[m] - 1) / 2;
          tD[m] = ivM (mpz_mod_ui (chtmp1, chtmps[m], CMD[m]), CMD[m], CMhD[m]);
        }
      else if (lb1 > lb0)
        {
          mpz_set_ui (chtmps[m], 1);
          CMhD[m] = (CMD[m] - 1) / 2;
        }
      else
        {
          continue;
        }

      count=1;
      k = m;
      for (i = m + 1; i < MAXthreads; i++)
        {
          if (lb[i * SPEED] == lb1)
            {
              count++;
              mpz_mul_ui (chtmps[i], chtmps[k], CMD[k]);
              k = i;
              CMhD[i] = (CMD[i] - 1) / 2;
              tD[i] = ivM (mpz_mod_ui (chtmp1, chtmps[i], CMD[i]), CMD[i], CMhD[i]);
            }
        }
      mpz_mul_ui (chtmp2, chtmps[k], CMD[k]);

      if (lb1 == lb0)
        {
          flag=0;
        }
      else
        {
          if(count!=1)
            flag=0;
          else
            flag=1;
        }

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

        mpz_init (chtmpu);

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

            if (lb1 == lb0)
              {
                i = m;
              }
            else
              {
                tmp = CM_result[m][j];
		if (tmp > CMhD[m])
                  tmp = tmp - CMD[m];
                else if (tmp < -CMhD[m])
                  tmp = tmp + CMD[m];
                mpz_set_si (resultZ[j], (long int) (int) tmp);
                i = m+1;
              }

            for (; i < MAXthreads; i++)
              {
                if (lb[i * SPEED] == lb1){
		  tmp = (CM_result[i][j] - (SI) mpz_mod_ui (chtmpu, resultZ[j],CMD[i])) % CMD[i];

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

        mpz_clear (chtmpu);

      }

      if (flag == 0)
        break;
      
      lb0 = lb1;
    }

  flag = 0;

  printf ("check\n");

#pragma omp parallel private(j,k,l,pp,chtmpy,chtmpz)
  {

    if ((chtmpy =
	 (mpz_t * restrict ) malloc (sizeof (mpz_t)*MAXmat)) == NULL)
      {
	printf ("out of memory.\n");
	flag=1;
      }

    if (flag==0){
      for (j = 0; j < MAXmat; j++)
	{
	  mpz_init (chtmpy[j]);
	}
    }

    if ((chtmpz =
	 (mpz_t * restrict ) malloc (sizeof (mpz_t)*MAXmat)) == NULL)
      {
	printf ("out of memory.\n");
	flag=1;
      }
    
    if (flag==0){
      for (j = 0; j < MAXmat; j++)
	{
	  mpz_init (chtmpz[j]);
	}
    }

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

	if(flag==1)
	  continue;

	for (j = 0; j < MAXmat; j++)
	  {
	    mpz_set_si (chtmpy[j], 0);
	  }
	mpz_set_si (chtmpy[i], 1);

	pp = 1;
	for (j = lb1 - 1; j >= 0; j--)
	  {

	    pp = 1 - pp;
	    if (pp == 0)
	      {

		for (k = 0; k < MAXmat; k++)
		  {
		    mpz_set_si (chtmpz[k], 0);
		  }

		for (k = 0; k < MAXmat; k++)
                  {
                    for(l = IAP[k]; l < IAP[k+1]; l++)
                      {
                        mpz_addmul (chtmpz[k], zmatrix[l],
                                    chtmpy[JA[l]]);
                      }
                  }
		mpz_add (chtmpz[i], resultZ[j], chtmpz[i]);
	      }
	    else
	      {

		for (k = 0; k < MAXmat; k++)
		  {
		    mpz_set_si (chtmpy[k], 0);
		  }

		for (k = 0; k < MAXmat; k++)
                  {
                    for(l = IAP[k]; l < IAP[k+1]; l++)
                      {
                        mpz_addmul (chtmpy[k], zmatrix[l], chtmpz[JA[l]]);
                      }
                  }
		mpz_add (chtmpy[i], resultZ[j], chtmpy[i]);
	      }
	  }

	if (pp == 0)
	  {
	    for (j = 0; j < MAXmat; j++)
	      {
		if (mpz_cmp_si (chtmpz[j], 0) != 0)
		  {
		    flag = 1;
		    printf ("error\n");
		    break;
		  }
	      }
	  }
	else
	  {
	    for (j = 0; j < MAXmat; j++)
	      {
		if (mpz_cmp_si (chtmpy[j], 0) != 0)
		  {
		    flag = 1;
		    printf ("error\n");
		    break;
		  }
	      }
	  }

      }

  }

  if (flag == 1)
    {
      return 0;
    }

  return lb1;
}

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

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

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

  mpz_init (chtmp1);
  mpz_init (chtmp2);

  srand (getpid ());

  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 ((JA =
               (int * restrict ) malloc (sizeof (int)*nonzero)) == NULL)
            {
              printf ("out of memory.\n");
              return 0;
            }

	  break;

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

	  g1 = sizeof (SI)*(MAXmat + 1);
	  g2 = sizeof (SI)*MAXmat;
	  g3 = sizeof (SLL)*MAXmat;
	  h = sizeof (SI)*2*MAXmat;

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

          IAP[0] = 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++){
              IAP[k] = l;
            }
            k=row0;
          }

          JA[l] = colm0;
          l++;

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

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

  key = setjmp (ma);

  if (key == -1)
    return 0;

  if ((t = BM_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;

}