123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391 |
- /**********************************************************
-
- 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
|