// ************************************
// Alessandro MIRONE 
// Peter      KRUEGER
//
// September 2002
// *************************************

#ifndef __HILBERTXX__
#define __HILBERTXX__





// #define first_bit_mask 00000000000010000000000

static  unsigned long   first_bit_mask_ =  010000000000 ;
static  unsigned long   first_bit_mask  =  2* first_bit_mask_ ;


#define last_bit_mask  000000000001

class  Hxx_TypeForBits {
// Class to represent second quantization states.
// 32 bits ( on intel ) was  enough at the beginning.
// The second quatization states were represented  with longs :
//     #define Hxx_TypeForBits long
// Now  we use a custom class containing a
// list of longs. Their number is N_Hxx_TypeForBits :
#define N_Hxx_TypeForBits 2

public:
  inline friend  Hxx_TypeForBits   operator >>(Hxx_TypeForBits  &a, int i   ); 
  inline friend  Hxx_TypeForBits   operator <<(Hxx_TypeForBits  &a, int i   ); 
  inline friend  Hxx_TypeForBits   operator & (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b); 
  inline friend  Hxx_TypeForBits   operator | (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b); 
  inline friend  Hxx_TypeForBits   operator ^ (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b); 
  inline friend  Hxx_TypeForBits   operator ~ (Hxx_TypeForBits  &a) ;

  inline friend  int   operator == (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b); 
  inline friend  int   operator != (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b); 
  inline friend  int   operator < (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b); 

  int countOne();

  void return_string(char*);

  inline void setTo0() {
    for (int i=0; i< N_Hxx_TypeForBits ; i++) {
      x[i]=0;
    }
  };

  void print_values() ;

  inline void setToOne(){
    for (int i=0; i< N_Hxx_TypeForBits ; i++) {
      x[i]=0;
    }
    x[N_Hxx_TypeForBits-1]=last_bit_mask  ;
  } 
  
  inline void setToValue(int v) { x[N_Hxx_TypeForBits-1]=v; } ;

  inline int isZero()  {
    for (int i=0; i< N_Hxx_TypeForBits ; i++) {
      if(x[i]) return 0;
    }
    return 1;
  };


  inline unsigned long   lastBit()  {  return x[N_Hxx_TypeForBits-1]&last_bit_mask ;  } ;
  inline unsigned long   firstBit() { return x[0]&first_bit_mask ; } ;
  unsigned long    x[N_Hxx_TypeForBits];
}  ;




inline Hxx_TypeForBits   operator >>(Hxx_TypeForBits  &a, int ishift   ) {
  Hxx_TypeForBits res = a;
  unsigned long    residual;
  unsigned long    new_residual;

  for(int ish =0; ish<ishift; ish++) {
    residual=0;

    for(int i=0; i< N_Hxx_TypeForBits ; i++) {

      new_residual = res.x[i] & last_bit_mask ;

      res.x[i] = res.x[i] >> 1;
      if(residual) {

	res.x[i] = res.x[i] |  first_bit_mask;
      }
      residual = new_residual;
    }
  }
  return res;
}


inline Hxx_TypeForBits   operator << (Hxx_TypeForBits  &a, int ishift   ) {

  Hxx_TypeForBits res = a;
  unsigned long    residual;
  unsigned long    new_residual;
  for(int ish =0; ish<ishift; ish++) {
    residual=0;

    for(int i= N_Hxx_TypeForBits-1; i>=0 ; i--) {
       
      new_residual = res.x[i] & first_bit_mask ;
      res.x[i] = res.x[i] << 1;
      if(residual) {
	res.x[i] = res.x[i] | last_bit_mask ;
      }
      residual = new_residual;
    }
  }
  return res;
}


inline   Hxx_TypeForBits   operator ~ (Hxx_TypeForBits  &a) {
  Hxx_TypeForBits res ;
  for(int i=0; i< N_Hxx_TypeForBits ; i++) {
    res.x[i]= ~  a.x[i];
  }
  return res;
} 



