123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499 |
- /**********************************************************
-
- This software is part of J.-S. Caux's ABACUS library.
-
- Copyright (c) J.-S. Caux.
-
- -----------------------------------------------------------
-
- File: ABACUS_util.h
-
- Purpose: Defines basic math functions.
-
- ***********************************************************/
-
- #ifndef ABACUS_UTIL_H
- #define ABACUS_UTIL_H
-
- #include "ABACUS.h"
-
-
- typedef double DP;
-
- // Global constants
-
- const double PI = 3.141592653589793238462643;
- const double sqrtPI = sqrt(PI);
- const double twoPI = 2.0*PI;
- const double logtwoPI = log(twoPI);
- const double Euler_Mascheroni = 0.577215664901532860606;
- const double Gamma_min_0p5 = -2.0 * sqrt(PI);
- const std::complex<double> II(0.0,1.0); // Shorthand for i
-
- const DP MACHINE_EPS = std::numeric_limits<DP>::epsilon();
- const DP MACHINE_EPS_SQ = pow(MACHINE_EPS, 2.0);
-
- // Now for some basic math utilities:
-
- namespace ABACUS {
-
- // Inexplicably missing string functions in standard library:
-
- inline std::string DP_to_string (DP value) {
- std::stringstream s;
- s << std::setprecision(16) << value;
- return s.str();
- }
-
- 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)
- repl.replace(start_pos, from.length(), to);
- return repl;
- }
-
- 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;
- size_t start_pos = 0;
- while((start_pos = repl.find(from, start_pos)) != std::string::npos) {
- repl.replace(start_pos, from.length(), to);
- start_pos += to.length();
- }
- return repl;
- }
-
-
- // File checks:
-
- inline unsigned int count_lines(std::string filename)
- {
- std::ifstream infile(filename);
- return(std::count(std::istreambuf_iterator<char>(infile),
- std::istreambuf_iterator<char>(), '\n'));
- }
-
- inline bool file_exists (const char* filename)
- {
- std::fstream file;
- file.open(filename);
- bool exists = !file.fail();
- file.close();
-
- return(exists);
- }
-
- // Error handler:
-
- inline void ABACUSerror (const std::string error_text)
- // my error handler
- {
- std::cerr << "Run-time error... " << std::endl;
- std::cerr << error_text << std::endl;
- std::cerr << "Exiting to system..." << std::endl;
- exit(1);
- }
-
- struct Divide_by_zero {};
-
-
- // Basics: min, max, fabs
-
- template<class T>
- inline const T max (const T& a, const T& b) { return a > b ? (a) : (b); }
-
- template<class T>
- inline const T min (const T& a, const T& b) { return a > b ? (b) : (a); }
-
- template<class T>
- inline const T fabs (const T& a) { return a >= 0 ? (a) : (-a); }
-
- inline long long int pow_lli (const long long int& base, const int& exp)
- {
- long long int answer = base;
- if (exp == 0) answer = 1LL;
- else for (int i = 1; i < exp; ++i) answer *= base;
- return(answer);
- }
-
- inline unsigned long long int pow_ulli (const unsigned long long int& base, const int& exp)
- {
- unsigned long long int answer = base;
- if (exp == 0) answer = 1ULL;
- for (int i = 1; i < exp; ++i) answer *= base;
- return(answer);
- }
-
- inline int fact (const int& N)
- {
- int ans = 0;
-
- if (N < 0) {
- std::cerr << "Error: factorial of negative number. Exited."
- << std::endl;
- exit(1);
- }
- else if ( N == 1 || N == 0) ans = 1;
- else ans = N * fact(N-1);
-
- return(ans);
- }
-
- inline DP ln_fact (const int& N)
- {
- DP ans = 0.0;
-
- if (N < 0) {
- std::cerr << "Error: factorial of negative number. Exited."
- << std::endl;
- exit(1);
- }
- else if ( N == 1 || N == 0) ans = 0.0;
- else ans = log(DP(N)) + ln_fact(N-1);
-
- return(ans);
- }
-
- inline long long int fact_lli (const int& N)
- {
- long long int ans = 0;
-
- if (N < 0) {
- std::cerr << "Error: factorial of negative number. Exited."
- << std::endl;
- exit(1);
- }
- else if ( N == 1 || N == 0) ans = 1;
- else ans = fact_lli(N-1) * N;
-
- return(ans);
- }
-
- inline long long int fact_ulli (const int& N)
- {
- unsigned long long int ans = 0;
-
- if (N < 0) {
- std::cerr << "Error: factorial of negative number. Exited."
- << std::endl;
- exit(1);
- }
- else if ( N == 1 || N == 0) ans = 1;
- else ans = fact_ulli(N-1) * N;
-
- return(ans);
- }
-
- inline int choose (const int& N1, const int& N2)
- {
- // returns N1 choose N2
-
- int ans = 0;
- if (N1 < N2) {
- 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(N1)/(fact(N2) * fact(N1 - N2));
- else {
- ans = 1;
- int mult = N1;
- while (mult > max(N2, N1 - N2)) ans *= mult--;
- ans /= fact(min(N2, N1 - N2));
- }
-
- return(ans);
- }
-
- inline DP ln_choose (const int& N1, const int& N2)
- {
- // returns the log of N1 choose N2
-
- DP ans = 0.0;
- if (N1 < N2) {
- std::cout << "Error: N1 smaller than N2 in choose. Exited."
- << std::endl;
- exit(1);
- }
- else if (N1 == N2) ans = 0.0;
- else ans = ln_fact(N1) - ln_fact(N2) - ln_fact(N1 - N2);
-
- return(ans);
- }
-
-
- inline long long int choose_lli (const int& N1, const int& N2)
- {
- // returns N1 choose N2
-
- long long int ans = 0;
- if (N1 < N2) {
- 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
- int N2_min = min(N2, N1 - N2);
-
- ans = 1;
- for (int i = 0; i < N2_min; ++i) {
- ans *= (N1 - i);
- ans /= i + 1;
- }
- }
-
- return(ans);
- }
-
- inline unsigned long long int choose_ulli (const int& N1, const int& N2)
- {
- // returns N1 choose N2
-
- unsigned long long int ans = 0;
- if (N1 < N2) {
- 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
- int N2_min = min(N2, N1 - N2);
-
- ans = 1;
- for (int i = 0; i < N2_min; ++i) {
- ans *= (N1 - i);
- ans /= i + 1;
- }
- }
-
- return(ans);
- }
-
- inline DP SIGN (const DP &a, const DP &b)
- {
- return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
- }
-
- inline DP sign_of (const DP& a)
- {
- return (a >= 0.0 ? 1.0 : -1.0);
- }
-
- inline int sgn_int (const int& a)
- {
- return (a >= 0) ? 1 : -1;
- }
-
- inline int sgn_DP (const DP& a)
- {
- return (a >= 0) ? 1 : -1;
- }
-
- template<class T>
- inline void SWAP (T& a, T& b) {T dum = a; a = b; b = dum;}
-
- inline int kronecker (int a, int b)
- {
- return a == b ? 1 : 0;
- }
-
- template<class T>
- inline bool is_nan (const T& a)
- {
- return(!((a < T(0.0)) || (a >= T(0.0))));
- }
-
- inline std::complex<DP> atan_cx(const std::complex<DP>& x)
- {
- 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)
- {
- // Implementation of Lanczos method with g = 9.
- // Coefficients from Godfrey 2001.
-
- if (real(z) < 0.5) return(log(PI/(sin(PI*z))) - ln_Gamma(1.0 - z));
-
- else {
-
- std::complex<double> series = 1.000000000000000174663
- + 5716.400188274341379136/z
- - 14815.30426768413909044/(z + 1.0)
- + 14291.49277657478554025/(z + 2.0)
- - 6348.160217641458813289/(z + 3.0)
- + 1301.608286058321874105/(z + 4.0)
- - 108.1767053514369634679/(z + 5.0)
- + 2.605696505611755827729/(z + 6.0)
- - 0.7423452510201416151527e-2 / (z + 7.0)
- + 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(log(0.0)); // never called
- }
-
- inline std::complex<double> ln_Gamma_old (std::complex<double> z)
- {
- // Implementation of Lanczos method with g = 9.
- // Coefficients from Godfrey 2001.
-
- if (real(z) < 0.5) return(log(PI/(sin(PI*z))) - ln_Gamma(1.0 - z));
-
- else {
-
- int g = 9;
-
- double p[11] = { 1.000000000000000174663,
- 5716.400188274341379136,
- -14815.30426768413909044,
- 14291.49277657478554025,
- -6348.160217641458813289,
- 1301.608286058321874105,
- -108.1767053514369634679,
- 2.605696505611755827729,
- -0.7423452510201416151527e-2,
- 0.5384136432509564062961e-7,
- -0.4023533141268236372067e-8 };
-
- 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)
- - (z_min_1 + std::complex<double>(g) + 0.5) + log(series));
- }
-
- return(log(0.0)); // never called
- }
-
- inline std::complex<double> ln_Gamma_2 (std::complex<double> z)
- {
- // Implementation of Lanczos method with g = 7.
-
- if (real(z) < 0.5) return(log(PI/(sin(PI*z)) - ln_Gamma(1.0 - z)));
-
- else {
-
- 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
- };
-
- 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)
- - (z_min_1 + std::complex<double>(g) + 0.5) + log(series));
- }
-
- return(log(0.0)); // never called
- }
-
- /********** Partition numbers **********/
-
- 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.
-
- if (n < 0) ABACUSerror("Calling Partition_Function for n < 0.");
- else if (n == 0 || n == 1) return(1LL);
- else if (n == 2) return(2LL);
- else if (n == 3) return(3LL);
-
- else { // do recursion using pentagonal numbers
- long long int pn = 0LL;
- int pentnrplus, pentnrmin; // pentagonal numbers
- for (int i = 1; true; ++i) {
- pentnrplus = (i * (3*i - 1))/2;
- pentnrmin = (i * (3*i + 1))/2;
- if (n - pentnrplus >= 0) pn += (i % 2 ? 1LL : -1LL) * Partition_Function (n - pentnrplus);
- if (n - pentnrmin >= 0) pn += (i % 2 ? 1LL : -1LL) * Partition_Function (n - pentnrmin);
- else break;
- }
- return(pn);
- }
- return(-1LL); // never called
- }
-
-
- /********** Sorting **********/
-
- template <class T>
- void QuickSort (T* V, int l, int r)
- {
- int i = l, j = r;
- T pivot = V[l + (r-l)/2];
-
- while (i <= j) {
- while (V[i] < pivot) i++;
- while (V[j] > pivot) j--;
- if (i <= j) {
- std::swap(V[i],V[j]);
- i++;
- j--;
- }
- };
-
- if (l < j) QuickSort(V, l, j);
- if (i < r) QuickSort(V, i, r);
- }
-
- template <class T>
- void QuickSort (T* V, int* index, int l, int r)
- {
- int i = l, j = r;
- T pivot = V[l + (r-l)/2];
-
- while (i <= j) {
- while (V[i] < pivot) i++;
- while (V[j] > pivot) j--;
- if (i <= j) {
- std::swap(V[i],V[j]);
- std::swap(index[i],index[j]);
- i++;
- j--;
- }
- };
-
- if (l < j) QuickSort(V, index, l, j);
- if (i < r) QuickSort(V, index, i, r);
- }
-
-
- } // namespace ABACUS
-
- #endif
|