git.gag.com Git - debian/atlc/blob - src/non_gui/gsl_sf_ellint_Kcomp.c

   1
   2
   3 #include "config.h"
   4
   5 #ifdef HAVE_MATH_H
   6 #include <math.h>
   7 #endif
   8
   9 #include "definitions.h"
  10 #include "gsl_types.h"
  11 #include "gsl_definitions.h"
  12 #include "exit_codes.h"
  13
  14 /* This part of atlc is a copy from the code in the GPL'ed
  15 GNU Scientific Library, gsl. By including this code, it saves
  16 the user having to like with gsl. */
  17
  18 const double gsl_prec_eps[3];
  19
  20
  21
  22 double gsl_sf_ellint_Kcomp(double k, gsl_mode_t mode)
  23 {
  24   gsl_sf_result result;
  25   gsl_sf_ellint_Kcomp_e(k, mode, &result);
  26   return result.val;
  27 }
  28
  29
  30 /* [Carlson, Numer. Math. 33 (1979) 1, (4.5)] */
  31 int gsl_sf_ellint_Kcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
  32 {
  33   int return_val=0;
  34   if(k*k >= 1.0) {
  35     exit_with_msg_and_exit_code("domain error in gsl_sf_Kcomp_e", DOMAIN_ERROR);
  36   }
  37   else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
  38     /* [Abramowitz+Stegun, 17.3.33] */
  39     const double y = 1.0 - k*k;
  40     const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };
  41     const double b[] = { 0.5, 0.12498593597, 0.06880248576 };
  42     const double ta = a[0] + y*(a[1] + y*a[2]);
  43     const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));
  44     result->val = ta + tb;
  45     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  46     return_val=GSL_SUCCESS;
  47   }
  48   else {
  49     /* This was previously computed as,
  50
  51          return gsl_sf_ellint_RF_e(0.0, 1.0 - k*k, 1.0, mode, result);
  52
  53        but this underestimated the total error for small k, since the
  54        argument y=1-k^2 is not exact (there is an absolute error of
  55        GSL_DBL_EPSILON near y=0 due to cancellation in the subtraction).
  56        Taking the singular behavior of -log(y) above gives an error
  57        of 0.5*epsilon/y near y=0. (BJG) */
  58
  59     double y = 1.0 - k*k;
  60     int status = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, result);
  61     result->err += 0.5 * GSL_DBL_EPSILON / y;
  62     return_val=status ;
  63   }
  64   return(return_val);
  65 }
  66
  67 int gsl_sf_ellint_RF_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
  68 {
  69   int return_val=0;
  70   const double lolim = 5.0 * GSL_DBL_MIN;
  71   const double uplim = 0.2 * GSL_DBL_MAX;
  72   const gsl_prec_t goal = GSL_MODE_PREC(mode);
  73   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  74   const double prec   = gsl_prec_eps[goal];
  75
  76   if(x < 0.0 || y < 0.0 || z < 0.0) {
  77     exit_with_msg_and_exit_code("domain error in gsl_sf_ellint_RF_e", DOMAIN_ERROR);
  78   }
  79   else if(x+y < lolim || x+z < lolim || y+z < lolim) {
  80     exit_with_msg_and_exit_code("domain error in gsl_sf_ellint_RF_e", DOMAIN_ERROR);
  81   }
  82   else if(locMAX3(x,y,z) < uplim) {
  83     const double c1 = 1.0 / 24.0;
  84     const double c2 = 3.0 / 44.0;
  85     const double c3 = 1.0 / 14.0;
  86     double xn = x;
  87     double yn = y;
  88     double zn = z;
  89     double mu, xndev, yndev, zndev, e2, e3, s;
  90     while(1) {
  91       double epslon, lamda;
  92       double xnroot, ynroot, znroot;
  93       mu = (xn + yn + zn) / 3.0;
  94       xndev = 2.0 - (mu + xn) / mu;
  95       yndev = 2.0 - (mu + yn) / mu;
  96       zndev = 2.0 - (mu + zn) / mu;
  97       epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));
  98       if (epslon < errtol) break;
  99       xnroot = sqrt(xn);
 100       ynroot = sqrt(yn);
 101       znroot = sqrt(zn);
 102       lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
 103       xn = (xn + lamda) * 0.25;
 104       yn = (yn + lamda) * 0.25;
 105       zn = (zn + lamda) * 0.25;
 106     }
 107     e2 = xndev * yndev - zndev * zndev;
 108     e3 = xndev * yndev * zndev;
 109     s = 1.0 + (c1 * e2 - 0.1 - c2 * e3) * e2 + c3 * e3;
 110     result->val = s / sqrt(mu);
 111     result->err = prec * fabs(result->val);
 112     return_val= GSL_SUCCESS;
 113   }
 114   else {
 115     exit_with_msg_and_exit_code("domain error in gsl_sf_ellint_RF_e", DOMAIN_ERROR);
 116   }
 117   return(return_val);
 118 }
 119
 120
 121 /* static double locMAX3(double x, double y, double z) */
 122 double locMAX3(double x, double y, double z)
 123 {
 124   double xy = GSL_MAX(x, y);
 125   return GSL_MAX(xy, z);
 126 }
 127
 128
 129 #define EVAL_RESULT(fn) \
 130    gsl_sf_result result; \
 131    int status = fn; \
 132    if (status != GSL_SUCCESS) { \
 133      GSL_ERROR_VAL(#fn, status, result.val); \
 134    } ; \
 135 return result.val;
 136