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

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  ABACUS_Spec_Fns.h

Purpose:  Defines special math functions.

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

#ifndef ABACUS_SPEC_FNS_H
#define ABACUS_SPEC_FNS_H

#include "ABACUS.h"

namespace ABACUS {

  inline DP Cosine_Integral (DP x)
  {

    // Returns the Cosine integral -\int_x^\infty dt \frac{\cos t}{t}
    // Refer to GR[6] 8.23

    if (x <= 0.0) {
      std::cout << "Cosine_Integral called with real argument "
		<< x << " <= 0, which is ill-defined because of the branch cut." << std::endl;
      ABACUSerror("");
    }

    else if (x < 15.0) { // Use power series expansion

      // Ci (x) = gamma + \ln x + \sum_{n=1}^\infty (-1)^n x^{2n}/(2n (2n)!).

      int n = 1;
      DP minonetothen = -1.0;
      DP logxtothetwon = 2.0 * log(x);
      DP twologx = 2.0 * log(x);
      DP logtwonfact = log(2.0);

      DP series = minonetothen * exp(logxtothetwon - log(2.0 * n) - logtwonfact);
      DP term_n;

      do {
	n += 1;
	minonetothen *= -1.0;
	logxtothetwon += twologx;
	logtwonfact += log((2.0 * n - 1.0) * 2.0 * n);
	term_n = minonetothen * exp(logxtothetwon - log(2.0 * n) - logtwonfact);
	series += term_n;

      } while (fabs(term_n) > 1.0e-16);

      return(Euler_Mascheroni + log(x) + series);
    }


    else { // Use high x power series

      // Ci (x) = \frac{\sin x}{x} \sum_{n=0}^\infty (-1)^n (2n)! x^{-2n}
      //        - \frac{\cos x}{x} \sum_{n=0}^\infty (-1)^n (2n+1)! x^{-2n-1}

      int n = 0;
      DP minonetothen = 1.0;
      DP logxtothetwon = 0.0;
      DP logxtothetwonplus1 = log(x);
      DP twologx = 2.0 * log(x);
      DP logtwonfact = 0.0;
      DP logtwonplus1fact = 0.0;

      DP series1 = minonetothen * exp(logtwonfact - logxtothetwon);
      DP series2 = minonetothen * exp(logtwonplus1fact - logxtothetwonplus1);

      do {
	n += 1;
	minonetothen *= -1.0;
	logxtothetwon += twologx;
	logxtothetwonplus1 += twologx;
	logtwonfact += log((2.0 * n - 1.0) * 2.0 * n);
	logtwonplus1fact += log(2.0 * n  * (2.0 * n + 1));

	series1 += minonetothen * exp(logtwonfact - logxtothetwon);
	series2 += minonetothen * exp(logtwonplus1fact - logxtothetwonplus1);

      } while (n < 12);

      return((sin(x)/x) * series1 - (cos(x)/x) * series2);

    }


    return(log(-1.0));
  }


  /*********** Jacobi Theta functions *********/

  inline DP Jacobi_Theta_1_q (DP u, DP q) {

    // Uses the summation formula.
    // theta_1 (x) = 2 \sum_{n=1}^\infty (-1)^{n+1} q^{(n-1/2)^2} \sin (2n-1)u
    // in which q is the nome. (GR 8.180.1)
    // We always evaluate to numerical accuracy.

    if (q >= 1.0) ABACUSerror("Jacobi_Theta_1_q function called with q > 1.");


    DP answer = 0.0;
    DP contrib = 0.0;
    DP qtonminhalfsq = pow(q, 0.25);  // this will be q^{(n-1/2)^2}
    DP qtotwon = pow(q, 2.0); // this will be q^{2n}
    DP qsq = q*q;
    int n = 1;

    do {
      contrib = (n % 2 ? 2.0 : -2.0) * qtonminhalfsq * sin((2.0*n - 1.0)*u);
      answer += contrib;
      qtonminhalfsq *= qtotwon;
      qtotwon *= qsq;
      n++;
    } while (fabs(contrib/answer) > MACHINE_EPS);

    return(answer);
  }

  inline std::complex<DP> ln_Jacobi_Theta_1_q (std::complex<DP> u, std::complex<DP> q) {

    // This uses the product representation
    // \theta_1 (x) = 2 q^{1/4} \sin{u} \prod_{n=1}^\infty (1 - 2 q^{2n} \cos 2u + q^{4n}) (1 - q^{2n})
    // (GR 8.181.2)

    std::complex<DP> contrib = 0.0;
    std::complex<DP> qtotwon = q*q; // this will be q^{2n}
    std::complex<DP> qsq = q*q;
    std::complex<DP> twocos2u = 2.0 * cos(2.0*u);
    int n = 1;
    std::complex<DP> answer = log(2.0 * sin(u)) + 0.25 * log(q);

    do {
      contrib = log((1.0 - twocos2u * qtotwon + qtotwon * qtotwon) * (1.0 - qtotwon));
      answer += contrib;
      qtotwon *= qsq;
      n++;
    } while (abs(contrib) > 1.0e-12);

    return(answer);
  }


  /************  Barnes function ************/

  inline DP ln_Gamma_for_Barnes_G_RE (Vect_DP args)
  {
    return(real(ln_Gamma(std::complex<double>(args[0]))));
  }

  inline DP ln_Barnes_G_RE (DP z)
  {
    //  Implementation according to equation (28) of 2004_Adamchik_CPC_157
    //  Restricted to real arguments.

    Vect_DP args (0.0, 2);

    DP req_rel_prec = 1.0e-6;
    DP req_abs_prec = 1.0e-6;
    int max_nr_pts = 10000;
    Integral_result integ_ln_Gamma = Integrate_optimal (ln_Gamma_for_Barnes_G_RE, args, 0, 0.0, z - 1.0, req_rel_prec, req_abs_prec, max_nr_pts);

    return(0.5 * (z - 1.0) * (2.0 - z + logtwoPI)
	   + (z - 1.0) * real(ln_Gamma(std::complex<double>(z - 1.0))) - integ_ln_Gamma.integ_est);
  }

} // namespace ABACUS

#endif