This file is the GSL bug tracking system.  The CVS version of this
file should be kept up-to-date.

----------------------------------------------------------------------
BUG#1 -- gsl_sf_hyperg_2F1_e fails for some arguments 

From: keith.briggs@bt.com
To: gsl-discuss@sources.redhat.com
Subject: gsl_sf_hyperg_2F1 bug report
Date: Thu, 31 Jan 2002 12:30:04 -0000

gsl_sf_hyperg_2F1_e fails with arguments (1,13,14,0.999227196008978,&r).
It should return 53.4645... .

#include <gsl/gsl_sf.h>
#include <stdio.h>

int main (void)
{
  gsl_sf_result r;
  gsl_sf_hyperg_2F1_e (1,13,14,0.999227196008978,&r);
  printf("r = %g %g\n", r.val, r.err);
}

NOTES: The program overflows the maximum number of iterations in
gsl_sf_hyperg_2F1, due to the presence of a nearby singularity at
(c=a+b,x=1) so the sum is slowly convergent.

The exact result is 53.46451441879150950530608621 as calculated by
gp-pari using sumpos(k=0,gamma(a+k)*gamma(b+k)*gamma(c)*gamma(1)/
(gamma(c+k)*gamma(1+k)*gamma(a)*gamma(b))*x^k)

The code needs to be extended to handle the case c=a+b. This is the
main problem. The case c=a+b is special and needs to be computed
differently.  There is a special formula given for it in Abramowitz &
Stegun 15.3.10

----------------------------------------------------------------------
BUG#3 -- gsl_sf_bessel_Jn fails for large argument 

From: Rodolphe Conan <rconan@eso.org>
To: gsl list <gsl-discuss@sources.redhat.com>
Subject: Re: Bessel function bug or limitation?
Date: Tue, 09 Apr 2002 09:52:19 +0100

Rodolphe Conan writes:
 >   cout << "besselj(2,15000)="  << gsl_sf_bessel_Jn(2,15000) << endl;
 > and obtained:
 > gsl: bessel.c:531: ERROR: error
 > Abort (core dumped)

One of the continued fraction computations is breaking down

----------------------------------------------------------------------
BUG#4 -- gsl_linalg_solve_symm_cyc_tridiag() crashes 

From: David Necas (Yeti) <yeti@physics.muni.cz>
To: Brian Gough <bjg@network-theory.co.uk>
Subject: Re: gsl_linalg_solve_symm_cyc_tridiag problem
Date: Sun, 14 Apr 2002 18:02:21 +0200

