[genius] Tue Nov 27 12:59:21 2012 Jiri (George) Lebl <jirka 5z com>

[George Lebl <jirka src gnome org>]

commit 137987632e967829fc53c253677b9948c1f6ace1
Author: Jiri (George) Lebl <jirka 5z com>
Date:   Tue Nov 27 12:59:26 2012 -0600

    Tue Nov 27 12:59:21 2012  Jiri (George) Lebl <jirka 5z com>
    
    	* configure.in: Require MPFR 2.3.0
    
    	* src/funclib.c: add Bessel functions
    	  (BesselJ0, BesselJ1, BesselJn, BesselY0, BesselY1, BesselYn)

 ChangeLog     |    7 ++
 configure.in  |    4 +-
 src/funclib.c |  271 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 280 insertions(+), 2 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index d1426c6..2bd2687 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+Tue Nov 27 12:59:21 2012  Jiri (George) Lebl <jirka 5z com>
+
+	* configure.in: Require MPFR 2.3.0
+
+	* src/funclib.c: add Bessel functions
+	  (BesselJ0, BesselJ1, BesselJn, BesselY0, BesselY1, BesselYn)
+
 Tue Nov 27 10:05:50 2012  Jiri (George) Lebl <jirka 5z com>
 
 	* atlapack/ atblas/: remove, this was just a pain and I don't have
diff --git a/configure.in b/configure.in
index 8afb1ae..a490cf7 100644
--- a/configure.in
+++ b/configure.in
@@ -80,8 +80,8 @@ AC_CHECK_LIB(gmp, __gmpz_init_set,[GMP_LIB=-lgmp], [
 		AC_MSG_ERROR(GMP Library not found))])
 AC_SUBST(GMP_LIB)
 
