[Bug-gsl] bug report on about the Mathieu function

[Takeshi Morita]

Dear Sirs/Madams,

My name is Takeshi Morita, a Japanese physics researcher.

I found a bug in GSL about the Mathieu characteristic function and Mathieu 
function.

I needed to calculate the Mathieu characteristic function;
gsl_sf_mathieu_a(r, q, &result)
gsl_sf_mathieu_b(r, q, &result)
at large q. 

Then I found that
gsl_sf_mathieu_a(r, q, &result) - gsl_sf_mathieu_b(r, q, &result)
and/or
gsl_sf_mathieu_b(r+1, q, &result) - gsl_sf_mathieu_a(r, q, &result)
becomes NEGATIVE, although they have to be always POSITIVE.
(See p.724 (a) in the Handbook of Mathematical Functions by Abramowitz and 
Stegun.)

The attached file is a sample code to see this bug.

For example, when q=2000.0, 
gsl_sf_mathieu_b(70, q, &result) - gsl_sf_mathieu_a(69, q, &result)
becomes -355.55....

I tried several q and r, and this error happenes when q is roughly larger than
1024 and r is around 70.


Related to this issue, gsl_sf_mathieu_ce and  gsl_sf_mathieu_se
do not work properly.



I compiled the code by using gcc (version 4.5.3) on Cygwin on Win 7.
I asked my friend and he also confirmed the same bug on his following PC 
environment;
Linux 3.2.0-4-amd64 #1 SMP Debian 3.2.32-1 x86_64
gcc version 4.7.2 (Debian 4.7.2-5)
gsl version 1.15

So I think this bug is independent of the PC environment.




I would appreciate it if you could confirm the bug and fix the problem.

Sincerely yours,


Takeshi Morita
-- 
=================================================================
Takeshi Morita

KEK Theory Center

Phone  : +81-29-864-5401 (Lab)
Fax    : +81-29-864-5755  (Lab)
E-mail : tmor...@post.kek.jp
=================================================================


#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_mathieu.h>

int main(void) {

  gsl_sf_result result;

  int n = 100;
  int n1 = 10000;
  int i;

  double q = 2000.0;

  double diff;

  gsl_sf_mathieu_workspace * workspace =  gsl_sf_mathieu_alloc(n1, q);

  gsl_sf_result old, new;

  gsl_sf_mathieu_a (0, q, &old);

  for (i = 1; i < n; ++i) {
    gsl_sf_mathieu_b (i, q, &new);

    diff = new.val - old.val;
    if (diff <= 0) {
      printf("b[%d] - a[%d-1] , %f\n", i, i, diff);
    }
    old = new;

    gsl_sf_mathieu_a (i, q, &new);

    diff = new.val - old.val;
    if (diff <= 0) {
      printf("a[%d] - b[%d], %f\n", i, i, diff);
    }
    old = new;

  }
  gsl_sf_mathieu_free(workspace);

  return 0;
}