The function gsl_linalg_solve_symm_cyc_tridiag() crashes when run with
N = 1 (and doesn't make sense even with N = 2) -- it should return
some error code instead (the same applies to its non-cyclic
counterpair).

Both cyclic and non-cyclic solvers return zero error code when feed
with singular matrices, or fail for another reason (so the results are
only bunches of NaN's). They should return some error code instead.

----------------------------------------------------------------------

BUG#5 -- broken error terms for implicit odes 

The error terms for the implicit ode integrators are broken.  They
have 'FIXME' entries still in the code.  Only the bsimp error term is
implemented.

----------------------------------------------------------------------

BUG#7 -- gsl_sf_coupling_3j loses accuracy for large arguments 

From: Filip Floegel <floegel@itp.uni-hannover.de>
Subject: gsl_sf_coupling_3j
Date: Wed, 2 Oct 2002 17:45:52 +0200

gsl_sf_coupling_3j doesn't work properly for large
angular mementa like j1=j2=j3=40 m1=m2=m3=0 !!!

gsl_sf_coupling_3j_e(80,80,80,0,0,0,&result)

which correspond to j1=j2=j3=40 m1=m2=m3=0

gives me:

      result.val = 0.0625
      result.err = 0.172766

mathematica gives me:

      N[ThreeJSymbol[{40,0},{40,0},{40,0}]] = 0.0149685

my guess is that gamma.c which is used for calculating fractional
numbers probably diverges for large integers. 

Note from Serge Winitzki <serge@cosmos.phy.tufts.edu>:

The package "matpack" (www.matpack.de) includes many special functions, 
also the 3j symbols. They refer to some quite complicated numerical 
methods using recursion relations to get the right answers for large 
momenta, and to 1975-1976 papers by Schulten and Gordon for the 
description of the algorithms. The papers can be downloaded for free at
http://www.ks.uiuc.edu/Publications/Papers/ 


----------------------------------------------------------------------

BUG#8 -- inexact coefficients in rk8pd.c 

From: Luc Maisonobe <Luc.Maisonobe@c-s.fr>
To: gsl-discuss@sources.redhat.com
Subject: further thoughts about Dormand-Prince 8 (RK8PD)
Date: Wed, 14 Aug 2002 10:50:49 +0200

I was looking for some references concerning Runge-Kutta methods when I
noticed GSL had an high order one. I also found a question in the list
archive (April 2002) about the references of this method which is
implemented in rk8pd.c. It was said the coefficients were taken from the
"Numerical Algorithms with C" book by Engeln-Mullges and Uhlig.

I have checked the coefficients somewhat with a little java tool I have
developped (see http://www.spaceroots.org/archive.htm#RKCheckSoftware)
and found they were not exact. I think this method is really the method
that is already in rksuite (http://www.netlib.org/ode/rksuite/) were the
coefficients are given as real values with 30 decimal digits. The
coefficients have probably been approximated as fractions later on.
However, these approximations are not perfect, they are good only for
the first 16 or 18 digits depending on the coefficient.

This has no consequence for practical purposes since they are stored in
double variables, but give a false impression of beeing exact
expressions. Well, there are even some coefficients that should really
be rational numbers but for which wrong numerators and denominators are
given. As an example, the first and fourth elements of the b7 array are
given as 29443841.0 / 614563906.0 and 77736538.0 / 692538347, hence the
sum off all elements of the b7 array (which should theoretically be
equal to ah[5]) only approximate this. For these two coefficients, this
could have been avoided using  215595617.0 / 4500000000.0 and
202047683.0 / 1800000000.0, which also looks more coherent with the
other coefficients.

The rksuite comments say this method is described in this paper :

   High Order Embedded Runge-Kutta Formulae
   P.J. Prince and J.R. Dormand
   J. Comp. Appl. Math.,7, pp. 67-75, 1981

It also says the method is an 8(7) method (i.e. the coefficients set
used to advance integration is order 8 and error estimation is order 7).
If I use my tool to check the order, I am forced to check the order
conditions numerically with a tolerance since I do not have an exact
expression of the coefficients. Since even if some conditions are not
mathematically met, the residuals are small and could be below the
tolerance. There are tolerance values for which such numerical test
dedeuce the method is of order 9, as is said in GSL. However, I am not
convinced, there are to few parameters for the large number of order
conditions needed at order 9.

I would suggest to correct the coefficients in rk8pd.c (just put the
literal constants of rksuite) and to add the reference to the article.

----------------------------------------------------------------------

BUG#9 -- problem with multimin example 

From: Thomas Kunert <kunert@physik.tu-dresden.de>
To: gsl-discuss@sources.redhat.com
Subject: Re: multimin example

 The example for multimin doesn't give me the expected output but a much
 worse result. Any idea about the reason? I have gsl 1.2 on Mandrake 8.2.

     1 4.99629 6.99072  687.84780
     2 4.98886 6.97215  683.55456
     3 4.97400 6.93501  675.01278
     4 4.94429 6.86073  658.10798

It's true that there is a problem there, but I don't know why.

----------------------------------------------------------------------
BUG#10 -- gsl_sf_fermi_dirac_int error estimates 

Some of the error estimates on gsl_sf_fermi_dirac_int are much too
large.  The value itself is pretty accurate.

In the test_sf_result you need to work in something like

  if(r.err > 1.0e5 * (fabs(r.val - val) + GSL_DBL_EPSILON * fabs(val)))
        s |= TEST_SF_INCONS;   /* error estimate too large */

in addition to the existing

   if(fabs(val - r.val) > 2.0*r.err) s |= TEST_SF_INCONS;

to catch those.  

#include<stdio.h>
#include<gsl/gsl_errno.h>
#include<gsl/gsl_sf.h>

int main() {
  gsl_sf_result r;
  int status;

  status = gsl_sf_fermi_dirac_int_e (9, 500.0, &r);
  printf("FD_9(500) = %.18e +- %.18e\n", r.val, r.err);
}

result,

FD_9(500) = 2.692738498426942915e+20 +- 2.629627439870118259e+47

Exact     = 2.692738498...e+20

----------------------------------------------------------------------

BUG#11 -- SVD does not handle N=1 

From: Tiago de Paula Peixoto <count0@fissionauthority.com>
To: gsl-discuss@sources.redhat.com
Subject: Problem with least-squares fitting.
Date: Mon, 4 Nov 2002 14:46:02 -0200

 I'm experiencing a weird error message using the GSL linear least squares
multi-parameter fitting. When I try to fit a function with only one
parameter, such as "y = a * x", I get the following error:

subvector_source.c:28: vector length n must be positive integer
(gsl_error: invalid argument supplied by user)

The fit results seem OK, and I don't know what I'm doing wrong. When I fit
with 2 parameters or more, everything works fine.

The calling stack is:

#0  least_gsl_error_handler (
    reason=0x4019be80 "vector length n must be positive integer",
    file=0x4019be4e "subvector_source.c", line=28, gsl_errno=4)
    at least_test.c:11
#1  0x4008fc92 in gsl_error () from /usr/lib/libgsl.so.0
#2  0x40169857 in gsl_vector_subvector () from /usr/lib/libgsl.so.0
#3  0x400ba859 in gsl_linalg_SV_decomp () from /usr/lib/libgsl.so.0
#4  0x400bb4c4 in gsl_linalg_SV_decomp_mod () from /usr/lib/libgsl.so.0
#5  0x400ea171 in gsl_multifit_linear () from /usr/lib/libgsl.so.0
#6  0x080489a5 in main () at least_test.c:50

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>
void least_gsl_error_handler( const char * reason,
			      const char * file,
			      int line,
			      int gsl_errno)
{
  printf("%s:%d: %s (gsl_error: %s)\n", 
	 file, line, reason, gsl_strerror(gsl_errno));
}
int main()
{
  gsl_vector *x,*y,*p;
  gsl_multifit_linear_workspace *workspace;
  gsl_matrix *X, *cov;
  double chi2;
  int i;
  double xi[5] = {1.0, 2.0, 3.0, 4.0, 5.0};
  double yi[5] = {4.0, 8.0, 12.0, 16.0, 20.0};
  gsl_set_error_handler( &least_gsl_error_handler );
  workspace = gsl_multifit_linear_alloc(5, 1);
  cov = gsl_matrix_alloc(1, 1);
  p = gsl_vector_alloc(1);
  x = gsl_vector_alloc(5);
  y = gsl_vector_alloc(5);
  X = gsl_matrix_alloc( 5, 1 );
  for (i = 0; i < 5; i++)
    {
      gsl_vector_set(x, i, xi[i]);
      gsl_vector_set(y, i, yi[i]);
      gsl_matrix_set(X, i, 0, xi[i]); /* straight line */
    }
  gsl_multifit_linear(X, y, p, cov, &chi2, workspace);
  printf(" a = %e +/- %e\n chi2 = %e\n", 
	 gsl_vector_get(p, 0), sqrt(gsl_matrix_get(cov,0,0)), chi2);
  return 0;
}

----------------------------------------------------------------------
