Re: problem with gsl_sf_lnpoch

[James Theiler <jt at lanl dot gov>]

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