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