===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/bfaux.c,v
retrieving revision 1.5
retrieving revision 1.10
diff -u -p -r1.5 -r1.10
--- OpenXM_contrib2/asir2000/builtin/bfaux.c	2015/08/07 05:30:35	1.5
+++ OpenXM_contrib2/asir2000/builtin/bfaux.c	2015/08/20 08:42:07	1.10
@@ -1,9 +1,15 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.4 2015/08/06 23:41:52 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.9 2015/08/17 05:18:36 noro Exp $ */
 #include "ca.h"
 #include "parse.h"
 
 void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
-void Pmpfr_gamma();
+void Pmpfr_ai();
+void Pmpfr_eint(), Pmpfr_erf(), Pmpfr_erfc(), Pmpfr_li2();
+void Pmpfr_zeta();
+void Pmpfr_j0(), Pmpfr_j1();
+void Pmpfr_y0(), Pmpfr_y1();
+void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma();
+void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil();
 
 struct ftab bf_tab[] = {
 	{"eval",Peval,-2},
@@ -11,7 +17,22 @@ struct ftab bf_tab[] = {
 	{"setbprec",Psetbprec,-1},
 	{"setround",Psetround,-1},
 	{"todouble",Ptodouble,1},
+	{"mpfr_ai",Pmpfr_ai,-2},
+	{"mpfr_zeta",Pmpfr_zeta,-2},
+	{"mpfr_j0",Pmpfr_j0,-2},
+	{"mpfr_j1",Pmpfr_j1,-2},
+	{"mpfr_y0",Pmpfr_y0,-2},
+	{"mpfr_y1",Pmpfr_y1,-2},
+	{"mpfr_eint",Pmpfr_eint,-2},
+	{"mpfr_erf",Pmpfr_erf,-2},
+	{"mpfr_erfc",Pmpfr_erfc,-2},
+	{"mpfr_li2",Pmpfr_li2,-2},
 	{"mpfr_gamma",Pmpfr_gamma,-2},
+	{"mpfr_lngamma",Pmpfr_gamma,-2},
+	{"mpfr_digamma",Pmpfr_gamma,-2},
+	{"mpfr_floor",Pmpfr_floor,-2},
+	{"mpfr_ceil",Pmpfr_ceil,-2},
+	{"mpfr_round",Pmpfr_round,-2},
 	{0,0,0},
 };
 
@@ -147,19 +168,20 @@ Num tobf(Num a,int prec);
 
 void mp_pi(NODE arg,BF *rp)
 {
-    int prec;
+  int prec;
 	BF r;
 
 	prec = arg ? QTOS((Q)ARG0(arg)) : 0;
 	NEWBF(r);
 	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
 	mpfr_const_pi(r->body,mpfr_roundmode);
-    *rp = r; 
+  if ( !cmpbf((Num)r,0) ) r = 0;
+  *rp = r; 
 }
 
 void mp_e(NODE arg,BF *rp)
 {
-    int prec;
+  int prec;
 	mpfr_t one;
 	BF r;
 
@@ -169,227 +191,355 @@ void mp_e(NODE arg,BF *rp)
 	mpfr_init(one);
 	mpfr_set_ui(one,1,mpfr_roundmode);
 	mpfr_exp(r->body,one,mpfr_roundmode);
-    *rp = r; 
+  if ( !cmpbf((Num)r,0) ) r = 0;
+  *rp = r; 
 }
 
-void mp_sin(NODE arg,BF *rp)
+void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp)
 {
 	Num a;
-    int prec;
-	BF r;
+  int prec;
+	BF r,re,im;
+  C c;
+  mpc_t mpc,a1;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
+	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
 	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_sin(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  if ( a && NID(a)==N_C ) {
+    mpc_init2(mpc,prec); mpc_init2(a1,prec);
+    re = (BF)((C)a)->r; im = (BF)((C)a)->i;
+    mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
+    (*mpc_f)(mpc,a1,mpfr_roundmode);
+    MPFRTOBF(mpc_realref(mpc),re);
+    MPFRTOBF(mpc_imagref(mpc),im);
+    if ( !cmpbf((Num)re,0) ) re = 0;
+    if ( !cmpbf((Num)im,0) ) im = 0;
+    if ( !im )
+      *rp = (Num)re;
+    else {
+      NEWC(c); c->r = (Num)re; c->i = (Num)im;
+      *rp = (Num)c; 
+    }
+  } else {
+	  NEWBF(r);
+	  mpfr_init2(r->body,prec);
+	  (*mpfr_f)(r->body,((BF)a)->body,mpfr_roundmode);
+    if ( !cmpbf((Num)r,0) ) r = 0;
+    *rp = (Num)r; 
+  }
 }
 
-void mp_cos(NODE arg,BF *rp)
+void mp_sin(NODE arg,Num *rp)
 {
+  mpfr_or_mpc(arg,mpfr_sin,mpc_sin,rp);
+}
+
+void mp_cos(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_cos,mpc_cos,rp);
+}
+
+void mp_tan(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_tan,mpc_tan,rp);
+}
+
+void mp_asin(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_asin,mpc_asin,rp);
+}
+
+void mp_acos(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_acos,mpc_acos,rp);
+}
+
+void mp_atan(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_atan,mpc_atan,rp);
+}
+
+void mp_sinh(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_sinh,mpc_sinh,rp);
+}
+
+void mp_cosh(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_cosh,mpc_cosh,rp);
+}
+
+void mp_tanh(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_tanh,mpc_tanh,rp);
+}
+
+void mp_asinh(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_asinh,mpc_asinh,rp);
+}
+
+void mp_acosh(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_acosh,mpc_acosh,rp);
+}
+
+void mp_atanh(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_atanh,mpc_atanh,rp);
+}
+
+void mp_exp(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_exp,mpc_exp,rp);
+}
+
+void mp_log(NODE arg,Num *rp)
+{
+  mpfr_or_mpc(arg,mpfr_log,mpc_log,rp);
+}
+
+void mp_pow(NODE arg,Num *rp)
+{
+	Num a,e;
+  int prec;
+	BF r,re,im;
+  C c;
+  mpc_t mpc,a1,e1;
+
+	prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
+	a = tobf(ARG0(arg),prec);
+	e = tobf(ARG1(arg),prec);
+  if ( NID(a) == N_C || NID(e) == N_C || MPFR_SIGN(((BF)a)->body) < 0 ) {
+    mpc_init2(mpc,prec); mpc_init2(a1,prec); mpc_init2(e1,prec);
+    if ( NID(a) == N_C ) {
+      re = (BF)((C)a)->r; im = (BF)((C)a)->i;
+      mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
+    } else {
+      re = (BF)a;
+      mpc_set_fr(a1,re->body,mpfr_roundmode);
+    }
+    if ( NID(e) == N_C ) {
+      re = (BF)((C)e)->r; im = (BF)((C)e)->i;
+      mpc_set_fr_fr(e1,re->body,im->body,mpfr_roundmode);
+    } else {
+      re = (BF)e;
+      mpc_set_fr(e1,re->body,mpfr_roundmode);
+    }
+    mpc_pow(mpc,a1,e1,mpfr_roundmode);
+    MPFRTOBF(mpc_realref(mpc),re);
+    MPFRTOBF(mpc_imagref(mpc),im);
+    if ( !cmpbf((Num)re,0) ) re = 0;
+    if ( !cmpbf((Num)im,0) ) im = 0;
+    if ( !im )
+      *rp = (Num)re;
+    else {
+      NEWC(c); c->r = (Num)re; c->i = (Num)im;
+      *rp = (Num)c; 
+    }
+  } else {
+	  NEWBF(r);
+	  mpfr_init2(r->body,prec);
+	  mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
+    *rp = (Num)r; 
+  }
+}
+
+#define SETPREC \
+ (prec)=NEXT(arg)?QTOS((Q)ARG1(arg)):0;\
+ (prec)*=3.32193;\
+ (a)=tobf(ARG0(arg),prec);\
+ NEWBF(r);\
+ prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
+
+
+void Pmpfr_gamma(NODE arg,BF *rp)
+{
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_cos(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_tan(NODE arg,BF *rp)
+void Pmpfr_lngamma(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_tan(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_asin(NODE arg,BF *rp)
+void Pmpfr_digamma(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_asin(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
-void mp_acos(NODE arg,BF *rp)
+
+void Pmpfr_zeta(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_acos(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
-void mp_atan(NODE arg,BF *rp)
+
+void Pmpfr_eint(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_atan(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_sinh(NODE arg,BF *rp)
+void Pmpfr_erf(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_sinh(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_cosh(NODE arg,BF *rp)
+void Pmpfr_erfc(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_cosh(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_tanh(NODE arg,BF *rp)
+void Pmpfr_j0(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_tanh(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_asinh(NODE arg,BF *rp)
+void Pmpfr_j1(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_asinh(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
-void mp_acosh(NODE arg,BF *rp)
+
+void Pmpfr_y0(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_acosh(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
-void mp_atanh(NODE arg,BF *rp)
+
+void Pmpfr_y1(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_atanh(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_exp(NODE arg,BF *rp)
+void Pmpfr_li2(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_exp(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_log(NODE arg,BF *rp)
+void Pmpfr_ai(NODE arg,BF *rp)
 {
 	Num a;
-    int prec;
+  int prec;
 	BF r;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_log(r->body,((BF)a)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
+  *rp = r; 
 }
 
-void mp_pow(NODE arg,BF *rp)
+void Pmpfr_floor(NODE arg,Q *rp)
 {
-	Num a,e;
-    int prec;
+	Num a;
+  int prec;
 	BF r;
+	mpz_t t;
+	GZ rz;
 
-	prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : 0;
-	a = tobf(ARG0(arg),prec);
-	e = tobf(ARG1(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
-    *rp = r; 
+  SETPREC
+	mpfr_floor(r->body,((BF)a)->body);
+	mpz_init(t);
+	mpfr_get_z(t,r->body,mpfr_roundmode);
+	MPZTOGZ(t,rz);
+	*rp = gztoz(rz); 
 }
 
-void Pmpfr_gamma(NODE arg,BF *rp)
+void Pmpfr_ceil(NODE arg,Q *rp)
 {
 	Num a;
   int prec;
 	BF r;
+	mpz_t t;
+	GZ rz;
 
-	prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
-  prec *= 3.32193;
-	a = tobf(ARG0(arg),prec);
-	NEWBF(r);
-	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
-	mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
-  *rp = r; 
+  SETPREC
+	mpfr_ceil(r->body,((BF)a)->body);
+	mpz_init(t);
+	mpfr_get_z(t,r->body,mpfr_roundmode);
+	MPZTOGZ(t,rz);
+	*rp = gztoz(rz); 
+}
+
+void Pmpfr_round(NODE arg,Q *rp)
+{
+	Num a;
+  int prec;
+	BF r;
+	mpz_t t;
+	GZ rz;
+
+  SETPREC
+	mpfr_round(r->body,((BF)a)->body);
+	mpz_init(t);
+	mpfr_get_z(t,r->body,mpfr_roundmode);
+	MPZTOGZ(t,rz);
+	*rp = gztoz(rz); 
 }