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
  17 #ifndef COMPLEXPART_REL_ZERO_LIMIT
  18 #define COMPLEXPART_REL_ZERO_LIMIT 1e-14
  19 #endif
  20
  21 gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
  22
  23 void IgnoreUnderflowsGSLErrorHandler (const char * reason,
  24                     const char * file,
  25                     const int line,
  26                     const int gsl_errno) {
  27   if (gsl_errno == GSL_EUNDRFLW)
  28     return;
  29
  30   gsl_stream_printf ("ERROR", file, line, reason);
  31
  32   fflush(stdout);
  33   fprintf (stderr, "Underflow-ignoring error handler invoked.\n");
  34   fflush(stderr);
  35
  36   abort();
  37 }
  38
  39 // DLMF 8.7.3 (latter expression) for complex second argument
  40 // BTW if a is large negative, it might take a while to evaluate.
  41 // This can't be used for non-positive integer a due to
  42 // Г(a) in the formula.
  43 int cx_gamma_inc_series_e(const double a, const complex double z, qpms_csf_result * result) {
  44   if (a <= 0 && a == (int) a) {
  45     result->val = NAN + NAN*I;
  46     result->err = NAN;
  47     GSL_ERROR("Undefined for non-positive integer values", GSL_EDOM);
  48   }
  49   gsl_sf_result fullgamma;
  50   int retval = gsl_sf_gamma_e(a, &fullgamma);
  51   if (GSL_EUNDRFLW == retval)
  52     result->err += DBL_MIN;
  53   else if (GSL_SUCCESS != retval){
  54     result->val = NAN + NAN*I; result->err = NAN;
  55     return retval;
  56   }
  57
  58   complex double sumprefac = cpow(z, a) * cexp(-z);
  59   double sumprefac_abs = cabs(sumprefac);
  60   complex double sum, sumc; ckahaninit(&sum, &sumc);
  61   double err, errc; kahaninit(&err, &errc);
  62
  63   bool breakswitch = false;
  64   for (int k = 0; (!breakswitch) && (a + k + 1 <= GSL_SF_GAMMA_XMAX); ++k) {
  65     gsl_sf_result fullgamma_ak;
  66     if (GSL_SUCCESS != (retval = gsl_sf_gamma_e(a+k+1, &fullgamma_ak))) {
  67       result->val = NAN + NAN*I; result->err = NAN;
  68       return retval;
  69     }
  70     complex double summand = - cpow(z, k) / fullgamma_ak.val; // TODO test branch selection here with cimag(z) = -0.0
  71     ckahanadd(&sum, &sumc, summand);
  72     double summanderr = fabs(fullgamma_ak.err * cabs(summand / fullgamma_ak.val));
  73     // TODO add also the rounding error
  74     kahanadd(&err, &errc, summanderr);
  75     // TODO put some smarter cutoff break here?
  76     if (a + k >= 18 && (cabs(summand) < err || cabs(summand) < DBL_EPSILON))
  77       breakswitch = true;
  78   }
  79   sum *= sumprefac; // Not sure if not breaking the Kahan summation here
  80   sumc *= sumprefac;
  81   err *= sumprefac_abs;
  82   errc *= sumprefac_abs;
  83   ckahanadd(&sum, &sumc, 1.);
  84   kahanadd(&err, &errc, DBL_EPSILON);
  85   result->err = cabs(sum) * fullgamma.err + err * fabs(fullgamma.val);
  86   result->val = sum * fullgamma.val; // + sumc*fullgamma.val???
  87   if (breakswitch)
  88     return GSL_SUCCESS;
  89   else GSL_ERROR("Overflow; the absolute value of the z argument is probably too large.", GSL_ETOL); // maybe different error code...
  90 }
  91
  92 /* Continued fraction which occurs in evaluation
  93  * of Q(a,z) or Gamma(a,z).
  94  * Borrowed from GSL and adapted for complex z.
  95  *
  96  *              1   (1-a)/z  1/z  (2-a)/z   2/z  (3-a)/z
  97  *   F(a,z) =  ---- ------- ----- -------- ----- -------- ...
  98  *             1 +   1 +     1 +   1 +      1 +   1 +
  99  *
 100  */
 101 static int
 102 cx_gamma_inc_F_CF(const double a, const complex double z, qpms_csf_result * result)
 103 {
 104   const int    nmax  =  5000;
 105   const double small =  DBL_EPSILON * DBL_EPSILON * DBL_EPSILON;
 106
 107   complex double hn = 1.0;           /* convergent */
 108   complex double Cn = 1.0 / small;
 109   complex double Dn = 1.0;
 110   int n;
 111
 112   /* n == 1 has a_1, b_1, b_0 independent of a,z,
 113      so that has been done by hand                */
 114   for ( n = 2 ; n < nmax ; n++ )
 115   {
 116     complex double an;
 117     complex double delta;
 118
 119     if(n % 2)
 120       an = 0.5*(n-1)/z;
 121     else
 122       an = (0.5*n-a)/z;
 123
 124     Dn = 1.0 + an * Dn;
 125     if ( cabs(Dn) < small )
 126       Dn = small;
 127     Cn = 1.0 + an/Cn;
 128     if ( cabs(Cn) < small )
 129       Cn = small;
 130     Dn = 1.0 / Dn;
 131     delta = Cn * Dn;
 132     hn *= delta;
 133     if(cabs(delta-1.0) < DBL_EPSILON) break;
 134   }
 135
 136   result->val = hn;
 137   result->err = 2.0*GSL_DBL_EPSILON * cabs(hn);
 138   result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * cabs(result->val);
 139
 140   if(n == nmax)
 141     GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER);
 142   else
 143     return GSL_SUCCESS;
 144 }
 145
 146 // Incomplete gamma fuction with complex second argument as continued fraction.
 147 int cx_gamma_inc_CF_e(const double a, const complex double z, qpms_csf_result *result)
 148 {
 149   qpms_csf_result F;
 150   gsl_sf_result pre;
 151   const complex double am1lgz = (a-1.0)*clog(z); // TODO check branches
 152   const int stat_F = cx_gamma_inc_F_CF(a, z, &F);
 153   const int stat_E = gsl_sf_exp_err_e(creal(am1lgz - z), GSL_DBL_EPSILON*cabs(am1lgz), &pre);
 154   complex double cpre = pre.val * cexp(I*cimag(am1lgz - z));// TODO add the error estimate for this
 155   //complex double cpre = cpow(z, a-1) * cexp(-z);
 156
 157
 158   result->val = F.val * cpre;
 159   result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err);
 160   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 161
 162   return GSL_ERROR_SELECT_2(stat_F, stat_E);
 163 }
 164
 165
 166 // Incomplete gamma function for complex second argument.
 167 int complex_gamma_inc_e(double a, complex double x, int m, qpms_csf_result *result) {
 168   int retval;
 169   if (creal(x) >= 0 &&
 170       (0 == fabs(cimag(x)) || // x is real positive; just use the real fun
 171       fabs(cimag(x)) < fabs(creal(x)) * COMPLEXPART_REL_ZERO_LIMIT)) {
 172     gsl_sf_result real_gamma_inc_result;
 173     retval = gsl_sf_gamma_inc_e(a, creal(x), &real_gamma_inc_result);
 174     result->val = real_gamma_inc_result.val;
 175     result->err = real_gamma_inc_result.err;
 176   } else if (creal(x) >= 0 && cabs(x) > 0.5)
 177     retval = cx_gamma_inc_CF_e(a, x, result);
 178   else if (QPMS_LIKELY(a > 0 || fmod(a, 1.0)))
 179     retval = cx_gamma_inc_series_e(a, x, result);
 180   else
 181   /* FIXME cx_gamma_inc_series_e() probably fails for non-positive integer a.
 182    * This does not matter for 2D lattices in 3D space,
 183    * but it might cause problems in the other cases.
 184    */
 185     QPMS_NOT_IMPLEMENTED("Incomplete Gamma function with non-positive integer a.");
 186   if (m) { // Non-principal branch.
 187     /* This might be sub-optimal, as Γ(a) has probably been already evaluated
 188      * somewhere in the functions called above. */
 189     gsl_sf_result fullgamma;
 190     int retval_fg = gsl_sf_gamma_e(a, &fullgamma);
 191     if (GSL_EUNDRFLW == retval_fg)
 192       fullgamma.err += DBL_MIN;
 193     else if (GSL_SUCCESS != retval_fg){
 194       result->val = NAN + NAN*I; result->err = NAN;
 195       return GSL_ERROR_SELECT_2(retval_fg, retval);
 196     }
 197     complex double f = cexp(2*m*M_PI*a*I);
 198     result->val *= f;
 199     f = -f + 1;
 200     result->err += cabs(f) * fullgamma.err;
 201     result->val += f * fullgamma.val;
 202   }
 203   return retval;
 204 }
 205
 206 // Exponential integral for complex argument; !UNTESTED! and probably not needed, as I expressed everything in terms of inc. gammas anyways.
 207 int complex_expint_n_e(int n, complex double x, qpms_csf_result *result) {
 208   if (creal(x) >= 0 &&
 209     (0 == fabs(cimag(x)) || // x is real positive; just use the real fun
 210     fabs(cimag(x)) < fabs(creal(x)) * COMPLEXPART_REL_ZERO_LIMIT)) {
 211     gsl_sf_result real_expint_result;
 212     int retval = gsl_sf_expint_En_e(n, creal(x), &real_expint_result);
 213     result->val = real_expint_result.val;
 214     result->err = real_expint_result.err;
 215     return retval;
 216   } else {
 217     int retval = complex_gamma_inc_e(-n+1, x, 0, result);
 218     complex double f = cpow(x, 2*n-2);
 219     result->val *= f;
 220     result->err *= cabs(f);
 221     return retval;
 222   }
 223 }
 224
 225 // inspired by GSL's hyperg_2F1_series
 226 int hyperg_2F2_series(const double a, const double b, const double c, const double d,
 227     const double x, gsl_sf_result *result
 228     )
 229 {
 230   double sum_pos = 1.0;
 231   double sum_neg = 0.0;
 232   double del_pos = 1.0;
 233   double del_neg = 0.0;
 234   double del = 1.0;
 235   double del_prev;
 236   double k = 0.0;
 237   int i = 0;
 238
 239   if(fabs(c) < GSL_DBL_EPSILON || fabs(d) < GSL_DBL_EPSILON) {
 240     result->val = NAN;
 241     result->err = INFINITY;
 242     GSL_ERROR ("error", GSL_EDOM);
 243   }
 244
 245   do {
 246     if(++i > 30000) {
 247       result->val  = sum_pos - sum_neg;
 248       result->err  = del_pos + del_neg;
 249       result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
 250       result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
 251       GSL_ERROR ("error", GSL_EMAXITER);
 252     }
 253     del_prev = del;
 254     del *= (a+k)*(b+k) * x / ((c+k) * (d+k) * (k+1.0));  /* Gauss series */
 255
 256     if(del > 0.0) {
 257       del_pos  =  del;
 258       sum_pos +=  del;
 259     }
 260     else if(del == 0.0) {
 261       /* Exact termination (a or b was a negative integer).
 262        */
 263       del_pos = 0.0;
 264       del_neg = 0.0;
 265       break;
 266     }
 267     else {
 268       del_neg  = -del;
 269       sum_neg -=  del;
 270     }
 271
 272     /*
 273      * This stopping criteria is taken from the thesis
 274      * "Computation of Hypergeometic Functions" by J. Pearson, pg. 31
 275      * (http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf)
 276      * and fixes bug #45926
 277      */
 278     if (fabs(del_prev / (sum_pos - sum_neg)) < GSL_DBL_EPSILON &&
 279         fabs(del / (sum_pos - sum_neg)) < GSL_DBL_EPSILON)
 280       break;
 281
 282     k += 1.0;
 283   } while(fabs((del_pos + del_neg)/(sum_pos-sum_neg)) > GSL_DBL_EPSILON);
 284
 285   result->val  = sum_pos - sum_neg;
 286   result->err  = del_pos + del_neg;
 287   result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
 288   result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);
 289
 290   return GSL_SUCCESS;
 291 }
 292
 293 // The Delta_n factor from [Kambe II], Appendix 3
 294 // \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f]
 295 void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, complex double x, complex double z) {
 296   complex double expfac = cexp(-x + 0.25 * z*z / x);
 297   complex double sqrtx = csqrt(x); // TODO check carefully, which branch is needed
 298   // These are used in the first two recurrences
 299   complex double w_plus  = Faddeeva_w(+z/(2*sqrtx) + I*sqrtx, 0);
 300   complex double w_minus = Faddeeva_w(-z/(2*sqrtx) + I*sqrtx, 0);
 301   QPMS_ASSERT(maxn >= 0);
 302   if (maxn >= 0) {
 303     target[0].val = 0.5 * M_SQRTPI * expfac * (w_minus + w_plus);
 304     target[0].err = NAN; //TODO
 305   }
 306   if (maxn >= 1) {
 307     target[1].val = I / z * M_SQRTPI * expfac * (w_minus - w_plus);
 308     target[1].err = NAN; //TODO
 309   }
 310   for(int n = 1; n < maxn; ++n) { // The rest via recurrence
 311     // TODO The cpow(x, 0.5 - n) might perhaps better be replaced with a recurrently computed variant
 312     target[n+1].val = -(4 / (z*z)) * (-(0.5 - n) * target[n].val + target[n-1].val - cpow(x, 0.5 - n) * expfac);
 313     target[n+1].err = NAN; //TODO
 314   }
 315 }
 316
 317