#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "basis.h"

/* 
         Subroutines for adding states to a sorted basis
         and finding a state in a sorted basis.
         Summer 2002 version.        
*/

#define PRINTIT 0 /*print basis size & timing information in makebasis */

/*
  Add ly terms y to a sorted list x of length lx and max length maxlx. 
*/

void addterms(element *x, int *lx, int maxlx, element *y, int ny){
  int i,j,k;
  for(i=0; i<ny; i++){
    for(j=0; j<*lx && x[j]<y[i]; j++);
    if(*lx<maxlx)(*lx)++;
    for(k=*lx-1; k>j; k--)x[k]=x[k-1];
    if(j<maxlx)x[j]=y[i];
  }
  return;
}

/* print a state to stderr without knowing what type it is */

void printanystate(FILE *fp, partype *state, int np){
  int j;
  for(j=0; j<np; j++){ 
    if(state[j]==HOLE)break;
    fprintf(fp," (%i,%i)",state[j]>>HSHIFT,state[j]&KMASK);
  }
}


/* Add state to sorted basis;
   the state is not added if it is a copy. Returns new length of basis.
   The speed of this routine is the limiting factor for
   the generation of large bases. 

   The order for sorting is a little peculular since it mixes
   up different particle numbers.

*/

size_t addstate(partype *basis, partype *state, size_t nst, int np){
  size_t i,l=0,u=nst;
  int length=np*sizeof(partype);
  /*int j;*/
  int order;
  partype *bp;
#if 0   /* debugging print */
   fprintf(stderr,__FILE__":      add state to %i:  ",(int) nst);
  for(j=0;j<np;j++) fprintf(stderr,__FILE__":  %i ",(int) state[j]);
   fprintf(stderr,__FILE__":  \n");  */
#endif
  while(u-l>0){
    i=((u-l)>>1)+l;
    bp=basis+i*np;
    order=order_states(state,bp,np);
    if(order==0)return nst;
    else if(order>0)
      u=i;
    else
      l=i+1; 
  }
/*  printf("values u=%i, l=%i\n",u,l); */
/* 
   The following line is the theoretical limiting time factor
   in the generation of a large basis.
*/
  memmove(basis+(u+1)*np,basis+u*np,(nst-u)*length);
  memcpy(basis+u*np,state,length);

  return nst+1;
}

/* 
   Return position of state in basis (length nst) assuming the
   state is in the basis and the basis is sorted.  

   It is not a good idea to use memcmp() in the following 
   because if partype spans two bytes, memcmp() can assign
   an incorrect order for the two states.    
*/

size_t findnext(partype *basis, partype *state, size_t nst, int np){
  size_t i,l=0,u=1;
  int j;
  int order;
  partype *bp;
  for(;u<nst;){
    bp=basis+u*np;
    order=order_states(state,bp,np);
    if(order==0)return u;
    else if(order<0){
      l=u+1; 
      u<<=1;
    }else
      break;
  }
  if(u>nst)u=nst;
  while(u-l>0){
    i=((u-l)>>1)+l;
    bp=basis+i*np;
    order=order_states(state,bp,np);
    if(order==0)return i;
    else if(order>0)
      u=i;
    else
      l=i+1;
  }
  fprintf(stderr,__FILE__":  findnext():  Can't find state in basis,"
	  " p=%i, exiting\n",np);
  fprintf(stderr,"  State: ");
  printanystate(stderr,state,np);
  fprintf(stderr,"\n");
  fprintf(stderr,"  in basis (first 50): \n");
  for(j=0; j<nst && j<50; j++){
    fprintf(stderr,"  ");
    printanystate(stderr,basis+j*np,np);
    fprintf(stderr,"\n");
  }
  return nst;
}

/*
  give a canonical order to two states.  
  It is not a good idea to use memcmp() in the following 
  because if partype spans two bytes, memcmp() can assign
  an incorrect order for the two states.
*/

int order_states(partype *s1, partype *s2, int np){
  int i;
  for(i=0; i<np; i++){
    if(s1[i]==s2[i])continue;
    if(s1[i]<s2[i])return 1;
    else return -1;
  }
  return 0;
}

