octave-gsl/elliptic-patch

   1 diff --git a/INDEX b/INDEX
   2 index 542b535..defc704 100644
   3 --- a/INDEX
   4 +++ b/INDEX
   5 @@ -52,6 +52,8 @@ Special functions
   6   debye_4
   7   ellint_Ecomp
   8   ellint_Kcomp
   9 + ellint_F
  10 + ellint_E
  11   erf_Q
  12   erf_Z
  13   erf_gsl
  14 diff --git a/PKG_ADD b/PKG_ADD
  15 index d11ee4a..af19737 100644
  16 --- a/PKG_ADD
  17 +++ b/PKG_ADD
  18 @@ -45,6 +45,8 @@ autoload ("debye_3", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct")
  19  autoload ("debye_4", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  20  autoload ("ellint_Ecomp", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  21  autoload ("ellint_Kcomp", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  22 +autoload ("ellint_F", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  23 +autoload ("ellint_E", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  24  autoload ("erfc_gsl", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  25  autoload ("erf_gsl", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  26  autoload ("erf_Q", fullfile (fileparts (mfilename ("fullpath")), "gsl_sf.oct"));
  27 diff --git a/src/README b/src/README
  28 index e93b138..14506f1 100644
  29 --- a/src/README
  30 +++ b/src/README
  31 @@ -43,6 +43,7 @@ As of this writing there are:
  32         double fn(int,double)
  33         double fn(int,int,double)
  34         double fn(double,mode)
  35 +       double fn(double,double,mode)
  36
  37  mode is a precision mode accepting a gsl_mode_t value.
  38
  39 diff --git a/src/buildgsl_sf.sh b/src/buildgsl_sf.sh
  40 index 5babc57..b6a93d9 100644
  41 --- a/src/buildgsl_sf.sh
  42 +++ b/src/buildgsl_sf.sh
  43 @@ -1076,7 +1076,23 @@ parameter @math{m = k^2}.
  44  EOF
  45  ./replace_template.sh double_mode_to_double.cc.template >> gsl_sf.cc
  46
  47 +# double double mode to double
  48
  49 +export octave_name=ellint_F
  50 +export    funcname=gsl_sf_ellint_F
  51 +cat << \EOF > docstring.txt
  52 +This routine computes the incomplete elliptic integral F(\phi,k) to
  53 +the accuracy specified by the mode variable mode.
  54 +EOF
  55 +./replace_template.sh double_double_mode_to_double.cc.template >> gsl_sf.cc
  56 +
  57 +export octave_name=ellint_E
  58 +export    funcname=gsl_sf_ellint_E
  59 +cat << \EOF > docstring.txt
  60 +This routine computes the incomplete elliptic integral E(\phi,k) to
  61 +the accuracy specified by the mode variable mode.
  62 +EOF
  63 +./replace_template.sh double_double_mode_to_double.cc.template >> gsl_sf.cc
  64
  65  # int to double
  66
  67 diff --git a/src/double_double_mode_to_double.cc.template b/src/double_double_mode_to_double.cc.template
  68 new file mode 100644
  69 index 0000000..e38a596
  70 --- /dev/null
  71 +++ b/src/double_double_mode_to_double.cc.template
  72 @@ -0,0 +1,147 @@
  73 +DEFUN_DLD(GSL_OCTAVE_NAME, args, nargout, "\
  74 +  -*- texinfo -*-\n\
  75 +@deftypefn {Loadable Function} {@var{z} =} GSL_OCTAVE_NAME (@var{x}, @var{y}, @var{mode})\n\
  76 +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} GSL_OCTAVE_NAME (@dots{})\n\
  77 +\n\
  78 +GSL_FUNC_DOCSTRING
  79 +\n\
  80 +The third argument @var{mode} must be an integer corresponding to\n\
  81 +\n\
  82 +@table @asis\n\
  83 +@item 0 = GSL_PREC_DOUBLE\n\
  84 +    Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\
  85 +@item 1 = GSL_PREC_SINGLE\n\
  86 +    Single-precision, a relative accuracy of approximately @code{10^-7}.\n\
  87 +@item 2 = GSL_PREC_APPROX\n\
  88 +    Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\
  89 +@end table\n\
  90 +\n\
  91 +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\
  92 +\n\
  93 +This function is from the GNU Scientific Library,\n\
  94 +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\
  95 +@end deftypefn\n\
  96 +")
  97 +{
  98 +  int i;
  99 +  dim_vector dv;
 100 +
 101 +  gsl_set_error_handler (octave_gsl_errorhandler);
 102 +
 103 +  if(!(args.length() == 3 || args.length() == 2)) {
 104 +    print_usage ();
 105 +    return octave_value();
 106 +  }
 107 +  if(!args(0).is_real_type() || !args(1).is_real_type() || (args.length() == 3 && !args(2).is_real_type())) {
 108 +    error("The arguments must be real.");
 109 +    print_usage ();
 110 +    return octave_value();
 111 +  }
 112 +  if(args.length() == 3 && !args(2).is_scalar_type()) {
 113 +    error("The mode must be scalar.");
 114 +    print_usage ();
 115 +    return octave_value();
 116 +  }
 117 +
 118 +  int mode = args.length() == 3 ? static_cast<int>((args(2).array_value())(0)) : 0;
 119 +  if(mode < 0)
 120 +    mode = 0;
 121 +  else if(mode > 2)
 122 +    mode = 2;
 123 +
 124 +  // Nice combinatorial explosion here
 125 +  NDArray n = args(0).array_value();
 126 +  NDArray x = args(1).array_value();
 127 +  if(n.length() == x.length()) {
 128 +    dv = x.dims();
 129 +    NDArray y(dv);
 130 +    int len = x.length();
 131 +    if(nargout < 2) {
 132 +      for(i = 0; i < len; i++) {
 133 +        y.xelem(i) = GSL_FUNC_NAME (n.xelem(i),
 134 +                                    x.xelem(i),
 135 +                                    mode);
 136 +      }
 137 +      return octave_value(y);
 138 +    } else {
 139 +      NDArray err(dv);
 140 +      gsl_sf_result result;
 141 +      octave_value_list retval;
 142 +      for(i = 0; i < len; i++) {
 143 +        GSL_FUNC_NAME_e (n.xelem(i),
 144 +                         x.xelem(i),
 145 +                         mode,
 146 +                         &result);
 147 +        y.xelem(i) = result.val;
 148 +        err.xelem(i) = result.err;
 149 +      }
 150 +      retval(1) = octave_value(err);
 151 +      retval(0) = octave_value(y);
 152 +      return retval;
 153 +    }
 154 +  } else if(n.length() == 1) {
 155 +    dv = x.dims();
 156 +    NDArray y(dv);
 157 +    int len = x.length();
 158 +    double ndouble = n.xelem(0);
 159 +    if(nargout < 2) {
 160 +      for(i = 0; i < len; i++) {
 161 +        y.xelem(i) = GSL_FUNC_NAME (ndouble,
 162 +                                    x.xelem(i),
 163 +                                    mode);
 164 +      }
 165 +      return octave_value(y);
 166 +    } else {
 167 +      NDArray err(dv);
 168 +      gsl_sf_result result;
 169 +      octave_value_list retval;
 170 +      for(i = 0; i < len; i++) {
 171 +        GSL_FUNC_NAME_e (ndouble, x.xelem(i), mode, &result);
 172 +        y.xelem(i) = result.val;
 173 +        err.xelem(i) = result.err;
 174 +      }
 175 +      retval(1) = octave_value(err);
 176 +      retval(0) = octave_value(y);
 177 +      return retval;
 178 +    }
 179 +  } else if(x.length() == 1) {
 180 +    dv = n.dims();
 181 +    NDArray y(dv);
 182 +    int len = n.length();
 183 +    double xdouble = x.xelem(0);
 184 +    if(nargout < 2) {
 185 +      for(i = 0; i < len; i++) {
 186 +        y.xelem(i) = GSL_FUNC_NAME (n.xelem(i),
 187 +                                    xdouble,
 188 +                                    mode);
 189 +      }
 190 +      return octave_value(y);
 191 +    } else {
 192 +      NDArray err(dv);
 193 +      gsl_sf_result result;
 194 +      octave_value_list retval;
 195 +      for(i = 0; i < len; i++) {
 196 +        GSL_FUNC_NAME_e (n.xelem(i),
 197 +                         xdouble,
 198 +                         mode,
 199 +                         &result);
 200 +        y.xelem(i) = result.val;
 201 +        err.xelem(i) = result.err;
 202 +      }
 203 +      retval(1) = octave_value(err);
 204 +      retval(0) = octave_value(y);
 205 +      return retval;
 206 +    }
 207 +  } else {
 208 +    error("First and second argument must either have the same size, or one of them must be scalar.");
 209 +    print_usage ();
 210 +  }
 211 +
 212 +  return octave_value();
 213 +}
 214 +
 215 +/*
 216 +  ;;; Local Variables: ***
 217 +  ;;; mode: C++ ***
 218 +  ;;; End: ***
 219 +*/