ABACUS/include/JSC_util.h

/**********************************************************

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_