bugs

math library from gnu

	strcpy(filename,"lpdf-");
	strcat(filename,argv[1]);
	f2=fopen(filename,"w+");
	l=(double*)calloc(rand_numbers,sizeof(double));
	lpdf=(double*)calloc(hist_size,sizeof(double));
	r=gsl_rng_alloc (gsl_rng_mt19937);
	gettimeofday(tv,tz);
	rnd_seed=(unsigned long int)tv->tv_usec;
	gsl_rng_set(r,rnd_seed);
	i=0;
	do
	{
		l[i]=0.0;
		for (j=1;j<n;j++) l[i]+=gsl_ran_levy_skew(r,1.0,alpha,0.0);
		l[i]=l[i]/pow(n,1.0/alpha);
		if (abs(l[i])<=20) i++;//picks only random numbers in the interval  
[a,b] to get good precision in the histogram
	}while(i<rand_numbers);
	gsl_histogram_set_ranges_uniform(h,a,b);
	for(i=0;i<rand_numbers;i++) gsl_histogram_increment(h,l[i]);
	gsl_histogram_scale(h,(double)hist_size/((b-a)*gsl_histogram_sum(h)));
	gsl_histogram_fprintf(f1,h,"%g","%g");
	lpdf=levy_pdf(alpha,hist_size,a,b);
	for (i=0;i<hist_size;i++) fprintf(f2,"%e\t%e\n",i*dx+a,lpdf[i]);
	for (i=0;i<hist_size;i++) errabs+=pow((gsl_histogram_get(h,i)-lpdf[i]),2.0);
	printf("%e\n",errabs/(double)hist_size);
	gsl_histogram_free(h);
	gsl_rng_free(r);
	fclose(f1);
	exit (0);
}

> At Mon, 26 Mar 2007 19:48:01 -0300,
> [e-mail] wrote:
>>
>> The Lev� skew random number generator (gsl_ran_levy_skew) does not
>> procuce a Lev� random number when beta=0 (symmetric case), and the
>> gsl_ran_levy function does not work as stated in the docs. I made some
>> histograms from 10^6 samples to check the accuracy of the algorithms,
>> by comparison agaisnt the numerical integration of the equation of
>> Lev�'s PDF. For the gsl_ran_levy function there is a good precison for
>> alpha [1,2], for alpha (0.3,1) you must sum a series of random numbers
>> to get the same precision (tipicaly 100 or more gsl_ran_levy numbers).
>> For alpha<=0.3 the algorithm does not work properly, even worse, the
>> error increases as you add more random numbers. This contradicts the
>> manual that says "the algoritm only works for alpha (0,2]". The
>> function gsl_ran_levy_skew does not produce levy random numbers when
>> beta=0, instead the pdf of the random numbers is a linear (?!?!) one.
>
> Thanks for your email.  Please can you send a small example program
> which demonstrates the problem.
>
> Note that the Levy skew generator is tested in the GSL test suite for
> several cases where beta=0 -- if you have not done so, can you run
> "make check" and confirm that it works for these cases.
>
> --
> Brian Gough
> (GSL Maintainer)
>
> Network Theory Ltd,
> Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/
>

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

_______________________________________________
Bug-gsl mailing list
[e-mail]
http://lists.gnu.org/mailman/listinfo/bug-gsl

-

------------------------------------------------------------------------
BUG-ID:   21832   
STATUS:   Open/Fixed      
CATEGORY: None
SUMMARY:  gsl infnan.c configure/build bug on solaris

From: Richard Smith <[e-mail]>
To: [e-mail]
Subject: [Bug-gsl] gsl infnan.c configure/build bug
Date: Tue, 24 Jul 2007 10:24:22 +1000

infnan.c currently fails to compile on Solaris with Sun Studio 12 compilers
when using default compiler options as generated by running configure with
no options. Here is an extract of the output when running make:

/bin/bash ../libtool --tag=CC --mode=compile cc -DHAVE_CONFIG_H  -I. 
-I../../sys
 -I.. -I..     -g -c -o infnan.lo ../../sys/infnan.c
 cc -DHAVE_CONFIG_H -I. -I../../sys -I.. -I.. -g -c ../../sys/infnan.c  
-KPIC -D
PIC -o .libs/infnan.o
"/usr/include/ieeefp.h", line 74: syntax error before or at: 
__builtin_isfinite
cc: acomp failed for ../../sys/infnan.c

