jscaux
/
ABACUS-v1
3
0
Fork 0
Code Issues 1 Pull Requests Releases Wiki Activity
ABACUS-v1/include/ABACUS_Utils.h

469 lines
11 KiB
C++
Raw Blame History

/**********************************************************
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 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
Reference in New Issue View Git Blame Copy Permalink