#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <setjmp.h> 
#include <gmp.h>
#include <unistd.h>
#include <omp.h>
#include <math.h>
#include "kimura_common.h"

#define NB2 64
#define UP0 10
#define UP1 10
#define SPEED 64

#pragma pack(1)

int *restrict IAP,*restrict IA,*restrict JA;

mpz_t chtmp1,chtmp2,chtmp3,chtmp4;

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

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

int MAXmat;
int nonzero;
int NB;

int MAXthreads;

jmp_buf ma;

UI invm2_(UI a,UI CM){

  UI d, q, r;
  int s, t, x;

  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");
  }

  return (x+CM);
}

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 A_mul_v_(double (* restrict A)[MAXmat],double (* restrict B)[MAXmat],double (* restrict C)[MAXmat],double (* restrict WORKD)[MAXmat],SI CM, SI CMh){

  int j,k;
  int tmp0,tmp1,tmp2;
  double ALPHA=1.0,BETA=0.0;
  SI v;

  for(tmp0=NB;tmp0<=MAXmat-NB;tmp0=tmp0+NB){

    dgemm_("N","N",&MAXmat,&NB,&MAXmat,&ALPHA,A,&MAXmat,&B[tmp0-NB][0],&MAXmat,&BETA,&B[tmp0][0],&MAXmat);

    for(k=0;k<NB;k++){
      for(j=0;j<MAXmat;j++){
	v=(SLL)B[tmp0+k][j] % CM;
        if (v>CMh)
          v=v-CM;
        else if (v<-CMh)
          v=v+CM;

        B[tmp0+k][j]=v;
      }
    }
  }

  if(tmp0==MAXmat){
    dgemm_("N","N",&MAXmat,&NB,&MAXmat,&ALPHA,A,&MAXmat,&B[tmp0-NB][0],&MAXmat,&BETA,C,&MAXmat);

    for(k=0;k<NB;k++){
      for(j=0;j<MAXmat;j++){
	v=(SLL)C[k][j] % CM;
        if (v>CMh)
          v=v-CM;
        else if (v<-CMh)
          v=v+CM;

        C[k][j]=v;
      }
    }

  }
  else{
    dgemm_("N","N",&MAXmat,&NB,&MAXmat,&ALPHA,A,&MAXmat,&B[tmp0-NB][0],&MAXmat,&BETA,WORKD,&MAXmat);

    tmp1=MAXmat-tmp0;

    for(k=0;k<tmp1;k++){
      for(j=0;j<MAXmat;j++){
	v=(SLL)WORKD[k][j] % CM;
        if (v>CMh)
          v=v-CM;
        else if (v<-CMh)
          v=v+CM;

        B[k+tmp0][j]=v;
      }
    }

    tmp2=NB-tmp1;

    for(k=0;k<tmp2;k++){
      for(j=0;j<MAXmat;j++){
	v=(SLL)WORKD[k+tmp1][j] % CM;
        if (v>CMh)
          v=v-CM;
        else if (v<-CMh)
          v=v+CM;

        C[k][j]=v;
      }
    }

    dgemm_("N","N",&MAXmat,&tmp1,&MAXmat,&ALPHA,A,&MAXmat,&B[tmp0][0],&MAXmat,&BETA,&C[tmp2][0],&MAXmat);

    for(k=0;k<tmp1;k++){
      for(j=0;j<MAXmat;j++){
	v=(SLL)C[k+tmp2][j] % CM;
        if (v>CMh)
          v=v-CM;
        else if (v<-CMh)
          v=v+CM;

        C[k+tmp2][j]=v;
      }
    }

  }

}

UI makesind(UI (* restrict vmatrix)[MAXmat],UI * restrict ut,UI * restrict vt,UI * restrict v2t,UI CM,ULL CMX){

  int i,j;
  ULL tmp;
  UI tmp0;
  int NN;

  NN=MAXmat-NB;

  for(j=0;j<NN;j++){
    v2t[j]=vt[NB+j];
  }

  for(i=0;i<NB;i++){
    tmp=0;
    for(j=0;j<MAXmat;j++){
      tmp=tmp+(ULL)vmatrix[i][j]*(ULL)vt[j];
    }
    usermod(tmp,v2t[NN+i],CM,CMX);
  }

  tmp=0;
  for(i=0;i<MAXmat;i++){
    tmp=tmp+(ULL)v2t[i]*(ULL)ut[i];
  }

  usermod(tmp,tmp0,CM,CMX);

  return(tmp0);
}

