Index: cephes/cephes.h
===================================================================
RCS file: /cvsroot/gretl/gretl/cephes/cephes.h,v
retrieving revision 1.7
diff -r1.7 cephes.h
37,42c37,42
< double j0 (double x); /* Bessel J_0(x) */
< double j1 (double x); /* Bessel J_1(x) */
< double y0 (double x); /* Bessel Y_0(x) */
< double y1 (double x); /* Bessel Y_1(x) */
< double i0 (double x); /* Bessel I_0(x) */
< double i1 (double x); /* Bessel I_1(x) */
---
> double hyperg( double a, double b, double x);
> double cephes_bessel_Kn (int nn, double x);
> double cephes_bessel_I0 (double x);
> double cephes_bessel_I1 (double x);
> double cephes_bessel_K0 (double x);
> double cephes_bessel_K1 (double x);
Index: cephes/i0.c
===================================================================
RCS file: /cvsroot/gretl/gretl/cephes/i0.c,v
retrieving revision 1.2
diff -r1.2 i0.c
153c153
< double i0 (double x)
---
> double cephes_bessel_I0 (double x)
Index: cephes/i1.c
===================================================================
RCS file: /cvsroot/gretl/gretl/cephes/i1.c,v
retrieving revision 1.2
diff -r1.2 i1.c
153c153
< double i1 (double x)
---
> double cephes_bessel_I1(double x)
Index: cephes/k0.c
===================================================================
RCS file: /cvsroot/gretl/gretl/cephes/k0.c,v
retrieving revision 1.2
diff -r1.2 k0.c
134c134
< double k0 (double x)
---
> double cephes_bessel_K0 (double x)
145c145
< 	y = chbevl(y, A, 10) - log(0.5 * x) * i0(x);
---
> 	y = chbevl(y, A, 10) - log(0.5 * x) * cephes_bessel_I0(x);
166c166
< 	y = chbevl(y, A, 10) - log(0.5 * x) * i0(x);
---
> 	y = chbevl(y, A, 10) - log(0.5 * x) * cephes_bessel_I0(x);
Index: cephes/k1.c
===================================================================
RCS file: /cvsroot/gretl/gretl/cephes/k1.c,v
retrieving revision 1.2
diff -r1.2 k1.c
135c135
< double k1 (double x)
---
> double cephes_bessel_K1 (double x)
147c147
< 	y =  log(z) * i1(x) + chbevl(y, A, 11) / x;
---
> 	y =  log(z) * cephes_bessel_I1(x) + chbevl(y, A, 11) / x;
165c165
< 	y =  log(0.5 * x) * i1(x) + chbevl(y, A, 11) / x;
---
> 	y =  log(0.5 * x) * cephes_bessel_I1(x) + chbevl(y, A, 11) / x;
Index: cephes/kn.c
===================================================================
RCS file: /cvsroot/gretl/gretl/cephes/kn.c,v
retrieving revision 1.2
diff -r1.2 kn.c
84c84
< double kn (int nn, double x)
---
> double cephes_bessel_Kn (int nn, double x)
Index: lib/src/geneval.c
===================================================================
RCS file: /cvsroot/gretl/gretl/lib/src/geneval.c,v
retrieving revision 1.418
diff -r1.418 geneval.c
31a32
> #include "../../cephes/cephes.h"
2482a2484,2492
> static int gretl_round(double x)
> {
> 	if (x < 0) {
> 	    return (x - floor(x) <= 0.5)? floor(x) : ceil(x);
> 	} else {
> 	    return (x - floor(x) < 0.5)? floor(x) : ceil(x);
> 	}
> }
> 
2517,2521c2527
< 	if (x < 0) {
< 	    return (x - floor(x) <= 0.5)? floor(x) : ceil(x);
< 	} else {
< 	    return (x - floor(x) < 0.5)? floor(x) : ceil(x);
< 	}
---
> 	return gretl_round(x);
4705a4712,4896
> /* evaluates sclar bessel function */
> static int eval_bessel_scalar_func(char *type, double lambda, double x, double *res)
> {
> 	if (*type=='J') {
> 		if ((x<0) && (!(gretl_round(lambda)==lambda))) {
> 			/* negative x and non-integer lambda is not supported */
> 			*res=NADBL;
> 			return -2;
> 		}
> 		*res=cephes_bessel_Jv(lambda,x);
> 		return 0;
> 	} else if (*type=='I') {
> 		if ((x<0) && (!(gretl_round(lambda)==lambda))) {
> 			/* negative x and non-integer lambda is not supported */
> 			*res=NADBL;
> 			return -2;
> 		}
> 		if (lambda<0) {
> 			if (fabs(((double)gretl_round(2*lambda))/2 - (lambda)) < 0.00001) {
> 				/*  interpolation around (half) integer number */
> 				double lambda2=((double)gretl_round(2*lambda))/2+0.0001;
> 				double up;
> 				double down;
> 				char typeI='I';
> 				eval_bessel_scalar_func(&typeI, lambda2, x, &up);
> 				lambda2=((double)gretl_round(2*lambda))/2-0.0001;
> 				eval_bessel_scalar_func(&typeI, lambda2, x, &down);
> 				*res=down+(up-down)*((lambda-lambda2)/0.0002);
> 			} else {
> 				double lambda2=floor(fabs(lambda));				
> 				if (
> 					( (lambda2==0) && (x>7.4 ) ) ||
> 					( (lambda2==1) && (x>8   ) ) ||
> 					( (lambda2==2) && (x>8.2 ) ) ||
> 					( (lambda2==3) && ((x>8.3 ) || (x<1.0) ) ) ||
> 					( (lambda2==4) && ((x>8.5 ) || (x<1.6) ) ) ||
> 					( (lambda2==5) && ((x>8.8 ) || (x<2.0) ) ) ||
> 					( (lambda2==6) && ((x>9.9 ) || (x<2.3) ) ) ||	
> 					( (lambda2==7) && ((x>11.6) || (x<3.5) ) ) ||
> 					( (lambda2==8) && ((x>13  ) || (x<4.1) ) ) ||
> 					( (lambda2==9) && ((x>13  ) || (x<5  ) ) ) ||
> 					( (lambda2 >9) && (x>13) )
> 					)
> 				{
> 					*res=exp(-x)*pow(0.5*x,lambda)*hyperg(lambda+0.5,1+2*lambda,2*x)/cephes_gamma(lambda+1);
> 				} else {		
> 					*res=exp(x)*pow(0.5*x,lambda)*hyperg(lambda+0.5,1+2*lambda,-2*x)/cephes_gamma(lambda+1);
> 				}
> 			}
> 		} else if (lambda==0) {
> 			*res=cephes_bessel_I0(x);
> 		} else if (lambda==1) {
> 			*res=cephes_bessel_I1(x);
> 		} else {
> 			*res=cephes_bessel_Iv(lambda,x);
> 		}
> 		return 0;
> 	} else if (*type=='Y') {
> 		if ((x<0) && (!(gretl_round(lambda)==lambda))) {
> 			/* negative x and non-integer lambda is not supported */
> 			*res=NADBL;
> 			return -2;
> 		}
> 		*res=cephes_bessel_Yv(lambda,x);
> 		return 0;
> 	} else if (*type=='K') {
> 		/* besselk is symmetric around lambda==0 */
> 		lambda=fabs(lambda);
> 		if (x<0) {
> 			/* negative x is not supported */
> 			*res=NADBL;
> 			return -2;
> 		} else if (lambda==1) {
> 			*res=cephes_bessel_K1(x);
> 		} else if (lambda==0) {
> 			*res=cephes_bessel_K0(x);
> 		} else if (gretl_round(lambda) == lambda) {
> 			*res=cephes_bessel_Kn(lambda,x);
> 		} else if (fabs(gretl_round(lambda) - lambda) < 0.000001) {
> 			/*  interpolation around integer number */
> 			char typeK='K';
> 			double lambda2=gretl_round(lambda)+0.00001;
> 			double up,down=0;
> 			eval_bessel_scalar_func(&typeK,lambda2,x,&up);
> 			lambda2=gretl_round(lambda)-0.00001;
> 			eval_bessel_scalar_func(&typeK,lambda2,x,&down);
> 			*res=down+(up-down)*((lambda-lambda2)/0.00002);
> 		} else {
> 			char typeI='I';
> 			double besselIp=0;
> 			double besselIm=0;
> 			eval_bessel_scalar_func(&typeI,-1*lambda,x,&besselIm);
> 			eval_bessel_scalar_func(&typeI,lambda,x,&besselIp);
> 			*res=0.5*PI*(besselIm-besselIp)/sin(lambda*PI);			
> 		}
> 		return 0;
> 	}
> 	/* unknown bessel function */
> 	return -1;
> }
> 
> /* generic bessel function handler */
> static NODE *eval_bessel_func (NODE *l, NODE *m, NODE *r, parser *p) 
> {
> 	char *besselFunction = l->v.str;
> 	int status=0;
> 	
> 	NODE *ret=NULL;
> 	double besselRet=0;
> 	if (r->t == NUM) {
> 		double lambda=m->v.xval;
> 		ret = aux_scalar_node(p);
> 	
> 		double x=r->v.xval;
> 
>     	if (ret != NULL) {
> 			status=eval_bessel_scalar_func(besselFunction,lambda,x,&(ret->v.xval));
>     	}
> 	} else if (r->t == MAT) {
> 		double lambda=m->v.xval;
> 		ret = aux_matrix_node(p);
> 
> 		if (ret != NULL && starting(p)) {
> 			const gretl_matrix *x = r->v.m;
> 			int i, n = x->rows * x->cols;
> 
> 			if (node_allocate_matrix(ret, x->rows, x->cols, p)) {
> 				free_tree(ret, p, "On error");
> 				return NULL;
> 			}
> 
> 			for (i=0; i<n && !p->err; i++) {
> 				/* FIXME error handling? */ 
> 				status=eval_bessel_scalar_func(besselFunction,lambda,x->val[i],&besselRet);
> 				ret->v.m->val[i] = besselRet;
> 				if (status==-1) {
> 					break;
> 				}						
> 			}
> 		}
>     } else if (r->t == VEC) {
> 		ret = aux_vec_node(p, p->dinfo->n);
> 		int t, t1, t2;
> 
> 		if (ret != NULL) {
> 			const double *x;
> 
> 			x = r->v.xvec;
> 
> 			if (!p->err) {
> 				t1 = (autoreg(p))? p->obs : p->dinfo->t1;
> 				t2 = (autoreg(p))? p->obs : p->dinfo->t2;
> 				if (scalar_node(m)){
> 					double lambda=m->v.xval;
> 					for (t=t1; t<=t2; t++) {
> 						status=eval_bessel_scalar_func(besselFunction,lambda,x[t],&besselRet);
> 						ret->v.xvec[t] = besselRet;
> 						if (status==-1) {
> 							break;
> 						}						
> 
> 					}
> 				} else {
> 					const double *lambda;
> 					lambda=m->v.xvec;
> 					for (t=t1; t<=t2; t++) {
> 						status=eval_bessel_scalar_func(besselFunction,lambda[t],x[t],&besselRet);
> 						ret->v.xvec[t] = besselRet;
> 						if (status==-1) {
> 							break;
> 						}						
> 
> 					}
> 				}
> 			}
> 		}
> 	}
> 
> 	if (status==-1) {
> 		p->err=E_ARGS;
> 	}
> 	return ret;
> }
> 
> 
6606a6798,6811
>     case F_BESSEL:
> 	/* functions taking one char, one scalar/series and one matrix/series/scalars as args */
> 		if (l->t != STR) {
> 			node_type_error(t->t, 1,
> 				STR, l, p);
> 		} else if (scalar_node(m) && ( scalar_node(r) || (r->t == MAT) || (r->t == VEC))) {
> 			ret = eval_bessel_func(l,m,r,p);
> 		} else if ( (m->t == VEC) && (r->t == VEC)) {
> 			ret = eval_bessel_func(l,m,r,p);		
> 		} else {
> 			node_type_error(t->t, 3,
> 				NUM, r, p);
> 		}		
> 		break;
Index: lib/src/genlex.c
===================================================================
RCS file: /cvsroot/gretl/gretl/lib/src/genlex.c,v
retrieving revision 1.197
diff -r1.197 genlex.c
321a322
>     { F_BESSEL,   "bessel" },
Index: lib/src/genparse.h
===================================================================
RCS file: /cvsroot/gretl/gretl/lib/src/genparse.h,v
retrieving revision 1.177
diff -r1.177 genparse.h
300a301
>     F_BESSEL,