-AC_CHECK_LIB(mpfr, mpfr_pow_z, [],
-	[AC_MSG_ERROR([New enough MPFR (2.2.0+) not found, see http://www.mpfr.org])],
+AC_CHECK_LIB(mpfr, mpfr_j0, [],
+	[AC_MSG_ERROR([New enough MPFR (2.3.0+) not found, see http://www.mpfr.org])],
 	"$GMP_LIB")
 
 # No schemas yet so ignore this section for now
diff --git a/src/funclib.c b/src/funclib.c
index 0b7d3d6..845ff7a 100644
--- a/src/funclib.c
+++ b/src/funclib.c
@@ -68,6 +68,12 @@ static GelEFunc *Im_function = NULL;
 /*static GelEFunc *ErrorFunction_function = NULL;*/
 static GelEFunc *RiemannZeta_function = NULL;
 static GelEFunc *GammaFunction_function = NULL;
+static GelEFunc *BesselJ0_function = NULL;
+static GelEFunc *BesselJ1_function = NULL;
+static GelEFunc *BesselY0_function = NULL;
+static GelEFunc *BesselY1_function = NULL;
+/*static GelEFunc *BesselJn_function = NULL;
+static GelEFunc *BesselYn_function = NULL;*/
 static GelEFunc *pi_function = NULL;
 static GelEFunc *e_function = NULL;
 static GelEFunc *GoldenRatio_function = NULL;
@@ -1538,6 +1544,253 @@ GammaFunction_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
 	return gel_makenum (retw);
 }
 
+static GelETree *
+BesselJ0_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+	mpfr_ptr num;
+	mpfr_t tmp;
+	mpfr_t ret;
+	mpw_t retw;
+
+	if (a[0]->type == GEL_FUNCTION_NODE ||
+	    a[0]->type == GEL_IDENTIFIER_NODE) {
+		return gel_function_from_function (BesselJ0_function, a[0]);
+	}
+
+	if (a[0]->type == GEL_MATRIX_NODE)
+		return gel_apply_func_to_matrix (ctx, a[0], BesselJ0_op, "BesselJ0", exception);
+
+	if G_UNLIKELY ( ! check_argument_number (a, 0, "BesselJ0"))
+		return NULL;
+	if G_UNLIKELY (mpw_is_complex (a[0]->val.value)) {
+		gel_errorout (_("%s: Not implemented (yet) for complex values"),
+			      "BesselJ0");
+		return NULL;
+	}
+
+	MPW_MPF_REAL (num, a[0]->val.value, tmp);
+
+	mpfr_init (ret);
+	mpfr_j0 (ret, num, GMP_RNDN);
+
+	MPW_MPF_KILL (num, tmp);
+
+	mpw_init (retw);
+	mpw_set_mpf_use (retw, ret);
+
+	return gel_makenum (retw);
+}
+
+static GelETree *
+BesselJ1_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+	mpfr_ptr num;
+	mpfr_t tmp;
+	mpfr_t ret;
+	mpw_t retw;
+
+	if (a[0]->type == GEL_FUNCTION_NODE ||
+	    a[0]->type == GEL_IDENTIFIER_NODE) {
+		return gel_function_from_function (BesselJ1_function, a[0]);
+	}
+
+	if (a[0]->type == GEL_MATRIX_NODE)
+		return gel_apply_func_to_matrix (ctx, a[0], BesselJ1_op, "BesselJ1", exception);
+
+	if G_UNLIKELY ( ! check_argument_number (a, 0, "BesselJ1"))
+		return NULL;
+	if G_UNLIKELY (mpw_is_complex (a[0]->val.value)) {
+		gel_errorout (_("%s: Not implemented (yet) for complex values"),
+			      "BesselJ1");
+		return NULL;
+	}
+
+	MPW_MPF_REAL (num, a[0]->val.value, tmp);
+
+	mpfr_init (ret);
+	mpfr_j1 (ret, num, GMP_RNDN);
+
+	MPW_MPF_KILL (num, tmp);
+
+	mpw_init (retw);
+	mpw_set_mpf_use (retw, ret);
+
+	return gel_makenum (retw);
+}
+
+/* FIXME: implement over matrices / functions */
+static GelETree *
+BesselJn_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+	mpfr_ptr num;
+	mpfr_t tmp;
+	mpfr_t ret;
+	mpw_t retw;
+	long n;
+
+	if (a[0]->type == GEL_MATRIX_NODE)
+		return gel_apply_func_to_matrix (ctx, a[0], BesselJn_op, "BesselJn", exception);
+
+	if G_UNLIKELY ( ! check_argument_integer (a, 0, "BesselJn"))
+		return NULL;
+	n = mpw_get_long(a[0]->val.value);
+	if G_UNLIKELY (gel_error_num != 0) {
+		gel_error_num = 0;
+		return NULL;
+	}
+
+	if G_UNLIKELY ( ! check_argument_number (a, 1, "BesselJn"))
+		return NULL;
+	if G_UNLIKELY (mpw_is_complex (a[1]->val.value)) {
+		gel_errorout (_("%s: Not implemented (yet) for complex values"),
+			      "BesselJn");
+		return NULL;
+	}
+
+	MPW_MPF_REAL (num, a[1]->val.value, tmp);
+
+	mpfr_init (ret);
+	mpfr_jn (ret, n, num, GMP_RNDN);
+
+	MPW_MPF_KILL (num, tmp);
+
+	mpw_init (retw);
+	mpw_set_mpf_use (retw, ret);
+
+	return gel_makenum (retw);
+}
+
+static GelETree *
+BesselY0_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+	mpfr_ptr num;
+	mpfr_t tmp;
+	mpfr_t ret;
+	mpw_t retw;
+
+	if (a[0]->type == GEL_FUNCTION_NODE ||
+	    a[0]->type == GEL_IDENTIFIER_NODE) {
+		return gel_function_from_function (BesselY0_function, a[0]);
+	}
+
+	if (a[0]->type == GEL_MATRIX_NODE)
+		return gel_apply_func_to_matrix (ctx, a[0], BesselY0_op, "BesselY0", exception);
+
+	if G_UNLIKELY ( ! check_argument_number (a, 0, "BesselY0"))
+		return NULL;
+	if G_UNLIKELY (mpw_is_complex (a[0]->val.value)) {
+		gel_errorout (_("%s: Not implemented (yet) for complex values"),
+			      "BesselY0");
+		return NULL;
+	}
+	if G_UNLIKELY (mpw_sgn (a[0]->val.value) <= 0) {
+		gel_errorout (_("%s: Bessel functions of second kind not defined for nonpositive real numbers"),
+			      "BesselY0");
+		return NULL;
+	}
+
+	MPW_MPF_REAL (num, a[0]->val.value, tmp);
+
+	mpfr_init (ret);
+	mpfr_y0 (ret, num, GMP_RNDN);
+
+	MPW_MPF_KILL (num, tmp);
+
+	mpw_init (retw);
+	mpw_set_mpf_use (retw, ret);
+
+	return gel_makenum (retw);
+}
+
+static GelETree *
+BesselY1_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+	mpfr_ptr num;
+	mpfr_t tmp;
+	mpfr_t ret;
+	mpw_t retw;
+
+	if (a[0]->type == GEL_FUNCTION_NODE ||
+	    a[0]->type == GEL_IDENTIFIER_NODE) {
+		return gel_function_from_function (BesselY1_function, a[0]);
+	}
+
+	if (a[0]->type == GEL_MATRIX_NODE)
+		return gel_apply_func_to_matrix (ctx, a[0], BesselY1_op, "BesselY1", exception);
+
+	if G_UNLIKELY ( ! check_argument_number (a, 0, "BesselY1"))
+		return NULL;
+	if G_UNLIKELY (mpw_is_complex (a[0]->val.value)) {
+		gel_errorout (_("%s: Not implemented (yet) for complex values"),
+			      "BesselY1");
+		return NULL;
+	}
+	if G_UNLIKELY (mpw_sgn (a[0]->val.value) <= 0) {
+		gel_errorout (_("%s: Bessel functions of second kind not defined for nonpositive real numbers"),
+			      "BesselY1");
+		return NULL;
+	}
+
+	MPW_MPF_REAL (num, a[0]->val.value, tmp);
+
+	mpfr_init (ret);
+	mpfr_y1 (ret, num, GMP_RNDN);
+
+	MPW_MPF_KILL (num, tmp);
+
+	mpw_init (retw);
+	mpw_set_mpf_use (retw, ret);
+
+	return gel_makenum (retw);
+}
+
+/* FIXME: implement over matrices / functions */
+static GelETree *
+BesselYn_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+	mpfr_ptr num;
+	mpfr_t tmp;
+	mpfr_t ret;
+	mpw_t retw;
+	long n;
+
+	if (a[0]->type == GEL_MATRIX_NODE)
+		return gel_apply_func_to_matrix (ctx, a[0], BesselYn_op, "BesselYn", exception);
+
+	if G_UNLIKELY ( ! check_argument_integer (a, 0, "BesselYn"))
+		return NULL;
+	n = mpw_get_long(a[0]->val.value);
+	if G_UNLIKELY (gel_error_num != 0) {
+		gel_error_num = 0;
+		return NULL;
+	}
+
+	if G_UNLIKELY ( ! check_argument_number (a, 1, "BesselYn"))
+		return NULL;
+	if G_UNLIKELY (mpw_is_complex (a[1]->val.value)) {
+		gel_errorout (_("%s: Not implemented (yet) for complex values"),
+			      "BesselYn");
+		return NULL;
+	}
+	if G_UNLIKELY (mpw_sgn (a[1]->val.value) <= 0) {
+		gel_errorout (_("%s: Bessel functions of second kind not defined for nonpositive real numbers"),
+			      "BesselYn");
+		return NULL;
+	}
+
+	MPW_MPF_REAL (num, a[1]->val.value, tmp);
+
+	mpfr_init (ret);
+	mpfr_yn (ret, n, num, GMP_RNDN);
+
+	MPW_MPF_KILL (num, tmp);
+
+	mpw_init (retw);
+	mpw_set_mpf_use (retw, ret);
+
+	return gel_makenum (retw);
+}
+
 
 static GelETree *
 IsNull_op(GelCtx *ctx, GelETree * * a, gboolean *exception)
