Partial cleanup of header files
This commit is contained in:
		
							parent
							
								
									83cad3587e
								
							
						
					
					
						commit
						bdf309d78c
					
				| 
						 | 
				
			
			@ -17,7 +17,7 @@ Purpose:  Core header file, includes all descendents.
 | 
			
		|||
 | 
			
		||||
// This core header file includes all the others
 | 
			
		||||
 | 
			
		||||
const char ABACUS_VERSION[20] = "ABACUS_0a";
 | 
			
		||||
const char ABACUS_VERSION[20] = "1.0.0";
 | 
			
		||||
 | 
			
		||||
// Standard includes
 | 
			
		||||
#include <cmath>
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -21,7 +21,7 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
  // First, some global constants...
 | 
			
		||||
 | 
			
		||||
  const long long int ID_UPPER_LIMIT = 10000000LL;  // max size of vectors we can define without seg fault
 | 
			
		||||
  const long long int ID_UPPER_LIMIT = 10000000LL;  // max size we can define without seg fault
 | 
			
		||||
  const int INTERVALS_SIZE = 100000;                // size of Scan_Intervals arrays
 | 
			
		||||
  const int NBASESMAX = 1000; // max number of bases kept
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			@ -30,7 +30,7 @@ namespace ABACUS {
 | 
			
		|||
  // Cutoffs on particle numbers
 | 
			
		||||
  const int MAXSTRINGS = 20;      // maximal number of particle types we allow in bases
 | 
			
		||||
 | 
			
		||||
  const int NEXC_MAX_HEIS = 16; // maximal number of excitations (string binding/unbinding, particle-hole) considered
 | 
			
		||||
  const int NEXC_MAX_HEIS = 16; // max nr of excitations (string binding/unbinding, particle-hole) considered
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
  //***********************************************************************
 | 
			
		||||
| 
						 | 
				
			
			@ -66,7 +66,8 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
  //****************************************************************************
 | 
			
		||||
 | 
			
		||||
  // Objects in class Heis_Base are a checked vector containing the number of rapidities of allowable types for a given state
 | 
			
		||||
  // Objects in class Heis_Base are a checked vector containing the number of rapidities
 | 
			
		||||
  // of allowable types for a given state
 | 
			
		||||
 | 
			
		||||
  class Heis_Base {
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -116,7 +116,8 @@ namespace ABACUS {
 | 
			
		|||
	if (i % 2 && bdry[i] < xmax_ref) xmax_reg++;
 | 
			
		||||
      }
 | 
			
		||||
 | 
			
		||||
      Vect<T> new_bdry(bdry.size() + 2 * (((xmin_reg + 1) % 2 && (xmax_reg + 1) % 2) - (xmax_reg - xmin_reg)/2));
 | 
			
		||||
      Vect<T> new_bdry(bdry.size()
 | 
			
		||||
		       + 2 * (((xmin_reg + 1) % 2 && (xmax_reg + 1) % 2) - (xmax_reg - xmin_reg)/2));
 | 
			
		||||
 | 
			
		||||
      int ishift = 0;
 | 
			
		||||
      for (int i = 0; i <= xmin_reg; ++i) new_bdry[i] = bdry[i];
 | 
			
		||||
| 
						 | 
				
			
			@ -201,18 +202,40 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
  // ******************************** Recursive integration functions ******************************
 | 
			
		||||
 | 
			
		||||
  DP Integrate_Riemann (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax, int Npts);
 | 
			
		||||
  DP Integrate_Riemann_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
 | 
			
		||||
  DP Integrate_Riemann (DP (*function) (Vect_DP),
 | 
			
		||||
			Vect_DP& args, int arg_to_integ,
 | 
			
		||||
			DP xmin, DP xmax,
 | 
			
		||||
			int Npts);
 | 
			
		||||
 | 
			
		||||
  DP Integrate_Riemann_using_table (DP (*function) (Vect_DP, I_table),
 | 
			
		||||
				    Vect_DP& args, int arg_to_integ,
 | 
			
		||||
				    I_table Itable,
 | 
			
		||||
				    DP xmin, DP xmax, int Npts);
 | 
			
		||||
 | 
			
		||||
  DP Integrate_rec (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax, DP req_prec, int max_rec_level);
 | 
			
		||||
  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
 | 
			
		||||
				DP xmin, DP xmax, DP req_prec, int max_rec_level);
 | 
			
		||||
  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
 | 
			
		||||
				DP xmin, DP xmax, DP req_prec, int max_rec_level, std::ofstream& outfile);
 | 
			
		||||
  DP Integrate_rec_using_table_and_file (DP (*function) (Vect_DP, I_table, std::ofstream&), Vect_DP& args,
 | 
			
		||||
					 int arg_to_integ, I_table Itable,
 | 
			
		||||
					 DP xmin, DP xmax, DP req_prec, int max_rec_level, std::ofstream& outfile);
 | 
			
		||||
  DP Integrate_rec (DP (*function) (Vect_DP),
 | 
			
		||||
		    Vect_DP& args, int arg_to_integ,
 | 
			
		||||
		    DP xmin, DP xmax,
 | 
			
		||||
		    DP req_prec, int max_rec_level);
 | 
			
		||||
 | 
			
		||||
  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table),
 | 
			
		||||
				Vect_DP& args, int arg_to_integ,
 | 
			
		||||
				I_table Itable,
 | 
			
		||||
				DP xmin, DP xmax,
 | 
			
		||||
				DP req_prec, int max_rec_level);
 | 
			
		||||
 | 
			
		||||
  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table),
 | 
			
		||||
				Vect_DP& args, int arg_to_integ,
 | 
			
		||||
				I_table Itable,
 | 
			
		||||
				DP xmin, DP xmax,
 | 
			
		||||
				DP req_prec, int max_rec_level,
 | 
			
		||||
				std::ofstream& outfile);
 | 
			
		||||
 | 
			
		||||
  DP Integrate_rec_using_table_and_file (DP (*function) (Vect_DP, I_table, std::ofstream&),
 | 
			
		||||
					 Vect_DP& args, int arg_to_integ,
 | 
			
		||||
					 I_table Itable,
 | 
			
		||||
					 DP xmin, DP xmax,
 | 
			
		||||
					 DP req_prec, int max_rec_level,
 | 
			
		||||
					 std::ofstream& outfile);
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			@ -248,34 +271,70 @@ namespace ABACUS {
 | 
			
		|||
    DP xmax;
 | 
			
		||||
 | 
			
		||||
  public:
 | 
			
		||||
    Integral_data (DP (*function_ref) (Vect_DP), Vect_DP& args, int arg_to_integ_ref, DP xmin_ref, DP xmax_ref);
 | 
			
		||||
    Integral_data (DP (*function_ref) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ_ref,
 | 
			
		||||
		   I_table Itable, DP xmin_ref, DP xmax_ref);
 | 
			
		||||
    Integral_data (DP (*function_ref) (Vect_DP, Integral_table), Vect_DP& args, int arg_to_integ_ref,
 | 
			
		||||
		   Integral_table Itable, DP xmin_ref, DP xmax_ref);
 | 
			
		||||
    Integral_data (DP (*function_ref) (Vect_DP),
 | 
			
		||||
		   Vect_DP& args, int arg_to_integ_ref,
 | 
			
		||||
		   DP xmin_ref, DP xmax_ref);
 | 
			
		||||
 | 
			
		||||
    Integral_data (DP (*function_ref) (Vect_DP, I_table),
 | 
			
		||||
		   Vect_DP& args, int arg_to_integ_ref,
 | 
			
		||||
		   I_table Itable,
 | 
			
		||||
		   DP xmin_ref, DP xmax_ref);
 | 
			
		||||
 | 
			
		||||
    Integral_data (DP (*function_ref) (Vect_DP, Integral_table),
 | 
			
		||||
		   Vect_DP& args, int arg_to_integ_ref,
 | 
			
		||||
		   Integral_table Itable,
 | 
			
		||||
		   DP xmin_ref, DP xmax_ref);
 | 
			
		||||
 | 
			
		||||
    void Save (std::ofstream& outfile);
 | 
			
		||||
    void Improve_estimate (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, int Npts_max);
 | 
			
		||||
    void Improve_estimate (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable, int Npts_max);
 | 
			
		||||
    void Improve_estimate (DP (*function) (Vect_DP, Integral_table), Vect_DP& args, int arg_to_integ,
 | 
			
		||||
			   Integral_table Itable, int Npts_max);
 | 
			
		||||
 | 
			
		||||
    void Improve_estimate (DP (*function) (Vect_DP),
 | 
			
		||||
			   Vect_DP& args, int arg_to_integ,
 | 
			
		||||
			   int Npts_max);
 | 
			
		||||
 | 
			
		||||
    void Improve_estimate (DP (*function) (Vect_DP, I_table),
 | 
			
		||||
			   Vect_DP& args, int arg_to_integ,
 | 
			
		||||
			   I_table Itable,
 | 
			
		||||
			   int Npts_max);
 | 
			
		||||
 | 
			
		||||
    void Improve_estimate (DP (*function) (Vect_DP, Integral_table),
 | 
			
		||||
			   Vect_DP& args, int arg_to_integ,
 | 
			
		||||
			   Integral_table Itable,
 | 
			
		||||
			   int Npts_max);
 | 
			
		||||
 | 
			
		||||
    ~Integral_data ();
 | 
			
		||||
 | 
			
		||||
  };
 | 
			
		||||
 | 
			
		||||
  Integral_result Integrate_optimal (DP (*function) (Vect_DP), Vect_DP& args,
 | 
			
		||||
				     int arg_to_integ, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
 | 
			
		||||
  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
 | 
			
		||||
						 I_table Itable, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
 | 
			
		||||
  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, Integral_table Itable), Vect_DP& args, int arg_to_integ,
 | 
			
		||||
						 Integral_table Itable, DP xmin, DP xmax, DP req_rel_prec,
 | 
			
		||||
						 DP req_abs_prec, int max_nr_pts);
 | 
			
		||||
  Integral_result Integrate_optimal (DP (*function) (Vect_DP),
 | 
			
		||||
				     Vect_DP& args, int arg_to_integ,
 | 
			
		||||
				     DP xmin, DP xmax,
 | 
			
		||||
				     DP req_rel_prec, DP req_abs_prec,
 | 
			
		||||
				     int max_nr_pts);
 | 
			
		||||
 | 
			
		||||
  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
 | 
			
		||||
						 I_table Itable, DP xmin, DP xmax, DP req_rel_prec,
 | 
			
		||||
						 DP req_abs_prec, int max_nr_pts, std::ofstream& outfile);
 | 
			
		||||
  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable),
 | 
			
		||||
						 Vect_DP& args, int arg_to_integ,
 | 
			
		||||
						 I_table Itable,
 | 
			
		||||
						 DP xmin, DP xmax,
 | 
			
		||||
						 DP req_rel_prec, DP req_abs_prec,
 | 
			
		||||
						 int max_nr_pts);
 | 
			
		||||
 | 
			
		||||
  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, Integral_table Itable),
 | 
			
		||||
						 Vect_DP& args, int arg_to_integ,
 | 
			
		||||
						 Integral_table Itable,
 | 
			
		||||
						 DP xmin, DP xmax,
 | 
			
		||||
						 DP req_rel_prec, DP req_abs_prec,
 | 
			
		||||
						 int max_nr_pts);
 | 
			
		||||
 | 
			
		||||
  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable),
 | 
			
		||||
						 Vect_DP& args, int arg_to_integ,
 | 
			
		||||
						 I_table Itable,
 | 
			
		||||
						 DP xmin, DP xmax,
 | 
			
		||||
						 DP req_rel_prec, DP req_abs_prec,
 | 
			
		||||
						 int max_nr_pts,
 | 
			
		||||
						 std::ofstream& outfile);
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
  // ******************************** Recursive version:  optimal, complex implementation ******************************
 | 
			
		||||
  // ********************** Recursive version:  optimal, complex implementation ********************
 | 
			
		||||
 | 
			
		||||
  // NB:  function returns complex values but takes real arguments
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			@ -307,23 +366,25 @@ namespace ABACUS {
 | 
			
		|||
    DP xmax;
 | 
			
		||||
 | 
			
		||||
  public:
 | 
			
		||||
    Integral_data_CX (std::complex<DP> (*function_ref) (Vect_DP), Vect_DP& args, int arg_to_integ_ref, DP xmin_ref, DP xmax_ref);
 | 
			
		||||
    Integral_data_CX (std::complex<DP> (*function_ref) (Vect_DP),
 | 
			
		||||
		      Vect_DP& args, int arg_to_integ_ref,
 | 
			
		||||
		      DP xmin_ref, DP xmax_ref);
 | 
			
		||||
 | 
			
		||||
    void Save (std::ofstream& outfile);
 | 
			
		||||
    void Improve_estimate (std::complex<DP> (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, int Npts_max);
 | 
			
		||||
 | 
			
		||||
    void Improve_estimate (std::complex<DP> (*function) (Vect_DP),
 | 
			
		||||
			   Vect_DP& args, int arg_to_integ,
 | 
			
		||||
			   int Npts_max);
 | 
			
		||||
 | 
			
		||||
    ~Integral_data_CX ();
 | 
			
		||||
 | 
			
		||||
  };
 | 
			
		||||
 | 
			
		||||
  Integral_result_CX Integrate_optimal (std::complex<DP> (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax,
 | 
			
		||||
					DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
 | 
			
		||||
  //Integral_result_CX Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
 | 
			
		||||
  //					 I_table Itable, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
 | 
			
		||||
  //Integral_result_CX Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
 | 
			
		||||
  //					 I_table Itable, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts, std::ofstream& outfile);
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
  Integral_result_CX Integrate_optimal (std::complex<DP> (*function) (Vect_DP),
 | 
			
		||||
					Vect_DP& args, int arg_to_integ,
 | 
			
		||||
					DP xmin, DP xmax,
 | 
			
		||||
					DP req_rel_prec, DP req_abs_prec,
 | 
			
		||||
					int max_nr_pts);
 | 
			
		||||
 | 
			
		||||
} // namespace ABACUS
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -8,7 +8,7 @@ Copyright (c) J.-S. Caux.
 | 
			
		|||
 | 
			
		||||
File:  ABACUS_Matrix.h
 | 
			
		||||
 | 
			
		||||
Purpose:  Declares square matrix class.
 | 
			
		||||
Purpose:  Declares square and rectangular matrix classes.
 | 
			
		||||
 | 
			
		||||
***********************************************************/
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -26,7 +26,8 @@ namespace ABACUS {
 | 
			
		|||
    // Refer to GR[6] 8.23
 | 
			
		||||
 | 
			
		||||
    if (x <= 0.0) {
 | 
			
		||||
      std::cout << "Cosine_Integral called with real argument " << x << " <= 0, which is ill-defined because of the branch cut." << std::endl;
 | 
			
		||||
      std::cout << "Cosine_Integral called with real argument "
 | 
			
		||||
		<< x << " <= 0, which is ill-defined because of the branch cut." << std::endl;
 | 
			
		||||
      ABACUSerror("");
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			@ -59,7 +60,8 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
    else { // Use high x power series
 | 
			
		||||
 | 
			
		||||
      // Ci (x) = \frac{\sin x}{x} \sum_{n=0}^\infty (-1)^n (2n)! x^{-2n} - \frac{\cos x}{x} \sum_{n=0}^\infty (-1)^n (2n+1)! x^{-2n-1}
 | 
			
		||||
      // Ci (x) = \frac{\sin x}{x} \sum_{n=0}^\infty (-1)^n (2n)! x^{-2n}
 | 
			
		||||
      //        - \frac{\cos x}{x} \sum_{n=0}^\infty (-1)^n (2n+1)! x^{-2n-1}
 | 
			
		||||
 | 
			
		||||
      int n = 0;
 | 
			
		||||
      DP minonetothen = 1.0;
 | 
			
		||||
| 
						 | 
				
			
			@ -167,7 +169,8 @@ namespace ABACUS {
 | 
			
		|||
    int max_nr_pts = 10000;
 | 
			
		||||
    Integral_result integ_ln_Gamma = Integrate_optimal (ln_Gamma_for_Barnes_G_RE, args, 0, 0.0, z - 1.0, req_rel_prec, req_abs_prec, max_nr_pts);
 | 
			
		||||
 | 
			
		||||
    return(0.5 * (z - 1.0) * (2.0 - z + logtwoPI) + (z - 1.0) * real(ln_Gamma(std::complex<double>(z - 1.0))) - integ_ln_Gamma.integ_est);
 | 
			
		||||
    return(0.5 * (z - 1.0) * (2.0 - z + logtwoPI)
 | 
			
		||||
	   + (z - 1.0) * real(ln_Gamma(std::complex<double>(z - 1.0))) - integ_ln_Gamma.integ_est);
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
} // namespace ABACUS
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -39,7 +39,9 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
  // Inexplicably missing string functions in standard library:
 | 
			
		||||
 | 
			
		||||
  inline std::string replace(const std::string& str, const std::string& from, const std::string& to) {
 | 
			
		||||
  inline std::string replace(const std::string& str,
 | 
			
		||||
			     const std::string& from,
 | 
			
		||||
			     const std::string& to) {
 | 
			
		||||
    std::string repl = str;
 | 
			
		||||
    size_t start_pos = repl.find(from);
 | 
			
		||||
    if(start_pos < std::string::npos)
 | 
			
		||||
| 
						 | 
				
			
			@ -47,7 +49,9 @@ namespace ABACUS {
 | 
			
		|||
    return repl;
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  inline std::string replace_all(const std::string& str, const std::string& from, const std::string& to) {
 | 
			
		||||
  inline std::string replace_all(const std::string& str,
 | 
			
		||||
				 const std::string& from,
 | 
			
		||||
				 const std::string& to) {
 | 
			
		||||
    std::string repl = str;
 | 
			
		||||
    if(from.empty())
 | 
			
		||||
      return repl;
 | 
			
		||||
| 
						 | 
				
			
			@ -125,7 +129,8 @@ namespace ABACUS {
 | 
			
		|||
    int ans = 0;
 | 
			
		||||
 | 
			
		||||
    if (N < 0) {
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if ( N == 1  ||  N == 0) ans = 1;
 | 
			
		||||
| 
						 | 
				
			
			@ -139,7 +144,8 @@ namespace ABACUS {
 | 
			
		|||
    DP ans = 0.0;
 | 
			
		||||
 | 
			
		||||
    if (N < 0) {
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if ( N == 1  ||  N == 0) ans = 0.0;
 | 
			
		||||
| 
						 | 
				
			
			@ -153,7 +159,8 @@ namespace ABACUS {
 | 
			
		|||
    long long int ans = 0;
 | 
			
		||||
 | 
			
		||||
    if (N < 0) {
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if ( N == 1  ||  N == 0) ans = 1;
 | 
			
		||||
| 
						 | 
				
			
			@ -167,7 +174,8 @@ namespace ABACUS {
 | 
			
		|||
    unsigned long long int ans = 0;
 | 
			
		||||
 | 
			
		||||
    if (N < 0) {
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
 | 
			
		||||
      std::cerr << "Error:  factorial of negative number.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if ( N == 1  ||  N == 0) ans = 1;
 | 
			
		||||
| 
						 | 
				
			
			@ -182,7 +190,8 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
    int ans = 0;
 | 
			
		||||
    if (N1 < N2) {
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if (N1 == N2) ans = 1;
 | 
			
		||||
| 
						 | 
				
			
			@ -203,7 +212,8 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
    DP ans = 0.0;
 | 
			
		||||
    if (N1 < N2) {
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if (N1 == N2) ans = 0.0;
 | 
			
		||||
| 
						 | 
				
			
			@ -219,13 +229,14 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
    long long int ans = 0;
 | 
			
		||||
    if (N1 < N2) {
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if (N1 == N2) ans = 1;
 | 
			
		||||
    else if (N1 < 12) ans = fact_lli(N1)/(fact_lli(N2) * fact_lli(N1 - N2));
 | 
			
		||||
    else {
 | 
			
		||||
      // Make sure that N2 is less than or equal to N1/2;  if not, just switch...
 | 
			
		||||
      // Make sure that N2 is less than or equal to N1/2;  if not, just switch
 | 
			
		||||
      int N2_min = min(N2, N1 - N2);
 | 
			
		||||
 | 
			
		||||
      ans = 1;
 | 
			
		||||
| 
						 | 
				
			
			@ -244,13 +255,14 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
    unsigned long long int ans = 0;
 | 
			
		||||
    if (N1 < N2) {
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
 | 
			
		||||
      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
 | 
			
		||||
		<< std::endl;
 | 
			
		||||
      exit(1);
 | 
			
		||||
    }
 | 
			
		||||
    else if (N1 == N2) ans = 1;
 | 
			
		||||
    else if (N1 < 12) ans = fact_ulli(N1)/(fact_ulli(N2) * fact_ulli(N1 - N2));
 | 
			
		||||
    else {
 | 
			
		||||
      // Make sure that N2 is less than or equal to N1/2;  if not, just switch...
 | 
			
		||||
      // Make sure that N2 is less than or equal to N1/2;  if not, just switch
 | 
			
		||||
      int N2_min = min(N2, N1 - N2);
 | 
			
		||||
 | 
			
		||||
      ans = 1;
 | 
			
		||||
| 
						 | 
				
			
			@ -302,6 +314,7 @@ namespace ABACUS {
 | 
			
		|||
    return(-0.5 * II * log((1.0 + II* x)/(1.0 - II* x)));
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
  /****************  Gamma function *******************/
 | 
			
		||||
 | 
			
		||||
  inline std::complex<double> ln_Gamma (std::complex<double> z)
 | 
			
		||||
| 
						 | 
				
			
			@ -325,7 +338,8 @@ namespace ABACUS {
 | 
			
		|||
	+ 0.5384136432509564062961e-7 / (z + 8.0)
 | 
			
		||||
	- 0.4023533141268236372067e-8 / (z + 9.0);
 | 
			
		||||
 | 
			
		||||
      return(0.5 * logtwoPI + (z - 0.5) * log(z + 8.5) - (z + 8.5) + log(series));
 | 
			
		||||
      return(0.5 * logtwoPI + (z - 0.5) * log(z + 8.5)
 | 
			
		||||
	     - (z + 8.5) + log(series));
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    return(log(0.0));  // never called
 | 
			
		||||
| 
						 | 
				
			
			@ -359,7 +373,8 @@ namespace ABACUS {
 | 
			
		|||
      for (int i = 1; i < g+2; ++i)
 | 
			
		||||
	series += p[i]/(z_min_1 + std::complex<double>(i));
 | 
			
		||||
 | 
			
		||||
      return(0.5 * logtwoPI + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
 | 
			
		||||
      return(0.5 * logtwoPI
 | 
			
		||||
	     + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
 | 
			
		||||
	     - (z_min_1 + std::complex<double>(g) + 0.5) + log(series));
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			@ -376,16 +391,25 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
      int g = 7;
 | 
			
		||||
 | 
			
		||||
      double p[9] = { 0.99999999999980993, 676.5203681218851, -1259.1392167224028,
 | 
			
		||||
		      771.32342877765313, -176.61502916214059, 12.507343278686905,
 | 
			
		||||
		      -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
 | 
			
		||||
      double p[9] = {
 | 
			
		||||
	0.99999999999980993,
 | 
			
		||||
	676.5203681218851,
 | 
			
		||||
	-1259.1392167224028,
 | 
			
		||||
	771.32342877765313,
 | 
			
		||||
	-176.61502916214059,
 | 
			
		||||
	12.507343278686905,
 | 
			
		||||
	-0.13857109526572012,
 | 
			
		||||
	9.9843695780195716e-6,
 | 
			
		||||
	1.5056327351493116e-7
 | 
			
		||||
      };
 | 
			
		||||
 | 
			
		||||
      std::complex<double> z_min_1 = z - 1.0;
 | 
			
		||||
      std::complex<double> series = p[0];
 | 
			
		||||
      for (int i = 1; i < g+2; ++i)
 | 
			
		||||
	series += p[i]/(z_min_1 + std::complex<double>(i));
 | 
			
		||||
 | 
			
		||||
      return(0.5 * logtwoPI + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
 | 
			
		||||
      return(0.5 * logtwoPI
 | 
			
		||||
	     + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
 | 
			
		||||
	     - (z_min_1 + std::complex<double>(g) + 0.5) + log(series));
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			@ -396,7 +420,8 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
  inline long long int Partition_Function (int n)
 | 
			
		||||
  {
 | 
			
		||||
    // Returns the value of the partition function p(n), giving the number of partitions of n into integers.
 | 
			
		||||
    // Returns the value of the partition function p(n),
 | 
			
		||||
    // giving the number of partitions of n into integers.
 | 
			
		||||
 | 
			
		||||
    if (n < 0) ABACUSerror("Calling Partition_Function for n < 0.");
 | 
			
		||||
    else if (n == 0 || n == 1) return(1LL);
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -45,7 +45,8 @@ namespace ABACUS {
 | 
			
		|||
  public:
 | 
			
		||||
    Young_Tableau ();    // empty constructor, does nothing
 | 
			
		||||
    Young_Tableau (int Nr, int Nc);  // constructs empty tableau
 | 
			
		||||
    Young_Tableau (int Nr, int Nc, long long int idnr);  // constructs the tableau corresponding to identification number idnr
 | 
			
		||||
    // constructs the tableau with identification number idnr:
 | 
			
		||||
    Young_Tableau (int Nr, int Nc, long long int idnr);
 | 
			
		||||
    Young_Tableau (const Young_Tableau& RefTableau);  // copy constructor
 | 
			
		||||
    Young_Tableau (int Nr, int Nc, const Young_Tableau& RefTableau);
 | 
			
		||||
    Young_Tableau& operator= (const Young_Tableau& RefTableau);  // assignment
 | 
			
		||||
| 
						 | 
				
			
			@ -55,11 +56,13 @@ namespace ABACUS {
 | 
			
		|||
    Young_Tableau& Compute_Map (long long int idnr_to_reach);  // fills the map vector
 | 
			
		||||
    Young_Tableau& Distribute_boxes (int nboxes_to_dist, int level);
 | 
			
		||||
    Young_Tableau& Set_to_id (long long int idnr); // sets the tableau to the one corresponding to idnr
 | 
			
		||||
    Young_Tableau& Set_to_id (long long int idnr, int option); // sets the tableau to the one corresponding to idnr according to rule option
 | 
			
		||||
    // sets the tableau to the one corresponding to idnr according to rule option:
 | 
			
		||||
    Young_Tableau& Set_to_id (long long int idnr, int option);
 | 
			
		||||
    Young_Tableau& Set_Row_L (Vect<int>& Row_Lengths);  // set row lengths
 | 
			
		||||
    Young_Tableau& Set_Col_L_given_Row_L ();   // sets the Col_L array self-consistently
 | 
			
		||||
    Young_Tableau& Set_Row_L_given_Col_L ();   // sets the Col_L array self-consistently
 | 
			
		||||
    long long int Compute_Descendent_id (int option, Vect<int>& Desc_Row_L, int Nrows_Desc, int Ncols_Desc,
 | 
			
		||||
    long long int Compute_Descendent_id (int option, Vect<int>& Desc_Row_L,
 | 
			
		||||
					 int Nrows_Desc, int Ncols_Desc,
 | 
			
		||||
					 const Young_Tableau& RefTableau);
 | 
			
		||||
    Young_Tableau& Compute_id();       // computes the id number of tableau
 | 
			
		||||
    Young_Tableau& Compute_id(int option);       // computes the id number of tableau according to rule option
 | 
			
		||||
| 
						 | 
				
			
			@ -69,7 +72,8 @@ namespace ABACUS {
 | 
			
		|||
    bool Raise_Row (int i);
 | 
			
		||||
    bool Lower_Col (int i);
 | 
			
		||||
    bool Raise_Col (int i);
 | 
			
		||||
    bool Raise_Lowest_Nonzero_Row(); // adds a box to the lowest nonzero length Row, recomputes id, returns true if tableau has changed
 | 
			
		||||
    // adds a box to the lowest nonzero length Row, recomputes id, returns true if tableau has changed:
 | 
			
		||||
    bool Raise_Lowest_Nonzero_Row();
 | 
			
		||||
    bool Raise_Next_to_Lowest_Nonzero_Row(); // same thing, but for Row under lowest nonzero length one.
 | 
			
		||||
    bool Move_Box_from_Col_to_Col (int ifrom, int ito);
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			@ -83,10 +87,12 @@ namespace ABACUS {
 | 
			
		|||
 | 
			
		||||
  inline int Nr_Nonzero_Rows (const Vect<Young_Tableau>& Tableau_ref)
 | 
			
		||||
  {
 | 
			
		||||
    // This function checks the number of rows containing at least one box
 | 
			
		||||
    // in the whole vector of Young tableaux given as argument.
 | 
			
		||||
    // The usefulness is to force descent of states in which only a few
 | 
			
		||||
    // excitations have started dispersing.
 | 
			
		||||
    /*
 | 
			
		||||
    This function checks the number of rows containing at least one box
 | 
			
		||||
    in the whole vector of Young tableaux given as argument.
 | 
			
		||||
    The usefulness is to force descent of states in which only a few
 | 
			
		||||
    excitations have started dispersing.
 | 
			
		||||
    */
 | 
			
		||||
 | 
			
		||||
    int nr_nonzero_rows = 0;
 | 
			
		||||
    for (int i = 0; i < Tableau_ref.size(); ++i)
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -819,7 +819,7 @@ namespace ABACUS {
 | 
			
		|||
		    << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
 | 
			
		||||
	LOG_outfile << "Resulting info: " << scan_info << endl;
 | 
			
		||||
      }
 | 
			
		||||
      LOG_outfile << "Code version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl << endl;
 | 
			
		||||
      LOG_outfile << "ABACUS version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl << endl;
 | 
			
		||||
      LOG_outfile.close();
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
| 
						 | 
				
			
			@ -408,7 +408,7 @@ namespace ABACUS {
 | 
			
		|||
    LOG_outfile << "Refining in parallel mode using " << nr_processors_at_newlevel << " processors."
 | 
			
		||||
		<< endl << "Refining info: " << scan_info_refinement
 | 
			
		||||
		<< endl << "Resulting info:  " << scan_info << endl;
 | 
			
		||||
    LOG_outfile << "Code version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl;
 | 
			
		||||
    LOG_outfile << "ABACUS version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl;
 | 
			
		||||
    LOG_outfile.close();
 | 
			
		||||
 | 
			
		||||
    return;
 | 
			
		||||
| 
						 | 
				
			
			
 | 
			
		|||
		Loading…
	
		Reference in New Issue