polylog.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file polylog.h
24  \brief File defining \ref o2scl::polylog
25 */
26 #ifndef O2SCL_POLYLOG_H
27 #define O2SCL_POLYLOG_H
28 
29 #include <iostream>
30 #include <fstream>
31 #include <string>
32 #include <cmath>
33 
34 #include <gsl/gsl_specfunc.h>
35 
36 #include <o2scl/constants.h>
37 #include <o2scl/err_hnd.h>
38 #include <o2scl/lib_settings.h>
39 //#include <gsl/gsl_sf_dilog.h>
40 #include <o2scl/inte_adapt_cern.h>
41 #include <o2scl/inte_double_exp_boost.h>
42 
43 #ifndef DOXYGEN_NO_O2NS
44 namespace o2scl {
45 #endif
46 
47  /** \brief Polylogarithms (approximate) \f$ Li_n(x)\f$
48 
49  This class is experimental.
50 
51  This gives an approximation to the polylogarithm functions.
52 
53  Only works at present for \f$n=0,1,...,6\f$. Uses GSL library
54  for n=2.
55 
56  Uses linear interpolation for \f$-1<x<0\f$
57  and a series expansion for \f$x<-1\f$
58 
59  \future
60  - Give error estimate?
61  - Improve accuracy?
62  - Use more sophisticated interpolation?
63  - Add the series \f$Li(n,x)=x+2^{-n} x^2+3^{-n} x^3+...\f$
64  for \f$ x \rightarrow 0\f$?
65  - Implement for positive arguments < 1.0
66  - Make another polylog class which implements series acceleration?
67 
68  For reference, there are exact relations
69  \f[
70  \mathrm{Li}_2 \left(\frac{1}{2}\right) =
71  \frac{1}{12}\left[\pi^2-6\left(\ln 2\right)^2\right]
72  \f]
73  \f[
74  \mathrm{Li}_3 \left(\frac{1}{2}\right) =
75  \frac{1}{24}\left[ 4\left(\ln 2\right)^3 - 2 \pi^2 \ln 2 +
76  21 \zeta (3) \right]
77  \f]
78  \f[
79  \mathrm{Li}_{-1} (x) = \frac{x}{\left(1-x\right)^2}
80  \f]
81  \f[
82  \mathrm{Li}_{-2} (x) = \frac{x\left(x+1\right)}{\left(1-x\right)^3}
83  \f]
84 
85  */
86  class polylog {
87 
88  public:
89 
90  /// 0-th order polylogarithm = \f$ x/(1-x)\f$
91  double li0(double x);
92 
93  /// 1-st order polylogarithm = \f$ -\ln(1-x) \f$
94  double li1(double x);
95 
96  /// 2-nd order polylogarithm
97  double li2(double x);
98 
99  /// 3-rd order polylogarithm
100  double li3(double x);
101 
102  /// 4-th order polylogarithm
103  double li4(double x);
104 
105  /// 5-th order polylogarithm
106  double li5(double x);
107 
108  /// 6-th order polylogarithm
109  double li6(double x);
110 
111  polylog();
112 
113  ~polylog();
114 
115  protected:
116 
117 #ifndef DOXYGEN_NO_O2NS
118 
119  double *arg;
120 
121  double *two;
122 
123  double *three;
124 
125  double *four;
126 
127  double *five;
128 
129  double *six;
130 
131  double li2neg1;
132 
133  double li4neg1;
134 
135  double li6neg1;
136 
137 #endif
138 
139  };
140 
141  /** \brief Fermi-Dirac integral
142 
143  This class performs direct computation of the
144  Fermi-Dirac integral
145  \f[
146  F_{a}(\mu) = \int_0^{\infty} \frac{x^a}{1+\exp^{x-\mu}}
147  \f]
148  using \ref o2scl::inte_adapt_cern .
149  */
150  template<class inte_t, class fp_t=double> class fermi_dirac_integ_tl {
151 
152  protected:
153 
154  /// Internal function type
155  typedef std::function<fp_t(fp_t)> func_t;
156 
157  /** \brief The Fermi-Dirac function
158  */
159  fp_t obj_func(fp_t x, fp_t a, fp_t mu) {
160  fp_t res;
161  if (x==0.0) res=0;
162  else res=pow(x,a)/(1.0+exp(x-mu));
163  if (!std::isfinite(res)) {
164  std::cout << x << " " << a << " " << mu << " " << x-mu << " "
165  << res << std::endl;
166  }
167  return res;
168  }
169 
170  public:
171 
172  /** \brief The integrator
173  */
174  inte_t iiu;
175 
176  /** \brief Compute the integral, storing the result in
177  \c res and the error in \c err
178  */
179  void calc_err(fp_t a, fp_t mu, fp_t &res, fp_t &err) {
180  func_t f=
181  std::bind(std::mem_fn<fp_t(fp_t,fp_t,fp_t)>
183  this,std::placeholders::_1,a,mu);
184  iiu.integ_err(f,res,err);
185  return;
186  }
187 
188  };
189 
190  /** \brief Double-precision version of \ref o2scl::fermi_dirac_integ_tl
191  */
192  typedef fermi_dirac_integ_tl
193  <inte_iu<funct,inte_adapt_cern
194  <funct,inte_gauss56_cern<funct,double,
195  inte_gauss56_coeffs_double>,
196  100,double>,double>,double> fermi_dirac_integ_double;
197 
198 #if defined(O2SCL_LD_TYPES) || defined(DOXYGEN)
199 
200  /** \brief Long double version of \ref o2scl::fermi_dirac_integ_tl
201  */
202  typedef fermi_dirac_integ_tl
204  <funct_ld,inte_gauss56_cern<funct_ld,long double,
206  1000,long double>,long double>,long double>
208 
209 #endif
210 
211  /** \brief Compute the fermion integrals for a non-relativistic
212  particle using the GSL functions
213  */
215 
216  public:
217 
218  /** \brief Fermi-Dirac integral of order \f$ 1/2 \f$
219  */
220  double calc_1o2(double y) {
221  return gsl_sf_fermi_dirac_half(y)*sqrt(o2scl_const::pi)/2.0;
222  }
223 
224  /** \brief Fermi-Dirac integral of order \f$ -1/2 \f$
225  */
226  double calc_m1o2(double y) {
227  return gsl_sf_fermi_dirac_mhalf(y)*sqrt(o2scl_const::pi);
228  }
229 
230  /** \brief Fermi-Dirac integral of order \f$ 3/2 \f$
231  */
232  double calc_3o2(double y) {
233  return gsl_sf_fermi_dirac_3half(y)*sqrt(o2scl_const::pi)*0.75;
234  }
235 
236  };
237 
238 #if defined(O2SCL_LD_TYPES) || defined(DOXYGEN)
239 
240 #ifdef O2SCL_NEW_BOOST_INTEGRATION
241 
242  /** \brief Compute the fermion integrals for a non-relativistic
243  particle by directly integrating in long double precision
244  */
245  class fermion_nr_integ_direct {
246 
247  protected:
248 
249  /** \brief The integrator
250  */
251  fermi_dirac_integ_tl<o2scl::inte_exp_sinh_boost
252  <funct_ld,15,long double>,long double> it;
253 
254  public:
255 
256  fermion_nr_integ_direct() {
257  it.iiu.tol_rel=1.0e-18;
258  }
259 
260  /** \brief Fermi-Dirac integral of order \f$ 1/2 \f$
261  */
262  double calc_1o2(double y) {
263  long double y2=y, res, err;
264  it.calc_err(0.5L,y2,res,err);
265  return ((double)res);
266  }
267 
268  /** \brief Fermi-Dirac integral of order \f$ -1/2 \f$
269  */
270  double calc_m1o2(double y) {
271  long double y2=y, res, err;
272  it.calc_err(-0.5L,y2,res,err);
273  return ((double)res);
274  }
275 
276  /** \brief Fermi-Dirac integral of order \f$ 3/2 \f$
277  */
278  double calc_3o2(double y) {
279  long double y2=y, res, err;
280  it.calc_err(1.5L,y2,res,err);
281  return ((double)res);
282  }
283 
284  /** \brief Polylogarithm
285 
286  \warning Doesn't seem to work yet.
287  */
288  double calc_polylog(double s, double y) {
289  long double a=s-1, mu=log(-y), res, err;
290  it.calc_err(a,mu,res,err);
291  return -((double)res);
292  }
293 
294  };
295 
296 #endif
297 
298 #endif
299 
300 #ifndef DOXYGEN_NO_O2NS
301 }
302 #endif
303 
304 #endif
o2scl::fermi_dirac_integ_tl::func_t
std::function< fp_t(fp_t)> func_t
Internal function type.
Definition: polylog.h:155
o2scl::fermi_dirac_integ_tl
Fermi-Dirac integral.
Definition: polylog.h:150
o2scl::polylog::li6
double li6(double x)
6-th order polylogarithm
o2scl::inte_iu
Integrate over .
Definition: inte.h:217
o2scl::inte_exp_sinh_boost
Exp-sinh integration class (Boost)
Definition: inte_double_exp_boost.h:103
o2scl_const::pi
const double pi
Definition: constants.h:65
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::fermion_nr_integ_gsl::calc_1o2
double calc_1o2(double y)
Fermi-Dirac integral of order .
Definition: polylog.h:220
o2scl::fermi_dirac_integ_tl::iiu
inte_t iiu
The integrator.
Definition: polylog.h:174
o2scl::inte_gauss56_cern
5,6-point Gaussian quadrature (CERNLIB)
Definition: inte_gauss56_cern.h:263
o2scl::fermi_dirac_integ_tl::obj_func
fp_t obj_func(fp_t x, fp_t a, fp_t mu)
The Fermi-Dirac function.
Definition: polylog.h:159
o2scl::fermi_dirac_integ_long_double
fermi_dirac_integ_tl< inte_iu< funct_ld, inte_adapt_cern< funct_ld, inte_gauss56_cern< funct_ld, long double, inte_gauss56_coeffs_long_double >, 1000, long double >, long double >, long double > fermi_dirac_integ_long_double
Long double version of o2scl::fermi_dirac_integ_tl.
Definition: polylog.h:207
o2scl::polylog::li3
double li3(double x)
3-rd order polylogarithm
o2scl::fermi_dirac_integ_double
fermi_dirac_integ_tl< inte_iu< funct, inte_adapt_cern< funct, inte_gauss56_cern< funct, double, inte_gauss56_coeffs_double >, 100, double >, double >, double > fermi_dirac_integ_double
Double-precision version of o2scl::fermi_dirac_integ_tl.
Definition: polylog.h:196
o2scl::fermi_dirac_integ_tl::calc_err
void calc_err(fp_t a, fp_t mu, fp_t &res, fp_t &err)
Compute the integral, storing the result in res and the error in err.
Definition: polylog.h:179
o2scl::fermion_nr_integ_gsl
Compute the fermion integrals for a non-relativistic particle using the GSL functions.
Definition: polylog.h:214
o2scl::polylog
Polylogarithms (approximate) .
Definition: polylog.h:86
o2scl::polylog::li4
double li4(double x)
4-th order polylogarithm
o2scl::funct_ld
std::function< long double(long double)> funct_ld
One-dimensional function typedef in src/base/funct.h.
Definition: funct.h:51
o2scl::fermion_nr_integ_gsl::calc_m1o2
double calc_m1o2(double y)
Fermi-Dirac integral of order .
Definition: polylog.h:226
o2scl::polylog::li2
double li2(double x)
2-nd order polylogarithm
o2scl::polylog::li5
double li5(double x)
5-th order polylogarithm
o2scl::inte_adapt_cern
Adaptive integration (CERNLIB)
Definition: inte_adapt_cern.h:243
o2scl::fermion_nr_integ_gsl::calc_3o2
double calc_3o2(double y)
Fermi-Dirac integral of order .
Definition: polylog.h:232
o2scl::funct
std::function< double(double)> funct
One-dimensional function typedef in src/base/funct.h.
Definition: funct.h:48
o2scl::polylog::li0
double li0(double x)
0-th order polylogarithm =
o2scl::polylog::li1
double li1(double x)
1-st order polylogarithm =
o2scl::inte_gauss56_coeffs_long_double
Integration weights and abcissas for o2scl::inte_gauss56_cern in long double precision.
Definition: inte_gauss56_cern.h:103

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).