problem with gsl_sf_lnpoch

problem with gsl_sf_lnpoch

[Giulio Bottazzi]

[-- Attachment #1: Type: text/plain, Size: 769 bytes --]

Hi,
It seems that the exponential of the logarithm of the Pochhammer symbol
with second argument equal to zero is NOT one. Consider the following
little program, which should print two times the number 1 separated by
a space

/* -- START test.c */
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
int main(){

  printf("exp(log(poch))=%g poch=%g\n",
        exp( gsl_sf_lnpoch(7,0)),gsl_sf_poch(7,0));

  return 0;
}
/* -- END test.c */

If I compile it with
#gcc -Wall test.c -lgsl -lgslcblas -lm -o test
and run
#./test
I get:
#exp(log(poch))=2.71828 poch=1

What's going on? I can't see any mistake in my code.

Best,
	Giulio.


--
Giulio Bottazzi                       PGP Key ID:BAB0A33F
giulio.bottazzi@libero.it   http://www.sssup.it/~bottazzi

[-- Attachment #2: Type: application/pgp-signature, Size: 189 bytes --]

Re: problem with gsl_sf_lnpoch

[James Theiler]

I think Giulio is right; I took a quick look at the code,
and the value when the second argument is zero is hardcoded;
ie:
 
int
gsl_sf_lnpoch_e(const double a, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(a <= 0.0 || a+x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 1.0;   <=================== should be 0.0
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    return lnpoch_pos(a, x, result);
  }
}

jt


On Fri, 29 Sep 2006, Giulio Bottazzi wrote:

] Hi,
] It seems that the exponential of the logarithm of the Pochhammer symbol
] with second argument equal to zero is NOT one. Consider the following
] little program, which should print two times the number 1 separated by
] a space
] 
] /* -- START test.c */
] #include <stdio.h>
] #include <math.h>
] #include <gsl/gsl_sf_gamma.h>
] int main(){
] 
]   printf("exp(log(poch))=%g poch=%g\n",
]         exp( gsl_sf_lnpoch(7,0)),gsl_sf_poch(7,0));
] 
]   return 0;
] }
] /* -- END test.c */
] 
] If I compile it with
] #gcc -Wall test.c -lgsl -lgslcblas -lm -o test
] and run
] #./test
] I get:
] #exp(log(poch))=2.71828 poch=1
] 
] What's going on? I can't see any mistake in my code.
] 
] Best,
] 	Giulio.
] 
] 
] --
] Giulio Bottazzi                       PGP Key ID:BAB0A33F
] giulio.bottazzi@libero.it   http://www.sssup.it/~bottazzi
] 

-- 
James Theiler
MS-B244, ISR-2, LANL; Los Alamos, NM 87544
Space and Remote Sensing Sciences; Los Alamos National Laboratory
http://public.lanl.gov/jt

Re: problem with gsl_sf_lnpoch

[Brian Gough]

Giulio Bottazzi wrote:
> Hi,
> It seems that the exponential of the logarithm of the Pochhammer symbol
> with second argument equal to zero is NOT one. Consider the following
> little program, which should print two times the number 1 separated by
> a space

Thanks for the bug report. I will fix that.

-- 
Brian Gough

Network Theory Ltd,
Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/