int ber(UI * restrict sind,UI * restrict polya,UI * restrict polyb,UI * restrict polyc,UI (* restrict vmatrix)[MAXmat],UI * restrict ut,UI * restrict vt,UI * restrict v2t,UI CM){

  UI TWO_CM=2*CM;
  ULL CMX=(ULL)(-1)/CM;
  int i,k;
  int lc=0;
  int lb=0;
  int t;
  UI D,D0;
  ULL tmp1;
  int pp;

  for(i=0;i<MAXmat;i++){
    tmp1=(ULL)rand()*(ULL)rand();
    usermod(tmp1,ut[i],CM,CMX);
    tmp1=(ULL)rand()*(ULL)rand();
    usermod(tmp1,vt[i],CM,CMX);
  }

  tmp1=0;
  for(i=0;i<MAXmat;i++){
    tmp1=tmp1+(ULL)ut[i]*(ULL)vt[i];
  }
  usermod(tmp1,sind[0],CM,CMX);
  polyc[0]=1;
  pp=0;

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

    if(i!=0){

      sind[i]=((pp=1-pp)==1) ? makesind(vmatrix,ut,vt,v2t,CM,CMX) :
	makesind(vmatrix,ut,v2t,vt,CM,CMX);
      
    }

    tmp1=0;
    for(k=0;k<lb;k++){
      tmp1=(ULL)polyb[k]*(ULL)sind[i-(k+1)]+tmp1;
    }
    tmp1=tmp1+sind[i];
    usermod(tmp1,D,CM,CMX);

    if(D==0 || D==CM){
      return 0;
    }

    D0=TWO_CM-D;

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

    lc=lc+1;

    if(lb<lc){
      D0=invm2_(D,CM);
      polyc[0]=D0;
      for(k=0;k<lb;k++){
	DMX(D0,polyb[k],CM,polyc[k+1],CMX);
      }
      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];
    }
  }

  return MAXmat;
}