The problem seems to occur as a result of the following combination of 
things:
1. "finite" function is used in various source files
2. configure.ac checks for <ieeefp.h> and finds it
3. On Solaris, "finite" is defined in <ieeefp.h>
4. configure.ac only checks <math.h> for "finite", and doesn't find it.
5. Since its a c99 environment, "isfinite" is found.
6. config.h ends up containing amongst other things the following lines:
#define HAVE_DECL_FINITE 0
#define HAVE_DECL_ISFINITE 1
#if !HAVE_DECL_FINITE
#if HAVE_DECL_ISFINITE
#define finite isfinite
#else
#define finite gsl_finite
#endif
#endif
   This means that "finite" macro has the value "isfinite"
7. "isfinite" is a macro, ultimately defined in <iso/math_c99.h>:
#define        isfinite(x)     __builtin_isfinite(x)
8. <ieeefp.h> declares "finite" function:
extern int      finite(double);
9. After macro expansion the compiler sees
extern int      __builtin_isfinite ( double ) ;
and complains with an error.

There's probably multiple ways of overcoming the problem. A workaround is
to force the compiler not to use a C99 environment e.g. -xc99=%none. However
since it would be desirable to have configure work well by default, maybe
the test in configure.ac for "finite" should be something like:
AC_CHECK_DECLS(finite,,,[#include <math.h>
#if HAVE_IEEEFP_H
# include <ieeefp.h>
#endif])

Alternatively all uses of "finite" could be changed to use "isfinite", and
the configuration/build be based around assuming a C99 environment, with
substitute functions and macros being provided where the necessary C99
facilities are missing. On Linux, the BSD floating point classification
functions are documented as being obsolete, so their use should probably
be avoided.

-- 
============================================================================
   ,-_|\   Richard Smith - Technical Specialist
  /     \  Sun Microsystems Australia         Phone : [phone]
[e-mail]                        Direct : [phone]
  \_,-._/  476 St Kilda Road                    Fax : [phone]
       v   Melbourne Vic 3004 Australia
=========================================================================== 

_______________________________________________
Bug-gsl mailing list
[e-mail]
http://lists.gnu.org/mailman/listinfo/bug-gsl

------------------------------------------------------------------------
BUG-ID:   21833   
STATUS:   Open            
CATEGORY: Performance
SUMMARY:  suboptimal performance of gsl permutation?

From: "djamel anonymous" <[e-mail]>
To: [e-mail]
Subject: gsl permutation
Date: Wed, 24 Jan 2007 07:42:06 +0000

Hi.
i am sending you this email about a possible issue in glibc permutation 
algorithm.i think it has qudratic worst case running time.for example if we 
have the permutation
1,2,3,4,,,,,n,0
we will permute all elements in first iteration then for each iteration we 
will traverse on average half the cycle (which is of size n) before we know 
that the elements of cycle have already been permuted.there is two possible 
solutions to make the algorithm linear:
1-at each step we traverse the full cycle and permute all elements in the 
cycle.for each permuted element i  we assign p[i]=i.the disavantage of this 
is that it will destroy the original permutation.
2-use a bit array of size n bits.each time we permute an element we set the 
relevant bit.if at iteration i we find that bit i is already set we skip to 
next iteration.the disavantage of this is that it equires additional storage 
allocation.
best regards.

_________________________________________________________________
MSN Hotmail sur i-mode™ : envoyez et recevez des e-mails depuis votre 
téléphone portable ! http://www.msn.fr/hotmailimode/

3 - Normal

------------------------------------------------------------------------
BUG-ID:   21835   
STATUS:   Open            
CATEGORY: Accuracy problem
SUMMARY:  gsl_sf_hyperg_2F1 problematic arguments

BUG#1 -- gsl_sf_hyperg_2F1_e fails for some arguments 

From: [e-mail]
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

As reported by Lee Warren <[e-mail]> another set of
arguments which fail are:

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

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

The correct value is -2.

See also, 

From: Olaf Wucknitz <[e-mail]>
To: [e-mail]
Subject: [Bug-gsl] gsl_sf_hyperg_2F1

Hi,

I am having a problem with gsl_sf_hyperg_2F1.
With the parameters (-0.5, 0.5, 1, x) the returned values show a jump at 
x=0.5. For x<0.5 the results seem to be correct, while for x>0.5 they 
aren't.
The function gsl_sf_hyperg_2F1_e calls hyperg_2F1_series for x<0.5, but
hyperg_2F1_reflect for x>0.5. The latter function checks for c-a-b being 
an integer (which it is in my case). If I change one of the parameters 
a,b,c by a small amount, the other branch of the function is taken and the 
results are correct again.
Unfortunately I know too little about the numerics of 2F1 to suggest a 
patch at the moment.

Regards,
Olaf Wucknitz
-- 
Joint Institute for VLBI in Europe                [e-mail]

------------------------------------------------------------------------
BUG-ID:   21836   
STATUS:   Open/Confirmed  
CATEGORY: Accuracy problem
SUMMARY:  gamma_inc_P and gamma_inc_Q only satisfy P+Q=1 within errors

BUG#44 -- gamma_inc_P and gamma_inc_Q only satisfy P+Q=1 within errors

The sum of gamma_inc_P and gamma_inc_Q doesn't always satisfy the
identity P+Q=1 exactly (although it is correct within errors), due the
slightly different branch conditions for the series and continued
fraction expansions.  These could be made identical so that P+Q=1 exactly.

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

int
main (void)
{
  gsl_sf_result r1, r2;
  double a = 0.3, x = 1.0;
  gsl_sf_gamma_inc_P_e (a, x, &r1);
  gsl_sf_gamma_inc_Q_e (a, x, &r2);
  printf("%.18e\n", r1.val);
  printf("%.18e\n", r2.val);
  printf("%.18e\n", r1.val + r2.val);
}

$ ./a.out
9.156741562411074842e-01
8.432584375889111417e-02
9.999999999999985567e-01

3 - NormalNone

------------------------------------------------------------------------
BUG-ID:   21837   
STATUS:   Open/Confirmed  
CATEGORY: Runtime error
SUMMARY:  gsl_linalg_solve_symm_tridiag requires positive definite matrix

A zero on the diagonal will cause NaNs even though a reasonable
solution could be computed in principle. 

#include <gsl/gsl_linalg.h>

int main (void)
{
  double d[] = { 0.00, 1.21, 0.80, 1.55, 0.76 } ;
  double e[] = { 0.82, 0.39, 0.09, 0.68 } ;
  double b[] = { 0.07, 0.62, 0.81, 0.11, 0.65} ;
  double x[] = { 0.00, 0.00, 0.00, 0.00, 0.00} ;

  gsl_vector_view dv = gsl_vector_view_array(d, 5);
  gsl_vector_view ev = gsl_vector_view_array(e, 4);
  gsl_vector_view bv = gsl_vector_view_array(b, 5);
  gsl_vector_view xv = gsl_vector_view_array(x, 5);

  gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector);
  gsl_vector_fprintf(stdout, &xv.vector, "% .5f");

  d[0] += 1e-5;
  gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector);
  gsl_vector_fprintf(stdout, &xv.vector, "% .5f");
}

