mas01mc@742: /*************************************
mas01mc@742:  An ANSI-C implementation of the digamma-function for real arguments based
mas01mc@742:  on the Chebyshev expansion proposed in appendix E of 
mas01mc@742:  http://arXiv.org/abs/math.CA/0403344 . For other implementations see
mas01mc@742:  the GSL implementation for Psi(Digamma) in
mas01mc@742:  http://www.gnu.org/software/gsl/manual/gsl-ref_toc.html
mas01mc@742: 
mas01mc@742: Richard J. Mathar, 2005-11-24
mas01mc@742: http://www.strw.leidenuniv.nl/~mathar/progs/digamma.c
mas01mc@742: **************************************/
mas01mc@742: #include <math.h>
mas01mc@742: 
mas01mc@742: #ifndef M_PIl
mas01mc@742: /** The constant Pi in high precision */
mas01mc@742: #define M_PIl 3.1415926535897932384626433832795029L
mas01mc@742: #endif
mas01mc@742: #ifndef M_GAMMAl
mas01mc@742: /** Euler's constant in high precision */
mas01mc@742: #define M_GAMMAl 0.5772156649015328606065120900824024L
mas01mc@742: #endif
mas01mc@742: #ifndef M_LN2l
mas01mc@742: /** the natural logarithm of 2 in high precision */
mas01mc@742: #define M_LN2l 0.6931471805599453094172321214581766L
mas01mc@742: #endif
mas01mc@742: 
mas01mc@742: /** The digamma function in long double precision.
mas01mc@742: * @param x the real value of the argument
mas01mc@742: * @return the value of the digamma (psi) function at that point
mas01mc@742: * @author Richard J. Mathar
mas01mc@742: * @since 2005-11-24
mas01mc@742: */
mas01mc@742: long double digammal(long double x)
mas01mc@742: {
mas01mc@742: 	/* force into the interval 1..3 */
mas01mc@742: 	if( x < 0.0L )
mas01mc@742: 		return digammal(1.0L-x)+M_PIl/tanl(M_PIl*(1.0L-x)) ;	/* reflection formula */
mas01mc@742: 	else if( x < 1.0L )
mas01mc@742: 		return digammal(1.0L+x)-1.0L/x ;
mas01mc@742: 	else if ( x == 1.0L)
mas01mc@742: 		return -M_GAMMAl ;
mas01mc@742: 	else if ( x == 2.0L)
mas01mc@742: 		return 1.0L-M_GAMMAl ;
mas01mc@742: 	else if ( x == 3.0L)
mas01mc@742: 		return 1.5L-M_GAMMAl ;
mas01mc@742: 	else if ( x > 3.0L)
mas01mc@742: 		/* duplication formula */
mas01mc@742: 		return 0.5L*(digammal(x/2.0L)+digammal((x+1.0L)/2.0L))+M_LN2l ;
mas01mc@742: 	else
mas01mc@742: 	{
mas01mc@742: 		/* Just for your information, the following lines contain
mas01mc@742: 		* the Maple source code to re-generate the table that is
mas01mc@742: 		* eventually becoming the Kncoe[] array below
mas01mc@742: 		* interface(prettyprint=0) :
mas01mc@742: 		* Digits := 63 :
mas01mc@742: 		* r := 0 :
mas01mc@742: 		* 
mas01mc@742: 		* for l from 1 to 60 do
mas01mc@742: 		* 	d := binomial(-1/2,l) :
mas01mc@742: 		* 	r := r+d*(-1)^l*(Zeta(2*l+1) -1) ;
mas01mc@742: 		* 	evalf(r) ;
mas01mc@742: 		* 	print(%,evalf(1+Psi(1)-r)) ;
mas01mc@742: 		*o d :
mas01mc@742: 		* 
mas01mc@742: 		* for N from 1 to 28 do
mas01mc@742: 		* 	r := 0 :
mas01mc@742: 		* 	n := N-1 :
mas01mc@742: 		*
mas01mc@742:  		*	for l from iquo(n+3,2) to 70 do
mas01mc@742: 		*		d := 0 :
mas01mc@742:  		*		for s from 0 to n+1 do
mas01mc@742:  		*		 d := d+(-1)^s*binomial(n+1,s)*binomial((s-1)/2,l) :
mas01mc@742:  		*		od :
mas01mc@742:  		*		if 2*l-n > 1 then
mas01mc@742:  		*		r := r+d*(-1)^l*(Zeta(2*l-n) -1) :
mas01mc@742:  		*		fi :
mas01mc@742:  		*	od :
mas01mc@742:  		*	print(evalf((-1)^n*2*r)) ;
mas01mc@742:  		*od :
mas01mc@742:  		*quit :
mas01mc@742: 		*/
mas01mc@742: 		static long double Kncoe[] = { .30459198558715155634315638246624251L,
mas01mc@742: 		.72037977439182833573548891941219706L, -.12454959243861367729528855995001087L,
mas01mc@742: 		.27769457331927827002810119567456810e-1L, -.67762371439822456447373550186163070e-2L,
mas01mc@742: 		.17238755142247705209823876688592170e-2L, -.44817699064252933515310345718960928e-3L,
mas01mc@742: 		.11793660000155572716272710617753373e-3L, -.31253894280980134452125172274246963e-4L,
mas01mc@742: 		.83173997012173283398932708991137488e-5L, -.22191427643780045431149221890172210e-5L,
mas01mc@742: 		.59302266729329346291029599913617915e-6L, -.15863051191470655433559920279603632e-6L,
mas01mc@742: 		.42459203983193603241777510648681429e-7L, -.11369129616951114238848106591780146e-7L,
mas01mc@742: 		.304502217295931698401459168423403510e-8L, -.81568455080753152802915013641723686e-9L,
mas01mc@742: 		.21852324749975455125936715817306383e-9L, -.58546491441689515680751900276454407e-10L,
mas01mc@742: 		.15686348450871204869813586459513648e-10L, -.42029496273143231373796179302482033e-11L,
mas01mc@742: 		.11261435719264907097227520956710754e-11L, -.30174353636860279765375177200637590e-12L,
mas01mc@742: 		.80850955256389526647406571868193768e-13L, -.21663779809421233144009565199997351e-13L,
mas01mc@742: 		.58047634271339391495076374966835526e-14L, -.15553767189204733561108869588173845e-14L,
mas01mc@742: 		.41676108598040807753707828039353330e-15L, -.11167065064221317094734023242188463e-15L } ;
mas01mc@742: 
mas01mc@742: 		register long double Tn_1 = 1.0L ;	/* T_{n-1}(x), started at n=1 */
mas01mc@742: 		register long double Tn = x-2.0L ;	/* T_{n}(x) , started at n=1 */
mas01mc@742: 		register long double resul = Kncoe[0] + Kncoe[1]*Tn ;
mas01mc@742: 		register int n;
mas01mc@742: 
mas01mc@742: 		x -= 2.0L ;
mas01mc@742: 
mas01mc@742: 		for(n = 2 ; n < sizeof(Kncoe)/sizeof(long double) ;n++)
mas01mc@742: 		{
mas01mc@742: 			const long double Tn1 = 2.0L * x * Tn - Tn_1 ;	/* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */
mas01mc@742: 			resul += Kncoe[n]*Tn1 ;
mas01mc@742: 			Tn_1 = Tn ;
mas01mc@742: 			Tn = Tn1 ;
mas01mc@742: 		}
mas01mc@742: 		return resul ;
mas01mc@742: 	}
mas01mc@742: }
mas01mc@742: 
mas01mc@742: /** The optional interface to CREASO's IDL is added if someone has defined
mas01mc@742: * the cpp macro export_IDL_REF, which typically has been done by including the
mas01mc@742: * files stdio.h and idl_export.h before this one here.
mas01mc@742: */
mas01mc@742: #ifdef export_IDL_REF
mas01mc@742: /** CALL_EXTERNAL interface.
mas01mc@742: * dg = CALL_EXTERNAL('digamma.so',X)
mas01mc@742: * @pararm argc the number of arguments. This is supposed to be 1 and not 
mas01mc@742: *    checked further because it might have negative impact on performance.
mas01mc@742: * @pararm argv the parameter list. The first element is the parameter x
mas01mc@742: *    supposed to be of type DOUBLE in IDL 
mas01mc@742: * @return the return value, again of IDL-type DOUBLE
mas01mc@742: * @since 2007-01-16
mas01mc@742: * @author Richard J. Mathar
mas01mc@742: */
mas01mc@742: double digammal_idl(int argc, void *argv[])
mas01mc@742: {
mas01mc@742: 	long double x = *(double*)argv[0] ;
mas01mc@742: 	return (double)digammal(x) ;
mas01mc@742: }
mas01mc@742: #endif /* export_IDL_REF */
mas01mc@742: 
mas01mc@742: 
mas01mc@742: #ifdef MAIN
mas01mc@742: 
mas01mc@742: #include <stdlib.h>
mas01mc@742: #include <stdio.h>
mas01mc@742: 
mas01mc@742: int main(int argc, char **argv){
mas01mc@742: 	if (argc < 2){
mas01mc@742: 		fprintf(stderr, "Syntax: %s x\n", argv[0]);
mas01mc@742: 		exit(1);
mas01mc@742: 	}
mas01mc@742: 	double x = atof(argv[1]);
mas01mc@742: 	double y = digammal(x);
mas01mc@742: 	printf("%g\n", y);
mas01mc@742: 	return 0;
mas01mc@742: }
mas01mc@742: #endif	/* MAIN */
mas01mc@742: 
mas01mc@742: #ifdef TEST
mas01mc@742: 
mas01mc@742: /* an alternate implementation for test purposes, using formula 6.3.16 of Abramowitz/Stegun with the 
mas01mc@742:    first n terms */
mas01mc@742: #include <stdio.h>
mas01mc@742: 
mas01mc@742: long double digammalAlt(long double x, int n)
mas01mc@742: {
mas01mc@742: 	/* force into the interval 1..3 */
mas01mc@742: 	if( x < 0.0L )
mas01mc@742: 		return digammalAlt(1.0L-x,n)+M_PIl/tanl(M_PIl*(1.0L-x)) ;	/* reflection formula */
mas01mc@742: 	else if( x < 1.0L )
mas01mc@742: 		return digammalAlt(1.0L+x,n)-1.0L/x ;
mas01mc@742: 	else if ( x == 1.0L)
mas01mc@742: 		return -M_GAMMAl ;
mas01mc@742: 	else if ( x == 2.0L)
mas01mc@742: 		return 1.0L-M_GAMMAl ;
mas01mc@742: 	else if ( x == 3.0L)
mas01mc@742: 		return 1.5L-M_GAMMAl ;
mas01mc@742: 	else if ( x > 3.0L)
mas01mc@742: 		return digammalAlt(x-1.0L,n)+1.0L/(x-1.0L) ;
mas01mc@742: 	else
mas01mc@742: 	{
mas01mc@742: 		x -= 1.0L ;
mas01mc@742: 		register long double resul = -M_GAMMAl ;
mas01mc@742: 
mas01mc@742: 		for( ; n >= 1 ;n--)
mas01mc@742: 			resul += x/(n*(n+x)) ;
mas01mc@742: 		return resul ;
mas01mc@742: 	}
mas01mc@742: }
mas01mc@742: int main(int argc, char *argv[])
mas01mc@742: {
mas01mc@742: 	long double x;
mas01mc@742: 	for(x=0.01 ; x < 5. ; x += 0.02)
mas01mc@742: 		printf("%.2Lf %.30Lf %.30Lf %.30Lf\n",x, digammal(x), digammalAlt(x,100), digammalAlt(x,200) ) ;
mas01mc@742: }
mas01mc@742: 
mas01mc@742: #endif /* TEST */