392 Zeilen
9.3 KiB
C++
392 Zeilen
9.3 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: ABACUS_Integ.h
|
|
|
|
Purpose: Declares combinatorics-related classes and functions.
|
|
|
|
***********************************************************/
|
|
|
|
#ifndef ABACUS_INTEG_H
|
|
#define ABACUS_INTEG_H
|
|
|
|
#include "ABACUS.h"
|
|
|
|
namespace ABACUS {
|
|
|
|
|
|
//********************** Class Domain ************************
|
|
|
|
template<class T>
|
|
class Domain {
|
|
|
|
private:
|
|
Vect<T> bdry;
|
|
|
|
public:
|
|
Domain () : bdry(Vect<T>(2))
|
|
{
|
|
bdry[0] = T(0);
|
|
bdry[1] = T(0);
|
|
}
|
|
|
|
public:
|
|
Domain (T xmin_ref, T xmax_ref) : bdry(Vect<T>(2))
|
|
{
|
|
if (xmax_ref < xmin_ref) ABACUSerror("Use xmax > xmin in Domain.");
|
|
|
|
bdry[0] = xmin_ref;
|
|
bdry[1] = xmax_ref;
|
|
}
|
|
|
|
public:
|
|
inline T xmin (int i)
|
|
{
|
|
if (i > bdry.size()/2) ABACUSerror("i index too high in Domain::xmin.");
|
|
|
|
return(bdry[2*i]);
|
|
}
|
|
|
|
public:
|
|
inline T xmax (int i)
|
|
{
|
|
if (i > bdry.size()/2) ABACUSerror("i index too high in Domain::xmax.");
|
|
|
|
return(bdry[2*i + 1]);
|
|
}
|
|
|
|
public:
|
|
inline int Ndomains ()
|
|
{
|
|
return(bdry.size()/2);
|
|
}
|
|
|
|
public:
|
|
void Include (T xmin_ref, T xmax_ref) {
|
|
|
|
// Determine the indices of xmin_ref & xmax_ref
|
|
int xmin_reg = -1;
|
|
int xmax_reg = -1;
|
|
for (int i = 0; i < bdry.size(); ++i) {
|
|
if ((i+1) % 2 && bdry[i] <= xmin_ref) xmin_reg++;
|
|
if (i % 2 && bdry[i] < xmin_ref) xmin_reg++;
|
|
if ((i+1) % 2 && bdry[i] <= xmax_ref) xmax_reg++;
|
|
if (i % 2 && bdry[i] < xmax_ref) xmax_reg++;
|
|
}
|
|
|
|
Vect<T> new_bdry(bdry.size() + 2 * ((xmin_reg % 2 && xmax_reg % 2) - (xmax_reg - xmin_reg)/2));
|
|
|
|
int ishift = 0;
|
|
for (int i = 0; i <= xmin_reg; ++i) new_bdry[i] = bdry[i];
|
|
if (xmin_reg % 2) {
|
|
new_bdry[xmin_reg + 1] = xmin_ref;
|
|
ishift++;
|
|
if (xmax_reg % 2) {
|
|
new_bdry[xmin_reg + 2] = xmax_ref;
|
|
ishift++;
|
|
}
|
|
}
|
|
else if ((xmin_reg + 1) % 2 && xmax_reg % 2) {
|
|
new_bdry[xmin_reg + 1] = xmax_ref;
|
|
ishift++;
|
|
}
|
|
for (int i = xmin_reg + ishift + 1; i < new_bdry.size(); ++i)
|
|
new_bdry[i] = bdry[xmax_reg - xmin_reg - ishift + i];
|
|
|
|
bdry = new_bdry;
|
|
|
|
return;
|
|
}
|
|
|
|
public:
|
|
void Exclude (T xmin_ref, T xmax_ref) {
|
|
|
|
// Determine the indices of xmin_ref & xmax_ref
|
|
int xmin_reg = -1;
|
|
int xmax_reg = -1;
|
|
for (int i = 0; i < bdry.size(); ++i) {
|
|
if ((i+1) % 2 && bdry[i] <= xmin_ref) xmin_reg++;
|
|
if (i % 2 && bdry[i] < xmin_ref) xmin_reg++;
|
|
if ((i+1) % 2 && bdry[i] <= xmax_ref) xmax_reg++;
|
|
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));
|
|
|
|
int ishift = 0;
|
|
for (int i = 0; i <= xmin_reg; ++i) new_bdry[i] = bdry[i];
|
|
if ((xmin_reg + 1) % 2) {
|
|
new_bdry[xmin_reg + 1] = xmin_ref;
|
|
ishift++;
|
|
if ((xmax_reg + 1) % 2) {
|
|
new_bdry[xmin_reg + 2] = xmax_ref;
|
|
ishift++;
|
|
}
|
|
}
|
|
else if (xmin_reg % 2 && (xmax_reg + 1) % 2) {
|
|
new_bdry[xmin_reg + 1] = xmax_ref;
|
|
ishift++;
|
|
}
|
|
for (int i = xmin_reg + ishift + 1; i < new_bdry.size(); ++i)
|
|
new_bdry[i] = bdry[xmax_reg - xmin_reg - ishift + i];
|
|
|
|
bdry = new_bdry;
|
|
|
|
return;
|
|
}
|
|
|
|
};
|
|
|
|
template<class T>
|
|
std::ostream& operator<< (std::ostream& s, Domain<T> dom)
|
|
{
|
|
for (int i = 0; i < dom.Ndomains(); ++i) {
|
|
if (i > 0) s << std::endl;
|
|
s << dom.xmin(i) << "\t" << dom.xmax(i);
|
|
}
|
|
return(s);
|
|
}
|
|
|
|
|
|
// ********************************* struct I_table ************************
|
|
|
|
struct I_table {
|
|
|
|
DP (*function) (DP, DP);
|
|
int Nvals;
|
|
DP rhomin;
|
|
DP rhomax;
|
|
DP alpha;
|
|
DP logalpha;
|
|
DP prec;
|
|
DP* rho_tbl;
|
|
DP* I_tbl;
|
|
|
|
I_table (DP (*function) (DP, DP), DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec);
|
|
DP Return_val (DP req_rho);
|
|
void Save ();
|
|
bool Load (DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec);
|
|
|
|
};
|
|
|
|
|
|
// ********************************* struct Integral_table ************************
|
|
|
|
struct Integral_table {
|
|
|
|
DP (*function) (DP, DP, int);
|
|
int Nvals;
|
|
DP rhomin;
|
|
DP rhomax;
|
|
DP alpha;
|
|
DP logalpha;
|
|
DP prec;
|
|
int maxnrpts;
|
|
DP* rho_tbl;
|
|
DP* I_tbl;
|
|
|
|
Integral_table (DP (*function) (DP, DP, int), const char* filenameprefix_ref, DP rhomin_ref,
|
|
DP rhomax_ref, int Nvals_ref, DP req_prec, int max_nr_pts);
|
|
DP Return_val (DP req_rho);
|
|
void Save (const char* filenameprefix);
|
|
bool Load (const char* filenameprefix, DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec, int max_nr_pts);
|
|
|
|
};
|
|
|
|
|
|
// ******************************** 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 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);
|
|
|
|
|
|
|
|
// ******************************** Recursive version: optimal ******************************
|
|
|
|
struct data_pt {
|
|
DP x;
|
|
DP f;
|
|
DP dx;
|
|
};
|
|
|
|
struct Integral_result {
|
|
DP integ_est;
|
|
DP abs_prec;
|
|
DP rel_prec;
|
|
int n_vals;
|
|
};
|
|
|
|
std::ostream& operator<< (std::ostream& s, const Integral_result& res);
|
|
|
|
class Integral_data {
|
|
|
|
private:
|
|
data_pt* data;
|
|
DP* abs_d2f_dx; // second derivative * dx
|
|
DP max_abs_d2f_dx; //
|
|
|
|
public:
|
|
Integral_result integ_res;
|
|
|
|
public:
|
|
DP xmin;
|
|
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);
|
|
|
|
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);
|
|
|
|
~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_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 ********************
|
|
|
|
// NB: function returns complex values but takes real arguments
|
|
|
|
struct data_pt_CX {
|
|
DP x;
|
|
std::complex<DP> f;
|
|
DP dx;
|
|
};
|
|
|
|
struct Integral_result_CX {
|
|
std::complex<DP> integ_est;
|
|
DP abs_prec;
|
|
DP rel_prec;
|
|
int n_vals;
|
|
};
|
|
|
|
class Integral_data_CX {
|
|
|
|
private:
|
|
data_pt_CX* data;
|
|
DP* abs_d2f_dx; // second derivative * dx
|
|
DP max_abs_d2f_dx; //
|
|
|
|
public:
|
|
Integral_result_CX integ_res;
|
|
|
|
public:
|
|
DP xmin;
|
|
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);
|
|
|
|
void Save (std::ofstream& outfile);
|
|
|
|
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);
|
|
|
|
} // namespace ABACUS
|
|
|
|
#endif
|