changeset 17:cb249ba61e9e

Moved from parent directory.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/general/numerical/scalar/hzeta.c	Thu Jan 17 13:21:46 2013 +0000
@@ -0,0 +1,96 @@
+/*=================================================================
+ * hzeta.c 
+ *
+ * Wrapper function for gsl_sf_hzeta_e
+ *
+ * This MEX function computes the Hurwitz-zeta-function by calling 
+ * the function gsl_sf_hzeta_e from the GNU Scientific Library. For 
+ * x>1 and a>=0 the Hurwitz-zeta-function is defined as
+ * 
+ *                        oo 
+ *                       -----
+ *                        \         1
+ *       hzeta(x, a) =     )     -------- .
+ *                        /             x
+ *                       -----   (i + a)
+ *                       i = 0
+ *
+ * The MEX function takes two arguments, x and a. If both arguments 
+ * are matrices of the same size, then the MEX function will compute
+ * the Hurwitz-zeta-function elementwise and the result of the MEX 
+ * function will be a matrix of the same size as the inputs. If one 
+ * input argument is a scalar s, then the function behaves as if this
+ * argument would be a matrix filled with the value s of the same 
+ * size as the other argument.
+ *
+ * Compile this MEX function with "mex hzeta.c -lgsl". See Matlab
+ * documentation for additional information on compiling MEX functions.
+ *
+ * (C) by Heiko Bauke, 2007, [e-mail]
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License in
+ * version 2 as published by the Free Software Foundation.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+ * 02111-1307, USA.
+ *
+ *=================================================================*/
+
+#include "mex.h"
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_sf_zeta.h>
+
+void mexFunction(int nlhs, mxArray *plhs[], 
+                 int nrhs, const mxArray *prhs[]) {
+  double *x, *y, *a;
+  int nrows_x, ncols_x, i_x, 
+      nrows_a, ncols_a, i_a, 
+      nrows_y, ncols_y, i_y, 
+      status;
+  gsl_sf_result result;
+
+  /* check input arguments */
+  if (nrhs!=2)
+    mexErrMsgTxt("Two input arguments required.\n");
+  if (!mxIsDouble(prhs[0]))
+    mexErrMsgTxt("1st argument not a matrix of double."); 
+  nrows_x=mxGetM(prhs[0]);
+  ncols_x=mxGetN(prhs[0]);
+  if (!mxIsDouble(prhs[1]))
+    mexErrMsgTxt("2nd argument not a matrix of double."); 
+  nrows_a=mxGetM(prhs[1]);
+  ncols_a=mxGetN(prhs[1]);
+  if (!( (nrows_a==nrows_x && ncols_a==ncols_x) ||
+	 (nrows_x==1 && ncols_x==1) ||
+	 (nrows_a==1 && ncols_a==1) ) ) 
+    mexErrMsgTxt("arguments must be matrices of the same size or scalars."); 
+  /* create matrix for the return argument */ 
+  nrows_y=nrows_x>nrows_a ? nrows_x : nrows_a;
+  ncols_y=ncols_x>ncols_a ? ncols_x : ncols_a;
+  plhs[0]=mxCreateDoubleMatrix(nrows_y, ncols_y, mxREAL); 
+  /* assign pointers to each input and output */ 
+  x=mxGetPr(prhs[0]); 
+  a=mxGetPr(prhs[1]);
+  y=mxGetPr(plhs[0]); 
+  /* compute zeta function */
+  gsl_set_error_handler_off();
+  for (i_x=i_a=i_y=0; i_y<nrows_y*ncols_y; ++i_y) {
+    if (gsl_sf_hzeta_e(x[i_x], a[i_a], &result)==GSL_SUCCESS)
+      y[i_y]=result.val;
+    else
+      y[i_y]=GSL_NAN;
+    if (!(nrows_x==1 && ncols_x==1))
+      ++i_x;
+    if (!(nrows_a==1 && ncols_a==1))
+      ++i_a;
+  }
+}

Binary file general/numerical/scalar/hzeta.mexmac has changed

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/general/numerical/scalar/safe_exp.m	Thu Jan 17 13:21:46 2013 +0000
@@ -0,0 +1,11 @@
+function [y,k]=safe_exp(x,l)
+% safe_exp - high dynamic range exponential
+%
+% safe_exp :: X:[[N,M]] -> [[N,M]], [[1,M]].
+%
+% returns y and k such that exp(X) = y * exp(k) and
+% maximum value in y is 1.
+
+k=max(x,[],1);
+y=exp(x-repmat(k,size(x,1),1));
+	

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/general/numerical/scalar/safe_exp2.m	Thu Jan 17 13:21:46 2013 +0000
@@ -0,0 +1,11 @@
+function [y,k]=safe_exp2(x,l)
+% safe_exp2 - high dynamic range exponential
+%
+% safe_exp2 :: X:[[N,M]], L:real -> [[N,M]], [[1,M]].
+%
+% returns y and k such that exp(X) = y * exp(k) and
+% maximum value in y is exp(L).
+% Maximum sensible value of L is log(realmax)=709.78.
+
+k=max(x,[],1)-l; 
+y=exp(x-repmat(k,size(x,1),1));