int BM_main(){

  UI tmp;
  int i,j,k=0,l,m,n;

  SI CM;
  SI CMh;
  SI v;
  UI * restrict tD,s1,s2;
  SI * restrict CMD;
  SI * restrict CMhD;

  int CM_flag;
  int num_CM;

  UI ** restrict CM_result;
  int INFO,IONE=1;

  mpz_t * restrict chtmps;
  mpz_t chtmpu;

  double (* restrict A)[MAXmat];
  double (* restrict B)[MAXmat];
  double (* restrict C)[MAXmat];
  SI (* restrict WORK)[MAXmat];
  double (* restrict WORKD3)[MAXmat];
  UI (* restrict WORK2)[MAXmat];
  UI * restrict sind;
  UI * restrict polya;
  UI * restrict polyb;
  UI * restrict polyc;
  UI * restrict ut;
  UI * restrict vt;
  UI * restrict v2t;
  int * restrict IPIV;
  double WORKD[NB0][NB0];
  double * restrict WORKD2;

  mpz_ui_pow_ui(chtmp2,2,53);
  mpz_tdiv_q_ui(chtmp2,chtmp2,MAXmat);
  mpz_sqrt(chtmp2,chtmp2);
  s1=2*mpz_get_ui(chtmp2)+1;

  mpz_ui_pow_ui(chtmp2,2,64);
  mpz_sub_ui(chtmp2,chtmp2,1);
  mpz_tdiv_q_ui(chtmp2,chtmp2,MAXmat+1);
  mpz_sqrt(chtmp2,chtmp2);
  s2=(mpz_get_ui(chtmp2)+1)/2;

  if(s2<s1) s1=s2;

  mpz_set_ui(chtmp2,1);

  CM=s1;

  m=0;

  for(;;CM-=2){

    mpz_set(chtmp4,chtmp2);

    CMD=NULL;

    num_CM=0;

    for(;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;

      mpz_mul_ui(chtmp4,chtmp4,CM);

      if(mpz_cmp(chtmp4,chtmp3)>=0)
	break;

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

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

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

    CM_flag=0;
    
#pragma omp parallel private(chtmpu,INFO,k,j,l,A,B,C,WORK,WORK2,WORKD,WORKD3,CM,CMh,v,polya,polyb,polyc,ut,vt,v2t,sind,IPIV,WORKD2)
    {

      mpz_init(chtmpu);

      if((A=(double (* restrict)[MAXmat])
	  malloc(sizeof(double)*MAXmat*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((B=(double (* restrict)[MAXmat])
	  malloc(sizeof(double)*MAXmat*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((C=(double (* restrict)[MAXmat])
	  malloc(sizeof(double)*MAXmat*NB+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }
  
      if((WORKD3=(double (* restrict)[MAXmat])
	  malloc(sizeof(double)*MAXmat*NB+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((WORK=(SI (*restrict)[MAXmat])
	  malloc(sizeof(SI)*MAXmat*NB+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((WORK2=(UI (*restrict)[MAXmat])
	  malloc(sizeof(UI)*MAXmat*NB+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((polya=(UI * restrict)malloc(sizeof(UI)*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((polyb=(UI * restrict)malloc(sizeof(UI)*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((polyc=(UI * restrict)malloc(sizeof(UI)*(MAXmat+1)+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((ut=(UI * restrict)malloc(sizeof(UI)*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((vt=(UI * restrict)malloc(sizeof(UI)*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((v2t=(UI * restrict)malloc(sizeof(UI)*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

      if((sind=(UI * restrict)malloc(sizeof(UI)*2*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

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

      if((WORKD2=(double * restrict)malloc(sizeof(double)*MAXmat+SPEED))==NULL){
	printf("out of memory.\n");
	CM_flag=1;
      }

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

	if(CM_flag==1){
	  continue;
	}

	if((CM_result[i]=(UI * restrict)malloc(sizeof(UI)*MAXmat))==NULL){
	  printf("Out of memory!\n");
	  fflush(stdout);
	  CM_flag=1;
	  continue;
	}

	CM=CMD[i];
	CMh=(CM-1)/2;
	CMhD[i]=CMh;

	if(NB<MAXmat){

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

	    if(CM_flag==1){
	      continue;
	    }
	
	    memset(A,0,sizeof(double)*MAXmat*MAXmat);

	    for (k = 0; k < MAXmat; k++){
	      for(j=IAP[k];j<IAP[k+1];j++){
		v = mpz_mod_ui (chtmpu, zmatrix[j],CM);
		if (v>CMh)
		  v=v-CM;

		A[k][JA[j]]=v;
	      }
	    }

	    for(j=0;j<NB;j++){
	      for(k=0;k<MAXmat;k++){
		v=(SLL)rand()*(SLL)rand() % CM;
		if (v>CMh)
		  v=v-CM;
		else if (v<-CMh)
		  v=v+CM;

		B[j][k]=v;
	      }
	    }
	
	    A_mul_v_(A,B,C,WORKD3,CM,CMh);
      
	    fgetrf_blas_inner_(MAXmat,MAXmat,MAXmat,B,
			       IPIV,WORK,WORKD,WORKD2,CM,CMh,&INFO);
	
	    if(INFO!=0){
	      printf("F1");
	      fflush(stdout);
	      continue;
	    }
	    
	    dlaswp_(&NB,C,&MAXmat,&IONE,&MAXmat,IPIV,&IONE);
	
	    ftrsmL_recursive_(MAXmat,NB,MAXmat,B,MAXmat,C,CM,CMh);
	    
	    ftrsmU_recursive_(MAXmat,NB,MAXmat,B,MAXmat,C,WORK,CM,CMh);
	    
	    for(k=0;k<NB;k++){
	      for(j=0;j<MAXmat;j++){
		WORK2[k][j]=(SI)C[k][j]+CM;
	      }   
	    }
	    
	    break;
	  }

	  if(l>=UP0){
	    free(CM_result[i]);
	    CM_result[i]=NULL;
	    continue;
	  }

	}
	else{
	  memset(WORK2,0,sizeof(UI)*MAXmat*MAXmat);
	  
	  for (k = 0; k < MAXmat; k++){
	    for(j=IAP[k];j<IAP[k+1];j++){
	      v = mpz_mod_ui (chtmpu, zmatrix[j],CM);
	      if (v>CMh)
		v=v-CM;
	      WORK2[k][JA[j]]=v;
	    }
	  }
	}

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

	  if(CM_flag==1){
	    continue;
	  }

	  if(ber(sind,polya,polyb,polyc,WORK2,ut,vt,v2t,CM)==0){
	    printf("F2");
	    fflush(stdout);
	    continue;
	  }

	  break;
	}

	if (l>=UP0){
	  free(CM_result[i]);
	  CM_result[i]=NULL;
	  continue;
	}

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

      mpz_clear(chtmpu);

      free(A);
      free(B);
      free(C);
      free(WORKD3);
      free(WORK);
      free(WORK2);

      free(polya);
      free(polyb);
      free(polyc);
      free(ut);
      free(vt);
      free(v2t);
      free(sind);
      free(IPIV);
      free(WORKD2);
    }

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

    n=0;
    for(i=0;i<num_CM;i++){
      if(CM_result[i]!=NULL){
	if(n==0){
	  k=i;
	}
	n=n+1;
      }
      else{
	m=m+1;
      }
    }

    if(m>=UP1){
      printf("FAULT!\n");
      fflush(stdout);
      return 0;
    }

    if(n>0){

      if((tD=(UI * restrict)malloc(sizeof(UI)*num_CM))==NULL){
	printf("Out of memory!\n");
	return 0;
      }
      
      if((chtmps=(mpz_t * restrict)malloc(sizeof(mpz_t)*num_CM))==NULL){
	printf("Out of memory!\n");
	return 0;
      }
      
      for(i=0;i<num_CM;i++){
	mpz_init(chtmps[i]);
      }
      
      mpz_set(chtmps[k],chtmp2);
      for(i=k+1;i<num_CM;i++){
	if(CM_result[i]!=NULL){
	  mpz_mul_ui(chtmps[i],chtmps[k],CMD[k]);
	  k=i;
	}
      }
      mpz_mul_ui(chtmp2,chtmps[k],CMD[k]);
      
      for(i=0;i<num_CM;i++){
	if(CM_result[i]!=NULL){
	  s1=mpz_mod_ui(chtmp1,chtmps[i],CMD[i]);
	  tD[i]=invm2_(s1,CMD[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++){
	    
	    if(CM_result[i]!=NULL){
	      
	      tmp=((UI)(CM_result[i][j]+CMD[i]-(UI)mpz_mod_ui(chtmpu,resultZ[j],CMD[i]))) % CMD[i];
	      
	      if(tmp!=0){
		
		tmp=((ULL)tmp*(ULL)tD[i]) % CMD[i];
		
		tmp=(tmp>CMhD[i]) ? tmp-CMD[i] : tmp;
		
		mpz_mul_si(chtmpu,chtmps[i],(long int)(int)tmp);
		
		mpz_add(resultZ[j],resultZ[j],chtmpu);
	      }
	      
	    }
	  }
	}
	
	mpz_clear(chtmpu);
	
      }
      
      if(n==num_CM)
	return MAXmat;
      
      for(j=0;j<num_CM;j++){
	if(CM_result[j]!=NULL){
	  free(CM_result[j]);
	}
      }

    }

    for(i=0;i<num_CM;i++){
      mpz_clear(chtmps[i]);
    }
    free(chtmps);

    free(tD);

    free(CMhD);

    free(CM_result);

    free(CMD);

  }

  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);
  mpz_init(chtmp4);

  MAXthreads = omp_get_max_threads ();

  srand(getpid());

  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;
	}

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

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

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

      IAP[0] = 0;

      NB=CMIN(MAXmat,NB2);

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

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

      break;
    }
  }
  for(k=k+1;k<=MAXmat;k++){
    IAP[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=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;

}