gsl_sf_hyperg_U_int_e10_e function

["Loredo, Elvira" <loredo at rand dot org>]

Below is a C program that uses the gsl_sf_hyperg_U_int_e10_e function.  The
program takes the user supplied parameter values, (n,k,z and r) and
generates
the correct solution over several ranges of these parameters when compared
to the results generated by Mathematica.  However, the C call begins to
break down when the combination of parameters exceeds a certain threshold,
while Mathematica continues to work.  

For example, I use Mathematica to evaluate the same function and compare the
results below:


					Mathematica         C
n = 500, k =8, z = 1800, r = 200               70.844	              70.844
n = 500, k =8, z = 1800, r = 250               56.7043             56.7043
n = 500, k= 8, z = 1800, r = 270	              52.117
52.117

Note that the function is evaluated  at k, 1+k+n and x where x = z*r/k  and
at k, k+n, x=z*r/k, e.g.. gsl_sf_hyperg_U_int_e10_e( k, 1+k+n, x, &thu2 )
When, x= z*r/k <52 the C function call breaks down, while the Mathematica
continues to provide a seemingly correct solution.  
For example
					Mathematica         C

n = 500, k= 8, z = 1800, r = 300              47.269              Illegal
Instruction
n = 500, k =6, z = 1800, r = 200              53.402              53.402
n = 500, k= 5, z = 1800, r = 000              47.269              Illegal
Instruction
n = 500, k =8, z = 1300, r = 200	             51.2012
51.2012
n = 500, k =8, z = 1299, r = 200	             51.2012
Illegal Instruction

The value of n also impacts the C programs ability to evaluate the function:

n = 500, k =8, z = 1299, r = 200	             51.2012
Illegal Instruction
n = 499, k =8, z = 1299, r = 200	             51.1601
51.1601

I am running Redhat Linux 7.2 on a Dell Optiplex GX-150 desktop with an i686
processor.

Any ideas as to why this is happening and any suggested solutions would be
appreciated.

Thanks,

Elvira Loredo



====================================CODE====================================
==========

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_math.h>

double thput( int, int, int, int );

int main( int argc, char *argv[] )
{
  int n,k,z,r; 
  
  n = atoi( argv[1] );
  k = atoi( argv[2] );
  z = atoi( argv[3] );
  r = atoi( argv[4] );

  printf( "%g\n", thput( n, k, z, r ) );
  return 0;
}

double thput( int n, int k, int z, int r )
{
   int status;
   double x, thu;
   gsl_sf_result_e10 thu1, thu2;

   x = (double) ( k * z ) / (double) r;
   status = gsl_sf_hyperg_U_int_e10_e( k, k+n, x, &thu1 );
   fprintf(stderr, "Status of Call 1 = %s\n", gsl_strerror(status));
   status = gsl_sf_hyperg_U_int_e10_e( k, 1+k+n, x, &thu2 );
   fprintf(stderr, "Status of Call 2 = %s\n", gsl_strerror(status));
   thu = (double) n * thu1.val / thu2.val * gsl_pow_int( 10.0, thu1.e10 -
         thu2.e10 );
   return thu;
}