inline   Hxx_TypeForBits   operator & (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b) {
  Hxx_TypeForBits res ;
  for(int i=0; i< N_Hxx_TypeForBits ; i++) {
    res.x[i]= a.x[i] & b.x[i];
  }
  return res;
} 

inline   Hxx_TypeForBits   operator | (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b) {
  Hxx_TypeForBits res ;
  for(int i=0; i< N_Hxx_TypeForBits ; i++) {
    res.x[i]= a.x[i] |  b.x[i];
  }
  return res;
} 

inline   Hxx_TypeForBits   operator ^ (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b) {
  Hxx_TypeForBits res ;
  for(int i=0; i< N_Hxx_TypeForBits ; i++) {
    res.x[i]= a.x[i] ^  b.x[i];
  }
  return res;
} 


inline  int   operator == (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b) {
  for(int i=0; i< N_Hxx_TypeForBits ; i++) {
    if(a.x[i] !=  b.x[i]) return 0;
  }
  return 1;
} 


inline  int   operator != (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b) {
  for(int i=0; i< N_Hxx_TypeForBits ; i++) {
    if(a.x[i] !=  b.x[i]) return 1;
  }
  return 0;
} 


inline  int   operator < (Hxx_TypeForBits  &a,  Hxx_TypeForBits  &b) {
  for(int i=0; i< N_Hxx_TypeForBits ; i++) {
    
    if(a.x[i] <  b.x[i])  return 1;
    else if(a.x[i] >  b.x[i]) return 0;
  }
  return 0;
} 

class Hxx_Filter {

public:
  void addCondition (Hxx_TypeForBits &mask, char operation, double factor);
  Hxx_Filter (int);
  ~Hxx_Filter ();
  int filtre(Hxx_TypeForBits &state);

  struct opFilter {
    Hxx_TypeForBits mask;
    char operation;
    double factor;
  } opFilter ;

  double min;
  double max;


private:

  struct opFilter* filtreOperation;
  int Hxx_Filter::validate(double res);
  int nop;
};

class Hxx_Utilities {
public:
  static int * SCRAMBLE_FOR_BITS ;
  static int * create_scrambler(int nstates);

};
// that because gnu-cc does not not like type casting from void *
#define Hxx_Object char


#define Hxx_FLOAT double

class Hxx_Chunker;

/**
 * provides a chained list for memory buffers
 */
class Hxx_Chunck {
public:
  char * ptr;
  Hxx_Chunck  * next;
};



/**
 * @short[ Hxx_Chunker provides  memory ]
 * Allocation is faster if done by bunches. Hxx_Chunck have
 * the get_space() that return a pointer to a bloc of memory whose
 * size is given by the stride parameter.
 * Allocation is not done on a one by one basis, but allocating
 * when necessary a buffer having size stride*chunck_size bytes
 */

class Hxx_Chunker {
public:

  /** 
   * Create a new object and set the stride and chunk_size parameters 
   */
  Hxx_Chunker(int stride, int chunck_size ) ;

  /** 
   * Give a pointer to N=stride free bytes 
   */
  inline Hxx_Object  * get_space();  

  /** 
   * When the memory given by the object is no more necessary, it can be freed by this function 
   */ 
  void free_space();

  inline ~Hxx_Chunker() {
    this->free_space();
  };


private:
  int stride;
  int chunck_size;
  int chunck_pos ; 
  Hxx_Chunck * top_chunk;
  Hxx_Chunck * actual_chunk ; 
};


class Hxx_1p_FermionicOperator_;

/**
 * @short[ Objects of this class represent a second quantisation state ]
 * Objects of this class represent a second quantisation state
 * as a sequence of 1 and 0, like for example
 *     100101011
 * Which are the occupancies. Together with the occupancies
 * the signes to extract a given fermion are given as the integral modulus 2 
 * of the occupancies
 *     011100110
 */
