1 #include "mconf.h"
2 
3 #define EPS 1.0e-17
4 
besselpoly(double a,double lambda,double nu)5 double besselpoly(double a, double lambda, double nu) {
6 
7   int m, factor=0;
8   double Sm, relerr, Sol;
9   double sum=0.0;
10 
11   /* Special handling for a = 0.0 */
12   if (a == 0.0) {
13     if (nu == 0.0) return 1.0/(lambda + 1);
14     else return 0.0;
15   }
16   /* Special handling for negative and integer nu */
17   if ((nu < 0) && (floor(nu)==nu)) {
18     nu = -nu;
19     factor = ((int) nu) % 2;
20   }
21   Sm = exp(nu*log(a))/(Gamma(nu+1)*(lambda+nu+1));
22   m = 0;
23   do {
24     sum += Sm;
25     Sol = Sm;
26     Sm *= -a*a*(lambda+nu+1+2*m)/((nu+m+1)*(m+1)*(lambda+nu+1+2*m+2));
27     m++;
28     relerr = fabs((Sm-Sol)/Sm);
29   } while (relerr > EPS && m < 1000);
30   if (!factor)
31     return sum;
32   else
33     return -sum;
34 }
35