@@ -6468,6 +6721,24 @@ gel_funclib_addall(void)
 	GammaFunction_function = f;
 	ALIAS (Gamma, 1, GammaFunction);
 
+	FUNC (BesselJ0, 1, "x", "functions", N_("The Bessel function of first kind of order 0"));
+	f->no_mod_all_args = 1;
+	BesselJ0_function = f;
+	FUNC (BesselJ1, 1, "x", "functions", N_("The Bessel function of first kind of order 1"));
+	f->no_mod_all_args = 1;
+	BesselJ1_function = f;
+	FUNC (BesselJn, 2, "n,x", "functions", N_("The Bessel function of first kind of order n"));
+	f->no_mod_all_args = 1;
+
+	FUNC (BesselY0, 1, "x", "functions", N_("The Bessel function of second kind of order 0"));
+	f->no_mod_all_args = 1;
+	BesselJ0_function = f;
+	FUNC (BesselY1, 1, "x", "functions", N_("The Bessel function of second kind of order 1"));
+	f->no_mod_all_args = 1;
+	BesselJ1_function = f;
+	FUNC (BesselYn, 2, "n,x", "functions", N_("The Bessel function of second kind of integer order n"));
+	f->no_mod_all_args = 1;
+
 	FUNC (sqrt, 1, "x", "numeric", N_("The square root"));
 	f->propagate_mod = 1;
 	sqrt_function = f;