class Hxx_FermionicState {
public:
  
  Hxx_FermionicState(){};
  Hxx_FermionicState(Hxx_TypeForBits occupancies);
  void initialise(Hxx_TypeForBits occupancies);

  static  inline int  are_equal(Hxx_FermionicState *a, Hxx_FermionicState *b ) ;
  static  inline int  less(Hxx_FermionicState *a, Hxx_FermionicState *b  ) ;
  friend class Hxx_1p_FermionicOperator_;


  // just for testing, never use this function ......
  inline Hxx_TypeForBits& get_value() { return occupancies ; }
  void print_values(int mask=0);
  int n_particles();
  int n_particles(int no, int *ops);
  int tag;
  Hxx_TypeForBits occupancies;
private:

  Hxx_TypeForBits signes;
};

class Hxx_1p_FermionicOperator_  {
public:
    
  
  inline Hxx_1p_FermionicOperator_(int pos=0) { 
this->initialize(pos); 
};
  void initialize( int pos);
  
  inline void Creation( Hxx_FermionicState * original , Hxx_FermionicState * result,
			int & sign ) ;
  
  inline void Destruction( Hxx_FermionicState * original , Hxx_FermionicState * result,
			   int & sign ) ;

  inline Hxx_TypeForBits & get_1p_state() { return this->c_d_mask; }
  
private:
  
  inline void operation( Hxx_FermionicState * original , Hxx_FermionicState * result,
			 int & sign) ;
  Hxx_TypeForBits c_d_mask ;
  Hxx_TypeForBits sign_mask;
  
};

class Hxx_1p_FermionicOperators  {
public:
  
  static Hxx_1p_FermionicOperator_ * operators ;
  
  static Hxx_1p_FermionicOperator_ * Initialise();
  
};

/**
 * @short[ Hxx_N_counter counts occupancy ]
 * Initialisation is done with an array of ints that refers to  Hxx_1p_FermionicOperators::operators
 * and a lenght
 */

class Hxx_N_counter  {

public:

  inline Hxx_N_counter( int n, int *ops, int reference_n=0) {
    this->n = n ;
    this->reference_n=reference_n;
    this->ops = new int [n];
    memcpy(this->ops, ops, n *sizeof(int) );
  };
  inline ~Hxx_N_counter() {
    delete[] this->ops;
  }


  int N(Hxx_FermionicState* orig);


  int n;
  int *ops;
  int reference_n;
};


typedef int  (Hxx_Generic_Comparator) (Hxx_Object   *a, Hxx_Object *b, Hxx_Object * user_data );


class Hxx_Tree ;



class Hxx_Normal_Operators_Collection {
public:

  Hxx_Normal_Operators_Collection( int Ncreation, int Ndestruction,
				   int chunck_size = 10000);

  void add_product(int * creators, int *destructors, Hxx_FLOAT coeff );
  void reorder();
  
  inline int get_n_items() { return this->n_items; }
  
  Hxx_FLOAT operate( int i, Hxx_FermionicState* orig, Hxx_FermionicState* result, int & next);
  void print();
  int Ncreation    ;
  int Ndestruction ; 

private:

  char * get_item_space();
  void clean_space();
  static  int  less(int  *a, int  *b , int *user_data ) ;
  void add_memorized(int *newTuple, int count);
  void set_next_indexes();


  int user_data[2] ;  //  Ncreation, Ndestruction
  

  int n_items;

  Hxx_Tree  * orderer ;
  
  Hxx_Chunker * chunker;

  int *indexes;
  Hxx_FLOAT * coeffs;
  int *next_indexes;

  int size_of_teeth;
  int op_order;

  double average_depth;
  double average_depth2;
  int size_of_creators;
  int size_of_destructors; 
};

// ---------------------------------------------------------------------
class Hxx_Tree_Node {
public:
  