$ ./a.out
 nan
 nan
 nan
 nan
 nan
 0.13626
 0.08536
 1.03840
-0.60009
 1.39219

AUG 2007: We now return an error code for this case.  To return a solution
we would need to do a permutation, see slatec/dgtsl.f

3 - NormalNone

------------------------------------------------------------------------
BUG-ID:   22478   
STATUS:   Open/Confirmed  
CATEGORY: None
SUMMARY:  Missing functions for complex vectors

From: Federico Zenith <[e-mail]>
To: [e-mail]
Subject: [Help-gsl] Missing functions for complex vectors
Date: Mon, 03 Mar 2008 13:50:43 +0100

Hi,
I am working on a C++ GSL wrapper (yes, there are not enough of them
already around :-), and I have been working on vectors last week; while
writing the wrapper, I noticed a couple of oddities I'd like to point
out before reporting them as a bug (I am using the GSL version 1.10):

1) complex vectors have defined all properties mentioned on this page:
http://www.gnu.org/software/gsl/manual/html_node/Vector-properties.html
except isnonneg(). I am not sure what's the meaning of ispos() or
isneg() for a vector of complex numbers, I suppose it means "both real
AND imaginary part are positive", but this should be documented;
furthermore, I do not understand why, if ispos() and isneg() could be
defined, why could not isnonneg() be as well.

2) Vector operations, defined on this page:
http://www.gnu.org/software/gsl/manual/html_node/Vector-operations.html
are not defined for complex types. As far as I can guess, all these
operations would still make sense for complex numbers. What is really
odd is, the same operations are defined for complex matrices.

Did I miss something?

Cheers,
-Federico

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