123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178 |
- /**********************************************************
-
- 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
|