Improved Mathieu equation eigenvalues and testing

[Lowell Johnson <ldj00 at sio dot midco dot net>]

Hi,

I've attached a patch for the Mathieu function eigenvalue accuracy and 
testing in the specfunc directory.  Based on input from Brian Gladman, 
I've increased the number of terms in the eigenvalue matrix to ensure 
double precision accuracy.

Precision comparison values with a good range of orders and arguments 
are still needed for improved Mathieu function tests.  Brian (Gladman), 
if you could send me some values, that would be great.  That goes for 
anyone else who has data to offer also.  Note that the values needed 
for test comparisons are for the actual Mathieu functions, not the 
characteristic values.

Regards,

    Lowell
-- 
Lowell D. Johnson

Index: gsl_sf_mathieu.h
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/gsl_sf_mathieu.h,v
retrieving revision 1.1
diff -u -r1.1 gsl_sf_mathieu.h
--- gsl_sf_mathieu.h	18 Apr 2006 17:59:46 -0000	1.1
+++ gsl_sf_mathieu.h	11 Jun 2006 20:05:15 -0000
@@ -57,7 +57,8 @@
 
 
 /* Allocate computational storage space for eigenvalue solution. */
-gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn);
+gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn,
+                                               const double qq);
 void gsl_sf_mathieu_free(gsl_sf_mathieu_workspace *workspace);
 
 /* Compute an angular Mathieu function. */
Index: mathieu_workspace.c
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/mathieu_workspace.c,v
retrieving revision 1.1
diff -u -r1.1 mathieu_workspace.c
--- mathieu_workspace.c	18 Apr 2006 17:59:46 -0000	1.1
+++ mathieu_workspace.c	11 Jun 2006 20:05:15 -0000
@@ -25,12 +25,16 @@
 #include "mathieu.h"
 
 
-gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn)
+gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn,
+                                               const double qq)
 {
   gsl_sf_mathieu_workspace *workspace;
   unsigned int even_order = nn/2 + 1, odd_order = (nn + 1)/2,
-               extra_values = 5;
+               extra_values;
 
+  /* Compute the maximum number of extra terms required for 10^-18 root
+     accuracy for a given value of q (contributed by Brian Gladman). */
+  extra_values = (int)(2.1*pow(qq, 0.37)) + 9;
   
   if (nn == 0)
   {
Index: test_mathieu.c
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/test_mathieu.c,v
retrieving revision 1.1
diff -u -r1.1 test_mathieu.c
--- test_mathieu.c	18 Apr 2006 17:59:46 -0000	1.1
+++ test_mathieu.c	11 Jun 2006 20:05:15 -0000
@@ -24,15 +24,17 @@
 #include "test_sf.h"
 #include "mathieu.h"
 
-/* static double ce[100]; */
-/* static double co[100]; */
-/* static double se[100]; */
-/* static double so[100]; */
+#define NVAL 100
+
+static double c[NVAL];
+static double s[NVAL];
 
 int test_mathieu(void)
 {
   gsl_sf_result r;
+  gsl_sf_mathieu_workspace *work;
   int s = 0;
+  int sa;
 
   TEST_SF(s, gsl_sf_mathieu_c, (0, 0.0, 0.0, &r),
           0.7071067811865475, TEST_SNGL, GSL_SUCCESS);
@@ -179,5 +181,34 @@
   TEST_SF(s, gsl_sf_mathieu_s, (15, 25.0, M_PI_2, &r),
           -0.9467086958780897, TEST_SNGL, GSL_SUCCESS);
 
+  work = gsl_sf_mathieu_alloc(NVAL, 20.0);
+  sa = 0;
+  gsl_sf_mathieu_c_array(0, 5, 0.0, M_PI_2, work, c);
+  sa += (test_sf_frac_diff(c[0], 0.7071067811865475) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[2], -1.0) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[4], 1.0) > TEST_SNGL);
+  gsl_test(sa, "gsl_sf_mathieu_c_array");
+  s += sa;
+
+  sa = 0;
+  gsl_sf_mathieu_c_array(0, 15, 20.0, 0.0, work, c);
+  sa += (test_sf_frac_diff(c[0], 0.0006037438292242197) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[1], 0.005051813764712904) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[2], 0.02864894314707431) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[5], 0.9365755314226215) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[10], 1.117788631259397) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[15], 1.047084344162887) > TEST_SNGL);
+  gsl_test(sa, "gsl_sf_mathieu_c_array");
+  s += sa;
+
+  sa = 0;
+  gsl_sf_mathieu_s_array(1, 15, 20.0, M_PI_2, work, c);
+  sa += (test_sf_frac_diff(c[0], 1.609891592603772) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[4], 0.8635431218533667) > TEST_SNGL);
+  sa += (test_sf_frac_diff(c[14], -0.9570452540612817) > TEST_SNGL);
+  gsl_test(sa, "gsl_sf_mathieu_s_array");
+  s += sa;
+
+  free(work);
   return s;
 }