#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fortran.h"
#include "basis.h"

int MAIN(){
  int i,j,k,np=6,kt,multi,o,charge;
  int hlist[HINDEX];
  reduced_basis *basis;
  partype p;
  latype g;
  enum state_type type=TYPE_SOURCE;

#if 0 /* loop thingy */
  for(i=0; i<7; i++){
    int shift,j;
    printf("starting length=%i\n",i);
    for(shift=1; i>=shift; shift*=2);
    while((shift>>=1) >= 1){
      for(j=shift-1; j<i; j+=2*shift)
	printf(" %i ",j);
      printf("\n");
    }
  }
#endif

#if 1  /* Generate basis */

  printf("Generate %s basis using makereduced.\n",
	 state_type_string(type));
  printf("Input ");
  if(type==TYPE_GLUE || type==TYPE_SOURCE)printf("ht[%i], ",HINDEX);
  printf("np, 2*K, charge, o, multi:\n");
  if(type==TYPE_GLUE || type==TYPE_SOURCE)
    for(i=0; i<HINDEX; i++)scanf("%i",hlist+i);
  scanf("%i %i %i %i %i",&np,&kt,&charge,&o,&multi);
  if((i=makereduced(&basis,type,hlist,np,kt,charge,o,subgroup[multi])))
    printf("Error code:  %i\n",i);
  printf("This includes the number of distinct states under\n"
	 " the given symmetry, including zero norm states.\n");
  
  if(basis->length<200){
    printf("(helicity,K) shown for each particle\n");
    for(i=0; i<basis->length; i++){
      for(j=0; j<np && basis->s[np*i+j]!=HOLE; j++);
      k=innerproduct(basis->s+np*i,basis->s+np*i,type,j,charge,o,multi);
      printf(" ");
      printstate(basis->s+np*i,type,np);
#if 1  /* print out inner product */
      printf(" %i",k);
#endif
      printf("\n");
    } 
  }
  printf("length = %i \n",basis->length); 
  free_reduced();

#if 0  /* Test action of the lattice operator on a state */

  printf("Transforms by statelat.  Input state number:  ");
  scanf("%i",&i);
  for(g=0; g<ORDER_GROUP; g++){
    partype out[NPMAX];
    printf("    group element %i:  ",g);
    statelat(out,basis->s+np*i,type,np,g);
    printstate(out,type,np);
    printf("\n");
  }
#endif

#elif 0  /* test Jacobi weights routine */

  doublereal xxx[4]={0.1,0.1,0.3,0.3},yyy[4]={0,0.2,0.1,0.2};
  integer kind=5,n=10,kpts=0;
  doublereal alpha=0.5,ends[2]={0.0,1.0},b[10],tt[10],w[10];
  printf("Try out bessel function J_n(x).\n");
  for(j=0;j<4;j++)
    printf("%g %g %.16f\n",xxx[j],yyy[j],MYBESSELJ(xxx+j,yyy+j));
  printf("Try out Jacobi weights, alpha=%f:\n",alpha);
  GAUSSQ(&kind,&n,&alpha,&alpha,&kpts,ends,b,tt,w);
  for(j=0, b[0]=0.0; j<n; j++){
    b[0]+=w[j];
    printf("  %.16f %.16f\n",0.5*(1.0+tt[j]),w[j]/pow(2.0,2.0*alpha+1.0));
  }
  printf("sum=%f\n",b[0]/pow(2.0,2.0*alpha+1.0) );

#elif 0

#define TOP 5
  element cc[TOP*TOP],beta=0.0;
  void inner22(element *, int top, element beta);
  printf("Test instant22\n");
  inner22(cc,TOP,beta);
#if 1
  for(i=0; i<TOP; i++)
    for(j=0; j<TOP; j++)
      printf(" m=%i n=%i %.16f \n",i,j,cc[i*TOP+j]);
#endif
  exit(1);

#elif 0

  element inner13(int,int,int,int,element);
  element b=0.1,pre,in,pi=3.1415926535897932385;
  printf("Test integrals for 1->3 particles inner13, beta=%f\n",b);
#if 1
  for(k=1; k<6; k++)
    for(i=0; i<k; i++)
      for(j=0; i+j<k && j<=i; j++){
	in=inner13(1,k,i,j,b);
	pre=1.0/(4.0*pi*sqrt((0.5+(element) k)*in));
	printf(" K=%i/2 k1=%i/2 k2=%i/2  %.12f %.12f %.12f\n",
	       2*k+1,2*i+1,2*j+1,in,
	       pre*inner13(2,k,i,j,b),
	       0.5*pre*inner13(3,k,i,j,b));
      }
#endif
  exit(1);

#elif 0

  printf("Try out permutations and factorial\n");
  for(j=0;j<4;j++)
    printf("%i! = %i\n",j,factorial(j));
  for(j=1;j<5;j++){
    printf("permutations of %i obects:\n",j);
    for(i=0;i<factorial(j);i++){
      printf("  ");
      for(k=0;k<j;k++)printf("%i ",perms(j)[i*j+k]);
      printf("\n");
    } 
  }  

#elif 0

  printf("Test isinvariantvector(...) for all multiplets.\n");
  printf("group of order %i (helicity,transform)\n",
	 (1<<HINDEX)*factorial(HINDEX));
  for(i=0; i<HINDEX; i++)hlist[i]=0;
  printf("test vector: (");
  for(i=0; i<HINDEX; i++)printf(" %i",hlist[i]); printf(")\n");
  for(i=0; i<NMULTIPLETS; i++)
    if(isinvariantvector(hlist,subgroup(i)))
      printf(" %i",i); printf("\n");
  hlist[0]=5; printf("test vector: (");
  for(i=0; i<HINDEX; i++)printf(" %i",hlist[i]); printf(")\n");
  for(i=0; i<NMULTIPLETS; i++)
    if(isinvariantvector(hlist,subgroup(i)))
      printf(" %i",i); printf("\n");
  if(HINDEX>1){
    hlist[1]=5; printf("test vector: (");
    for(i=0; i<HINDEX; i++)printf(" %i",hlist[i]); printf(")\n");
    for(i=0; i<NMULTIPLETS; i++)
      if(isinvariantvector(hlist,subgroup(i)))
	printf(" %i",i); printf("\n");
    hlist[1]=-5; printf("test vector: (");
    for(i=0; i<HINDEX; i++)printf(" %i",hlist[i]); printf(")\n");
    for(i=0; i<NMULTIPLETS; i++)
      if(isinvariantvector(hlist,subgroup(i)))
	printf(" %i",i); printf("\n");
  } 

#elif 0

  integer nsects=3,th=2;
  /* {{1,2,3},{2,4,5}} */
  doublereal w[100]={1.,2.,3.,8.,4.,5.};
  printf("Test pseudoinverse rountine.\n");
  PSEUDOINVERSE(w,&th,&nsects);
  for(i=0; i<th; i++){
    for(j=0; j<nsects; j++)printf(" %f",w[i*nsects+j]);
    printf("\n");
  }

#elif 0 /* test wellorder and addstate */

  partype basis[] = {4,211,88,33,4,5, 111,2,3,4,6,5};
  partype state[] = {33,22,1,4,5,6};
  for(j=0;j<3;j++){
    for(i=0;i<6;i++) printf("%i ",(int) basis[j*np+i]);
    printf("\n");
  } 
  for(i=0;i<2;i++)wellorder(&basis[i*np],6);
  wellorder(state,6);
  printf("nst after= %i \n",addstate(basis,state,2,6));
  for(j=0;j<3;j++){
    for(i=0;i<6;i++) printf("%i ",(int) basis[j*np+i]);
    printf("\n");
  } 

#elif 0 /* Test action of group on a helicity index */
 
  printf("Action of group on a helicity index\n");
  for(g=0; g<ORDER_GROUP; g++){
    /*printf("g=%i  ",g);*/
    for(i=0; i<2*HINDEX; i++){
      printf("%c%i",i==0?'{':',',gluelat(i<<HSHIFT,g)>>HSHIFT);
    }
    printf("},\n");
  }  


#endif


#if 1==0

 printf("try out resthk; to make a basis, input"
               " ht[%i], np, kt, cyclic, o, multi:\n",HINDEX);
  for(i=0;i<HINDEX;i++)scanf("%i",hlist+i);
  scanf("%i %i %i %i %i",&np,&kt,&cyclic,&o,&multi);
  s=mahk(hlist,np,kt,cyclic,o,subgroup(multi));
  if(s->length<200)
    for(i=0;i< s->length;i++){
      printf(" ");
      for(j=0;j<np;j++) 
	printf("%3i ",(int) (((s->locate)[np*i+j]&~(~0<<HSHIFT))<<1)+1);
      printf(" ");
      for(j=0; j<np; j++) printf("%3i ",(int) (s->locate)[np*i+j]>>HSHIFT);
      printf("\n");
    }    
  printf("loc, length = %p %i \n",s->locate,s->length);
  printf("Input which state to apply resthk to:  ");
  scanf("%i",&k);
  printf("\n");

  /*  Make associated basis and get state from basis*/
  makebasis(&x,hlist,np,kt,cyclic,o,multi);
  t=resthk(&x,np,s->locate+np*k,np);
  printf("%i-particle states\n",np); 
  if(t->length<200)
  for(i=0;i< t->length;i++){
    printf(" ");
    for(j=0;j<np;j++) 
      printf("%3i ",(int) (((t->locate)[np*i+j]&~(~0<<HSHIFT))<<1)+1);
    printf(" ");
    for(j=0; j<np; j++) printf("%3i ",(int) (t->locate)[np*i+j]>>HSHIFT);
    printf("\n");
  }  
  printf("length = %i \n",t->length); 
  t=resthk(&x,np+2,s->locate+np*k,np); 
  printf("%i-particle states\n",np+2);
  if(t->length<200)
  for(i=0;i< t->length;i++){
    printf(" ");
    for(j=0;j<np+2;j++) 
      printf("%3i ",(int) ((t->locate[(np+2)*i+j]&~(~0<<HSHIFT))<<1)+1);
    printf(" ");
    for(j=0; j<np+2; j++) printf("%3i ",(int) (t->locate[(np+2)*i+j]>>HSHIFT));
    printf("\n");
  }                                
  printf("length = %i \n",t->length); 

  if(s->length<50){
    printf("helicity map:\n");
    for(i=0; i < 2*HINDEX; i++){
      p=i<<HSHIFT;
      gluedecode(p,hlist);
      printf("  %i: ",i);
      for(j=0; j<HINDEX; j++) printf("%i ",hlist[j]);
      printf("\n");
    }
  }


#endif


  return 0;
}