  Hxx_Tree_Node *left;
  Hxx_Tree_Node *right;
  Hxx_Tree_Node *previous;
  
  Hxx_Object * Leaf;
};

class Hxx_Tree : public Hxx_Tree_Node  {
public:

  Hxx_Tree(Hxx_Generic_Comparator * Less , Hxx_Object * user_data=NULL , int chunck_size=10000) ;
  
  Hxx_Tree_Node * TopNode;
  Hxx_Tree_Node * ActualNode;
  
  Hxx_Object *  Add_Leaf(Hxx_Object * leaf, int &depth);
  
  Hxx_Object * getLowest(int &depth);
  Hxx_Object * getLowest_and_remove(int &depth);
  Hxx_Object * Again_getLowest_and_remove(int &depth);

  void Free_Chunkes();

  inline int get_n_items() { return n_items; }

  inline ~Hxx_Tree() {
    this->Free_Chunkes();
  };
  

private:

  Hxx_Generic_Comparator * Less   ;

  Hxx_Chunker * chunker;

  Hxx_Object * user_data;

  int n_items;
};

// -------------------------------------------------------------------------------------------
typedef int (CounterForHxxBasisFilter) (Hxx_FermionicState *);    
typedef int (*Callback) (Hxx_FermionicState *, void *);

class Hxx_basis {
public:
  Hxx_basis(int initialdimension, int chunck_size=10000);
  void add_state(Hxx_FermionicState * state);

int operate_andputin_Tree( Hxx_Normal_Operators_Collection * collection , int Ncolls=1, Hxx_Filter* filtre=0 );

  void reconstruct_basis();
  void set_filter(Callback filtre, void *clientdata);

  void reconstruct_basis(CounterForHxxBasisFilter *);
  void reconstruct_basis(Hxx_Filter* Hxfiltre);	



  inline int get_dimension() { return this->nstates; }

  int find_state(Hxx_FermionicState * state  );

  inline Hxx_FermionicState *get_basis() { return this->basis;}

private:
  int nstates;
  int Ncreated;
  Hxx_FermionicState *basis ;
  Callback  filtre;
  void *clientdata;
  Hxx_Tree  * orderer ;
  Hxx_Chunker * states_chunker;
  double average_depth;
  double average_depth2;
};

// -----------------------------------------------------------------------------------------

class Hxx_transition {
public:
  
  inline Hxx_transition(int from , int to, Hxx_FLOAT coeff) {
    this->initialise(from,to,coeff);
  }
  
  inline void initialise(int from , int to, Hxx_FLOAT coeff) {   
    this->from=from;
    this->to=to;
    this->coeff=coeff;
  }
  
  inline void add_coeff(Hxx_FLOAT coeff) {
    this->coeff +=coeff;
  }
  
  static  int less(Hxx_transition *a, Hxx_transition *b) ;

  inline int get_from() { return this->from;}
  inline int get_to() { return this->to;}
  inline Hxx_FLOAT  get_coeff() { return this->coeff;}
  

private:
  int from;
  int to;
  Hxx_FLOAT coeff;
};

inline int Hxx_transition::less(Hxx_transition *a, Hxx_transition *b) {
  if(a->from < b->from) {
    return 1;
  } 
  if(a->from > b->from) {
    return 0;
  } 
  
  if(a->to < b->to) {
    return 1;
  } 
  if(a->to > b->to) {
    return 0;
  } 
  return 0;
}

class Hxx_TransitionMatrix {
public:

  Hxx_TransitionMatrix();
  ~Hxx_TransitionMatrix();

  void add_contribution(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			Hxx_Normal_Operators_Collection  * Norm_Op, Hxx_FLOAT coeff=1.0, int complete_basis=1);

  void add_contribution_debug(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			Hxx_Normal_Operators_Collection  * Norm_Op, Hxx_FLOAT coeff=1.0, int complete_basis=1);

