[Help-gsl] Different value for mathieu_ce in Mathematica and GSL

[Phyks]

Hi,

I have some code that I prototyped in Mathematica and am now writing in 
C using GSL, that makes use of Mathieu functions. I have different 
results between the two of them, and I cannot figure out whether this is 
a bug in GSL, Mathematica or simply some misunderstanding from my part.

I am using `MathieuC` function in latest Mathematica 
(http://reference.wolfram.com/language/ref/MathieuC.html) which should 
be the same function as `gsl_sf_mathieu_ce` 
(https://www.gnu.org/software/gsl/manual/html_node/Angular-Mathieu-Functions.html#Angular-Mathieu-Functions) 
except that the former one takes a single `a` argument being the 
characteristic value whereas the GSL 2.3 implementation takes the order 
`n` and the `q` parameter directly.

So, I guess,
```
> N[MathieuC[MathieuCharacteristicA[0, -1], -1, 2*Pi/180]]
1.41071
```

should be equivalent to
```
gsl_sf_mathieu_ce(0, -1.0, 2.0 * M_PI / 180.0)
```
which gives a totally different value: 0.99751942347886335.


I tried to debug with different values, and the discrepancies between 
Mathematica and GSL seems to appear only when the `q` parameter (-1.0 
here) is negative. If I take 1.0 instead, I get values in agreement. I 
tried to find yet another implementation to debug it, and found Scipy 
(https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.mathieu_cem.html#scipy.special.mathieu_cem) 
which relies on Fortran SPECFUN library apparently, and is in agreement 
with GSL.

I am missing something? Thanks!
--
Phyks