/*
  In the following routines, we represent a sorted basis 
  as a binary tree. 
  This allows one to add states to the basis witout
  moving lots of memory around.
*/

/*
  allocate memory for the tree_basis structure.
*/

tree_basis *tree_realloc(tree_basis *tree, size_t length, int np){
  if(np>NPMAX){
    fprintf(stderr,__FILE__":%i  np=%i > NPMAX=%i in tree_realloc \n",
	    __LINE__,np,NPMAX);
    return NULL;
  }

  tree=realloc(tree,length*sizeof(tree_basis));
  if(length>0 && tree==NULL)
    fprintf(stderr,__FILE__":%i  Out of memory in tree_realloc\n",__LINE__);

  return tree;
}

/*
  add a state to the basis.  This routine assumes the tree structure has
  sufficient memory available.  It returns the new size of the
  basis.

  Normally this takes order log(N) operations, for N states.
  If the states are already sorted, then this takes order N operations.
*/

#define LEAF -1  /* the end of a branch */

size_t tree_addstate(tree_basis *tree, partype *state, size_t length, int np){
  int order, k;
  size_t branch=0;
  int check=length;

  while(check--){ /* for safety.  Also, don't do length=0 case */
    order=order_states(state,tree[branch].s,np);
    if(order==0)
      return length;
    if(order<0){       /*  the sign is chosen to match addstate */ 
      if(tree[branch].upper==LEAF){
	tree[branch].upper=length;
	break;
      }
      branch=tree[branch].upper;
    } else {
      if(tree[branch].lower==LEAF){
	tree[branch].lower=length;
	break;
      }
      branch=tree[branch].lower;
    }
  }

#if 0  /* debugging print */
  printf("length=%i,  adding state:  ",length);
  printanystate(stdout,state,np);
  printf("\n");
#endif

  for(k=0; k<np; k++)tree[length].s[k]=state[k];
  tree[length].upper=LEAF;
  tree[length].lower=LEAF;
  return length+1;
}

/*
  turn the tree structure into an ordinary basis.
  it deallocates the ordinary basis and returns a pointer
  to the ordinary basis.
*/

partype *tree_sort(tree_basis *tree, size_t length, int np){
  size_t branch,last_branch=-1;
  int k;
  partype *bp;
  partype *basis=NULL;
  size_t trunk=0;

#if 0  /* debugging print: print basis */
  printf("basis in tree form:\n");
  for(k=0; k<length; k++){
    printf("  %i:  ",k);
    printanystate(stdout,tree[k].s,np);
    printf("  " SIZE_T_PRINT "  " SIZE_T_PRINT "\n",
	   tree[k].lower,tree[k].upper);
  }
#endif

  if(length*np==0)goto end;
  basis=malloc(length*np*sizeof(partype));
  if(basis==NULL){
    fprintf(stderr,__FILE__":%i  Out of memory in tree_sorted\n",__LINE__);
    return basis;
  }
  
  bp=basis;
  while(length>0){
    branch=trunk;
    while(tree[branch].lower!=LEAF){
      last_branch=branch;
      branch=tree[branch].lower;
    }
    length--;  /* safety check */
    for(k=0; k<np; k++)*(bp++)=tree[branch].s[k]; /* copy state to basis */
    if(branch==trunk){
      if(tree[branch].upper==LEAF)break; /* no branches left */
      trunk=tree[branch].upper;          /* move trunk */
    }
    else if(tree[branch].upper==LEAF){
      tree[last_branch].lower=LEAF;          /* remove branch */
    } else 
      tree[last_branch].lower=tree[branch].upper; /* move branch */
  }

  /* check for success */
  if(tree[trunk].lower!=LEAF || tree[trunk].upper!=LEAF){
    fprintf(stderr,__FILE__":%i  tree cutting was unsuccessful\n",__LINE__);
    free(basis);
    basis=NULL;
  }
  
 end:
  free(tree);

  return basis;
}