486 строки
11 KiB
C++
486 строки
11 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c).
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: JSC_util.h
|
|
|
|
Purpose: Defines basic math functions.
|
|
|
|
***********************************************************/
|
|
|
|
#ifndef _JSC_UTIL_H_
|
|
#define _JSC_UTIL_H_
|
|
|
|
#include "JSC.h"
|
|
|
|
using namespace std;
|
|
|
|
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 complex<double> II(0.0,1.0); // Shorthand for i
|
|
|
|
const DP MACHINE_EPS = numeric_limits<DP>::epsilon();
|
|
const DP MACHINE_EPS_SQ = pow(MACHINE_EPS, 2.0);
|
|
|
|
// Now for some basic math utilities:
|
|
|
|
namespace JSC {
|
|
|
|
// File checks:
|
|
|
|
inline bool file_exists (const char* filename)
|
|
{
|
|
fstream file;
|
|
file.open(filename);
|
|
bool exists = !file.fail();
|
|
file.close();
|
|
|
|
return(exists);
|
|
}
|
|
|
|
// Error handler:
|
|
|
|
inline void JSCerror (const string error_text)
|
|
// my error handler
|
|
{
|
|
cerr << "Run-time error... " << endl;
|
|
cerr << error_text << endl;
|
|
cerr << "Exiting to system..." << 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) {
|
|
cerr << "Error: factorial of negative number. Exited." << 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) {
|
|
cerr << "Error: factorial of negative number. Exited." << 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) {
|
|
cerr << "Error: factorial of negative number. Exited." << 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) {
|
|
cerr << "Error: factorial of negative number. Exited." << 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) {
|
|
cout << "Error: N1 smaller than N2 in choose. Exited." << 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) {
|
|
cout << "Error: N1 smaller than N2 in choose. Exited." << 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) {
|
|
cout << "Error: N1 smaller than N2 in choose. Exited." << 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) {
|
|
cout << "Error: N1 smaller than N2 in choose. Exited." << 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 complex<DP> atan_cx(const complex<DP>& x)
|
|
{
|
|
return(-0.5 * II * log((1.0 + II* x)/(1.0 - II* x)));
|
|
}
|
|
|
|
/**************** Gamma function *******************/
|
|
|
|
inline complex<double> ln_Gamma (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 {
|
|
|
|
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 complex<double> ln_Gamma_old (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 };
|
|
|
|
complex<double> z_min_1 = z - 1.0;
|
|
complex<double> series = p[0];
|
|
for (int i = 1; i < g+2; ++i)
|
|
series += p[i]/(z_min_1 + complex<double>(i));
|
|
|
|
return(0.5 * logtwoPI + (z_min_1 + 0.5) * log(z_min_1 + complex<double>(g) + 0.5) - (z_min_1 + complex<double>(g) + 0.5) + log(series));
|
|
}
|
|
|
|
return(log(0.0)); // never called
|
|
}
|
|
|
|
inline complex<double> ln_Gamma_2 (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};
|
|
|
|
complex<double> z_min_1 = z - 1.0;
|
|
complex<double> series = p[0];
|
|
for (int i = 1; i < g+2; ++i)
|
|
series += p[i]/(z_min_1 + complex<double>(i));
|
|
|
|
return(0.5 * logtwoPI + (z_min_1 + 0.5) * log(z_min_1 + complex<double>(g) + 0.5) - (z_min_1 + 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) JSCerror("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;
|
|
//cout << "\ti = " << i << "pnrp = " << pentnrplus << "\tpnrm = " << pentnrmin << endl;
|
|
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 item_type>
|
|
void QuickSort (item_type* a, int l, int r)
|
|
{
|
|
static item_type m;
|
|
static int j;
|
|
int i;
|
|
|
|
if (r > l) {
|
|
m = a[r]; i = l-1; j = r;
|
|
for (;;) {
|
|
while (a[++i] < m);
|
|
while (a[--j] > m);
|
|
if (i >= j) break;
|
|
std::swap(a[i], a[j]);
|
|
}
|
|
std::swap(a[i], a[r]);
|
|
QuickSort(a, l, i-1);
|
|
QuickSort(a, i+1, r);
|
|
}
|
|
}
|
|
*/
|
|
|
|
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 item_type>
|
|
void QuickSort (item_type* a, int* idx, int l, int r)
|
|
{
|
|
static item_type m;
|
|
static int j;
|
|
int i;
|
|
|
|
if (r > l) {
|
|
m = a[r]; i = l-1; j = r;
|
|
for (;;) {
|
|
while (a[++i] < m);
|
|
while (a[--j] > m);
|
|
if (i >= j) break;
|
|
std::swap(a[i], a[j]);
|
|
std::swap(idx[i], idx[j]);
|
|
}
|
|
std::swap(a[i], a[r]);
|
|
std::swap(idx[i], idx[r]);
|
|
QuickSort(a, idx, l, i-1);
|
|
QuickSort(a, idx, i+1, 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 JSC
|
|
|
|
#endif // _JS_UTIL_H_
|