  void add_contribution(Hxx_basis * a_basis, Hxx_basis *b_basis ,
			int n_counters, Hxx_N_counter  ** diag_pro, Hxx_FLOAT  coeff);

  void Write_and_Clean(char *nome_file, int binary=1);


  static void  Hxx_TransitionMatrix::Read_arrays( char * nome_file, int &Nels , int *&from, int *&to, double *&coeffs, int binary) ;

     
  static void  Hxx_TransitionMatrix::Read_4arrays( char * nome_file, int &Nels , int *&from, 
						   int *&to, double *&coeffs,int * & qcol) ;



  void Read_add(char *nome_file, Hxx_FLOAT, int binary=1);


  void BuildAux();
  void AddFrom(Hxx_TransitionMatrix * term, Hxx_FLOAT fact, int dosquare=0);


  void create_3Array_and_Clean(int *& from, int *& to, Hxx_FLOAT *& coeff, int &ntrans);

  void free_space();


private:

  Hxx_Tree *Transitions_Orderer;
  Hxx_Chunker * chunker;
  int n_basis_from;
  int n_basis_to;

  int   aux_total_transitions;
  int * aux_i_from;
  int * aux_i_to;
  Hxx_FLOAT * aux_coeff;
};







// ############################################################################################
//                                      END OF INTERFACE
// ############################################################################################


inline Hxx_FermionicState::Hxx_FermionicState(Hxx_TypeForBits occupancies) {
  this->initialise(occupancies);
}

inline void Hxx_FermionicState::initialise(Hxx_TypeForBits occupancies) {
  this->occupancies=occupancies;
  this->signes.setTo0();
  int n = 8*sizeof(Hxx_TypeForBits);
  Hxx_TypeForBits x;

  x = this->occupancies;

  x = x>>1;

  for(int k=0; k<n;k++) {
    this->signes = this->signes ^ x;
    x = x>>1;
  }
  tag=0;

}



inline int  Hxx_FermionicState::are_equal(Hxx_FermionicState *a, Hxx_FermionicState *b ) {
  return a->occupancies == b->occupancies;
}

inline int  Hxx_FermionicState::less(Hxx_FermionicState *a, Hxx_FermionicState *b ) {
  return a->occupancies  < b->occupancies;
}


inline void Hxx_1p_FermionicOperator_::Creation( Hxx_FermionicState * original , Hxx_FermionicState * result,
						 int & sign ) {    
  sign = (  (original->occupancies        &  this->c_d_mask).isZero()) ? sign : 0  ;
  if(sign==0) return;
  this->operation(original,result,sign);
}

inline void Hxx_1p_FermionicOperator_::Destruction( Hxx_FermionicState * original , Hxx_FermionicState * result,
						    int & sign ) {    
  sign = (  (original->occupancies        &  this->c_d_mask).isZero()) ? 0 : sign  ;
  if(sign==0) return;
  this->operation(original,result,sign);
}


inline void Hxx_1p_FermionicOperator_::operation( Hxx_FermionicState * original , Hxx_FermionicState * result,
						  int & sign) {
  sign  = (  (original->signes        &  this->c_d_mask).isZero()) ? sign : -sign ;
  result->occupancies = original->occupancies  ^  (this->c_d_mask) ;
  result->signes      = original->signes       ^  (this->sign_mask) ;
}





inline char *  Hxx_Chunker::get_space() {
  if(this->chunck_pos < this->chunck_size-1) {
    this->chunck_pos++;
    return   this->actual_chunk->ptr + (this->chunck_pos  -  1) * this->stride    ; 
  } else {
    this->actual_chunk->next  = new Hxx_Chunck;
    this->actual_chunk  =this->actual_chunk->next  ;
    this->chunck_pos=0;
    this->actual_chunk->next  = NULL ; 
    this->actual_chunk->ptr = new char [(this->chunck_size)*(this->stride)];
    return this->get_space();
  }
};  

void Inizializza_Hilbertxx();


#endif
