qpms/ewaldsf.c

   1 #include "ewald.h"
   2 #include <gsl/gsl_sf_gamma.h>
   3 #include <gsl/gsl_sf_expint.h>
   4 #include <gsl/gsl_sf_exp.h>
   5 #include <gsl/gsl_sf_result.h>
   6 #include <gsl/gsl_errno.h>
   7 #include <gsl/gsl_machine.h> // Maybe I should rather use DBL_EPSILON instead of GSL_DBL_EPSILON.
   8 #include "kahansum.h"
   9 #include <math.h>
  10 #include <complex.h>
  11 //#include <gsl/gsl_integration.h>
  12 #include <gsl/gsl_errno.h>
  13 #include <float.h>
  14 #include <stdbool.h>
  15 #include <Faddeeva.h>
  16 #include "tiny_inlines.h"
  17
  18 #ifndef COMPLEXPART_REL_ZERO_LIMIT
  19 #define COMPLEXPART_REL_ZERO_LIMIT 1e-14
  20 #endif
  21
  22 gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
  23
  24 void IgnoreUnderflowsGSLErrorHandler (const char * reason,
  25                     const char * file,
  26                     const int line,
  27                     const int gsl_errno) {
  28   if (gsl_errno == GSL_EUNDRFLW)
  29     return;
  30
  31   gsl_stream_printf ("ERROR", file, line, reason);
  32
  33   fflush(stdout);
  34   fprintf (stderr, "Underflow-ignoring error handler invoked.\n");
  35   fflush(stderr);
  36
  37   abort();
  38 }
  39
  40 // DLMF 8.7.3 (latter expression) for complex second argument
  41 // BTW if a is large negative, it might take a while to evaluate.
  42 // This can't be used for non-positive integer a due to
  43 // Г(a) in the formula.
  44 int cx_gamma_inc_series_e(const double a, const complex double z, qpms_csf_result * result) {
  45   if (a <= 0 && a == (int) a) {
  46     result->val = NAN + NAN*I;
  47     result->err = NAN;
  48     GSL_ERROR("Undefined for non-positive integer values", GSL_EDOM);
  49   }
  50   gsl_sf_result fullgamma;
  51   int retval = gsl_sf_gamma_e(a, &fullgamma);
  52   if (GSL_EUNDRFLW == retval)
  53     result->err += DBL_MIN;
  54   else if (GSL_SUCCESS != retval){
  55     result->val = NAN + NAN*I; result->err = NAN;
  56     return retval;
  57   }
  58
  59   complex double sumprefac = cpow(z, a) * cexp(-z);
  60   double sumprefac_abs = cabs(sumprefac);
  61   complex double sum, sumc; ckahaninit(&sum, &sumc);
  62   double err, errc; kahaninit(&err, &errc);
  63
  64   bool breakswitch = false;
  65   for (int k = 0; (!breakswitch) && (a + k + 1 <= GSL_SF_GAMMA_XMAX); ++k) {
  66     gsl_sf_result fullgamma_ak;
  67     if (GSL_SUCCESS != (retval = gsl_sf_gamma_e(a+k+1, &fullgamma_ak))) {
  68       result->val = NAN + NAN*I; result->err = NAN;
  69       return retval;
  70     }
  71     complex double summand = - cpow(z, k) / fullgamma_ak.val; // TODO test branch selection here with cimag(z) = -0.0
  72     ckahanadd(&sum, &sumc, summand);
  73     double summanderr = fabs(fullgamma_ak.err * cabs(summand / fullgamma_ak.val));
  74     // TODO add also the rounding error
  75     kahanadd(&err, &errc, summanderr);
  76     // TODO put some smarter cutoff break here?
  77     if (a + k >= 18 && (cabs(summand) < err || cabs(summand) < DBL_EPSILON))
  78       breakswitch = true;
  79   }
  80   sum *= sumprefac; // Not sure if not breaking the Kahan summation here
  81   sumc *= sumprefac;
  82   err *= sumprefac_abs;
  83   errc *= sumprefac_abs;
  84   ckahanadd(&sum, &sumc, 1.);
  85   kahanadd(&err, &errc, DBL_EPSILON);
  86   result->err = cabs(sum) * fullgamma.err + err * fabs(fullgamma.val);
  87   result->val = sum * fullgamma.val; // + sumc*fullgamma.val???
  88   if (breakswitch)
  89     return GSL_SUCCESS;
  90   else GSL_ERROR("Overflow; the absolute value of the z argument is probably too large.", GSL_ETOL); // maybe different error code...
  91 }
  92
  93 /* Continued fraction which occurs in evaluation
  94  * of Q(a,z) or Gamma(a,z).
  95  * Borrowed from GSL and adapted for complex z.
  96  *
  97  *              1   (1-a)/z  1/z  (2-a)/z   2/z  (3-a)/z
  98  *   F(a,z) =  ---- ------- ----- -------- ----- -------- ...
  99  *             1 +   1 +     1 +   1 +      1 +   1 +
 100  *
 101  */
 102 static int
 103 cx_gamma_inc_F_CF(const double a, const complex double z, qpms_csf_result * result)
 104 {
 105   const int    nmax  =  5000;
 106   const double small =  DBL_EPSILON * DBL_EPSILON * DBL_EPSILON;
 107
 108   complex double hn = 1.0;           /* convergent */
 109   complex double Cn = 1.0 / small;
 110   complex double Dn = 1.0;
 111   int n;
 112
 113   /* n == 1 has a_1, b_1, b_0 independent of a,z,
 114      so that has been done by hand                */
 115   for ( n = 2 ; n < nmax ; n++ )
 116   {
 117     complex double an;
 118     complex double delta;
 119
 120     if(n % 2)
 121       an = 0.5*(n-1)/z;
 122     else
 123       an = (0.5*n-a)/z;
 124
 125     Dn = 1.0 + an * Dn;
 126     if ( cabs(Dn) < small )
 127       Dn = small;
 128     Cn = 1.0 + an/Cn;
 129     if ( cabs(Cn) < small )
 130       Cn = small;
 131     Dn = 1.0 / Dn;
 132     delta = Cn * Dn;
 133     hn *= delta;
 134     if(cabs(delta-1.0) < DBL_EPSILON) break;
 135   }
 136
 137   result->val = hn;
 138   result->err = 2.0*GSL_DBL_EPSILON * cabs(hn);
 139   result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * cabs(result->val);
 140
 141   if(n == nmax)
 142     GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER);
 143   else
 144     return GSL_SUCCESS;
 145 }
 146
 147 // Incomplete gamma fuction with complex second argument as continued fraction.
 148 int cx_gamma_inc_CF_e(const double a, const complex double z, qpms_csf_result *result)
 149 {
 150   qpms_csf_result F;
 151   gsl_sf_result pre;
 152   const complex double am1lgz = (a-1.0)*clog(z); // TODO check branches
 153   const int stat_F = cx_gamma_inc_F_CF(a, z, &F);
 154   const int stat_E = gsl_sf_exp_err_e(creal(am1lgz - z), GSL_DBL_EPSILON*cabs(am1lgz), &pre);
 155   complex double cpre = pre.val * cexp(I*cimag(am1lgz - z));// TODO add the error estimate for this
 156   //complex double cpre = cpow(z, a-1) * cexp(-z);
 157
 158
 159   result->val = F.val * cpre;
 160   result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err);
 161   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 162
 163   return GSL_ERROR_SELECT_2(stat_F, stat_E);
 164 }
 165
 166
 167 // Incomplete gamma function for complex second argument.
 168 int complex_gamma_inc_e(double a, complex double x, int m, qpms_csf_result *result) {
 169   int retval;
 170   if (creal(x) >= 0 &&
 171       (0 == fabs(cimag(x)) || // x is real positive; just use the real fun
 172       fabs(cimag(x)) < fabs(creal(x)) * COMPLEXPART_REL_ZERO_LIMIT)) {
 173     gsl_sf_result real_gamma_inc_result;
 174     retval = gsl_sf_gamma_inc_e(a, creal(x), &real_gamma_inc_result);
 175     result->val = real_gamma_inc_result.val;
 176     result->err = real_gamma_inc_result.err;
 177   } else if (creal(x) >= 0 && cabs(x) > 0.5)
 178     retval = cx_gamma_inc_CF_e(a, x, result);
 179   else if (QPMS_LIKELY(a > 0 || fmod(a, 1.0)))
 180     retval = cx_gamma_inc_series_e(a, x, result);
 181   else
 182   /* FIXME cx_gamma_inc_series_e() probably fails for non-positive integer a.
 183    * This does not matter for 2D lattices in 3D space,
 184    * but it might cause problems in the other cases.
 185    */
 186     QPMS_NOT_IMPLEMENTED("Incomplete Gamma function with non-positive integer a.");
 187   if (m) { // Non-principal branch.
 188     /* This might be sub-optimal, as Γ(a) has probably been already evaluated
 189      * somewhere in the functions called above. */
 190     gsl_sf_result fullgamma;
 191     int retval_fg = gsl_sf_gamma_e(a, &fullgamma);
 192     if (GSL_EUNDRFLW == retval_fg)
 193       fullgamma.err += DBL_MIN;
 194     else if (GSL_SUCCESS != retval_fg){
 195       result->val = NAN + NAN*I; result->err = NAN;
 196       return GSL_ERROR_SELECT_2(retval_fg, retval);
 197     }
 198     complex double f = cexp(2*m*M_PI*a*I);
 199     result->val *= f;
 200     f = -f + 1;
 201     result->err += cabs(f) * fullgamma.err;
 202     result->val += f * fullgamma.val;
 203   }
 204   return retval;
 205 }
 206
 207 // Exponential integral for complex argument; !UNTESTED! and probably not needed, as I expressed everything in terms of inc. gammas anyways.
 208 int complex_expint_n_e(int n, complex double x, qpms_csf_result *result) {
 209   if (creal(x) >= 0 &&
 210     (0 == fabs(cimag(x)) || // x is real positive; just use the real fun
 211     fabs(cimag(x)) < fabs(creal(x)) * COMPLEXPART_REL_ZERO_LIMIT)) {
 212     gsl_sf_result real_expint_result;
 213     int retval = gsl_sf_expint_En_e(n, creal(x), &real_expint_result);
 214     result->val = real_expint_result.val;
 215     result->err = real_expint_result.err;
 216     return retval;
 217   } else {
 218     int retval = complex_gamma_inc_e(-n+1, x, 0, result);
 219     complex double f = cpow(x, 2*n-2);
 220     result->val *= f;
 221     result->err *= cabs(f);
 222     return retval;
 223   }
 224 }
 225
 226 // inspired by GSL's hyperg_2F1_series
 227 int hyperg_2F2_series(const double a, const double b, const double c, const double d,
 228     const double x, gsl_sf_result *result
 229     )
 230 {
 231   double sum_pos = 1.0;
 232   double sum_neg = 0.0;
 233   double del_pos = 1.0;
 234   double del_neg = 0.0;
 235   double del = 1.0;
 236   double del_prev;
 237   double k = 0.0;
 238   int i = 0;
 239
 240   if(fabs(c) < GSL_DBL_EPSILON || fabs(d) < GSL_DBL_EPSILON) {
 241     result->val = NAN;
 242     result->err = INFINITY;
 243     GSL_ERROR ("error", GSL_EDOM);
 244   }
 245
 246   do {
 247     if(++i > 30000) {
 248       result->val  = sum_pos - sum_neg;
 249       result->err  = del_pos + del_neg;
 250       result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
 251       result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
 252       GSL_ERROR ("error", GSL_EMAXITER);
 253     }
 254     del_prev = del;
 255     del *= (a+k)*(b+k) * x / ((c+k) * (d+k) * (k+1.0));  /* Gauss series */
 256
 257     if(del > 0.0) {
 258       del_pos  =  del;
 259       sum_pos +=  del;
 260     }
 261     else if(del == 0.0) {
 262       /* Exact termination (a or b was a negative integer).
 263        */
 264       del_pos = 0.0;
 265       del_neg = 0.0;
 266       break;
 267     }
 268     else {
 269       del_neg  = -del;
 270       sum_neg -=  del;
 271     }
 272
 273     /*
 274      * This stopping criteria is taken from the thesis
 275      * "Computation of Hypergeometic Functions" by J. Pearson, pg. 31
 276      * (http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf)
 277      * and fixes bug #45926
 278      */
 279     if (fabs(del_prev / (sum_pos - sum_neg)) < GSL_DBL_EPSILON &&
 280         fabs(del / (sum_pos - sum_neg)) < GSL_DBL_EPSILON)
 281       break;
 282
 283     k += 1.0;
 284   } while(fabs((del_pos + del_neg)/(sum_pos-sum_neg)) > GSL_DBL_EPSILON);
 285
 286   result->val  = sum_pos - sum_neg;
 287   result->err  = del_pos + del_neg;
 288   result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
 289   result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);
 290
 291   return GSL_SUCCESS;
 292 }
 293
 294 // Complex square root with branch selection
 295 static inline complex double csqrt_branch(complex double x, int xbranch) {
 296   return csqrt(x) * min1pow(xbranch);
 297 }
 298
 299 // The Delta_n factor from [Kambe II], Appendix 3
 300 // \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f]
 301 void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *err,
 302     int maxn, complex double x, int xbranch, complex double z) {
 303   complex double expfac = cexp(-x + 0.25 * z*z / x);
 304   complex double sqrtx = csqrt_branch(x, xbranch); // TODO check carefully, which branch is needed
 305   // These are used in the first two recurrences
 306   complex double w_plus  = Faddeeva_w(+z/(2*sqrtx) + I*sqrtx, 0);
 307   complex double w_minus = Faddeeva_w(-z/(2*sqrtx) + I*sqrtx, 0);
 308   QPMS_ASSERT(maxn >= 0);
 309   if (maxn >= 0)
 310     target[0] = 0.5 * M_SQRTPI * expfac * (w_minus + w_plus);
 311   if (maxn >= 1)
 312     target[1] = I / z * M_SQRTPI * expfac * (w_minus - w_plus);
 313   for(int n = 1; n < maxn; ++n) { // The rest via recurrence
 314     // TODO The cpow(x, 0.5 - n) might perhaps better be replaced with a recurrently computed variant
 315     target[n+1] = -(4 / (z*z)) * (-(0.5 - n) * target[n] + target[n-1] - sqrtx * cpow(x, -n) * expfac);
 316   }
 317   if (err) {
 318     // The error estimates for library math functions are based on
 319     // https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html
 320     // and are not guaranteed to be extremely precise.
 321     // The error estimate for Faddeeva's functions is based on the package's authors at
 322     // http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package
 323     // "we find that the accuracy is typically at at least 13 significant digits in both the real and imaginary parts"
 324     // FIXME the error estimate seems might be off by several orders of magnitude (try parameters x = -3, z = 0.5, maxn=20)
 325     double w_plus_abs = cabs(w_plus), w_minus_abs = cabs(w_minus), expfac_abs = cabs(expfac);
 326     double w_plus_err = w_plus_abs * 1e-13, w_minus_err = w_minus_abs * 1e-13; // LPTODO argument error contrib.
 327     double expfac_err = expfac_abs * (4 * DBL_EPSILON); // LPTODO add argument error contrib.
 328     double z_abs = cabs(z);
 329     double z_err = z_abs * DBL_EPSILON;
 330     double x_abs = cabs(x);
 331     if (maxn >= 0)
 332       err[0] = 0.5 * M_SQRTPI * (expfac_abs * (w_minus_err + w_plus_err) + (w_minus_abs + w_plus_abs) * expfac_err);
 333     if (maxn >= 1)
 334       err[1] = 2 * err[0] / z_abs + cabs(target[1]) * z_err / (z_abs*z_abs);
 335     for(int n = 1; n < maxn; ++n) {
 336       err[n+1] = (2 * cabs(target[n+1]) / z_abs + 4 * ((0.5+n) * err[n] + err[n-1] +
 337             pow(x_abs, 0.5 - n) * (2*DBL_EPSILON * expfac_abs + expfac_err))  // LPTODO not ideal, pow() call is an overkill
 338             ) * z_err / (z_abs*z_abs);
 339     }
 340   }
 341 }
 342
 343
 344 void ewald3_2_sigma_long_Delta_series(complex double *target, double *err,
 345     int maxn, complex double x, int xbranch, complex double z) {
 346   complex double w = 0.25*z*z;
 347   double w_abs = cabs(w);
 348   int maxk;
 349   if (w_abs == 0)
 350     maxk = 0; // Delta is equal to the respective incomplete Gamma functions
 351   else {
 352     // Estimate a suitable maximum k, using Stirling's formula, so that w**maxk / maxk! is less than DBL_EPSILON
 353     // This implementation is quite stupid, but it is still cheap compared to the actual computation, so LPTODO better one
 354     maxk = 1;
 355     double log_w_abs = log(w_abs);
 356     while (maxk * (log_w_abs - log(maxk) + 1) >= -DBL_MANT_DIG)
 357       ++maxk;
 358   }
 359   // TODO asserts on maxn, maxk
 360
 361   complex double *Gammas;
 362   double *Gammas_err = NULL, *Gammas_abs = NULL;
 363   QPMS_CRASHING_CALLOC(Gammas, maxk+maxn+1, sizeof(*Gammas));
 364   if(err) {
 365     QPMS_CRASHING_CALLOC(Gammas_err, maxk+maxn+1, sizeof(*Gammas_err));
 366     QPMS_CRASHING_MALLOC(Gammas_abs, (maxk+maxn+1) * sizeof(*Gammas_abs));
 367   }
 368
 369   for(int j = 0; j <= maxn+maxk; ++j) {
 370     qpms_csf_result g;
 371     QPMS_ENSURE_SUCCESS(complex_gamma_inc_e(0.5-j, x, xbranch, &g));
 372     Gammas[j] = g.val;
 373     if(err) {
 374       Gammas_abs[j] = cabs(g.val);
 375       Gammas_err[j] = g.err;
 376     }
 377   }
 378
 379   for(int n = 0; n <= maxn; ++n) target[n] = 0;
 380   if(err) for(int n = 0; n <= maxn; ++n) err[n] = 0;
 381
 382   complex double wpowk_over_fack = 1.;
 383   double wpowk_over_fack_abs = 1.;
 384   for(int k = 0; k <= maxk; ++k, wpowk_over_fack *= w/k) { // TODO? Kahan sum, Horner's method?
 385     // Also TODO? for small n, continue for higher k if possible/needed
 386     for(int n = 0; n <= maxn; ++n) {
 387       target[n] += Gammas[n+k] * wpowk_over_fack;
 388       if(err) {
 389         // DBL_EPSILON might not be the best estimate here, but...
 390         err[n] += wpowk_over_fack_abs * Gammas_err[n+k] + DBL_EPSILON * Gammas_abs[n+k];
 391         wpowk_over_fack_abs *= w_abs / (k+1);
 392       }
 393     }
 394   }
 395
 396   // TODO add an error estimate for the k-cutoff!!!
 397
 398   free(Gammas);
 399   free(Gammas_err);
 400   free(Gammas_abs);
 401 }
 402
 403
 404 void ewald3_2_sigma_long_Delta(complex double *target, double *err,
 405     int maxn, complex double x, int xbranch, complex double z) {
 406   double absz = cabs(z);
 407   if (absz < 2.) // TODO take into account also the other parameters
 408     ewald3_2_sigma_long_Delta_series(target, err, maxn, x, xbranch, z);
 409   else
 410     ewald3_2_sigma_long_Delta_recurrent(target, err, maxn, x, xbranch, z);
 411 }
 412