/************************************************************************

  Copyright
  Alessandro MIRONE
  mirone@esrf.fr

  Copyright 2002  by European Synchrotron Radiation Facility, Grenoble, 
                  France

                               ----------
 
                           All Rights Reserved
 
                               ----------

Permission to use, copy, modify, and distribute this software and its
documentation for any purpose and without fee is hereby granted,
provided that the above copyright notice appear in all copies and that
both that copyright notice and this permission notice appear in
supporting documentation, and that the names of European Synchrotron
Radiation Facility or ESRF or SCISOFT not be used in advertising or 
publicity pertaining to distribution of the software without specific, 
written prior permission.

EUROPEAN SYNCHROTRON RADIATION FACILITY DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL EUROPEAN SYNCHROTRON
RADIATION FACILITY OR ESRF BE LIABLE FOR ANY SPECIAL, INDIRECT OR 
CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, 
DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.

**************************************************************************/

#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;}
  inline Hxx_FermionicState *get_basis(int n=0) { return this->basis+n;}


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  Read_add_forhops( char * nome_file, Hxx_FLOAT coeff, int binary,
			  int * cc, double *hopfacts) ;


  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
