===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/bfaux.c,v
retrieving revision 1.2
retrieving revision 1.13
diff -u -p -r1.2 -r1.13
--- OpenXM_contrib2/asir2000/builtin/bfaux.c	2015/08/05 01:10:38	1.2
+++ OpenXM_contrib2/asir2000/builtin/bfaux.c	2017/03/09 00:46:44	1.13
@@ -1,16 +1,45 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.1 2015/08/04 06:55:02 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.12 2016/03/14 04:15:05 noro Exp $ */
 #include "ca.h"
 #include "parse.h"
 
-void Peval(), Psetprec(), Ptodouble();
+void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
+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();
+void Prk_ratmat();
 
 struct ftab bf_tab[] = {
 	{"eval",Peval,-2},
 	{"setprec",Psetprec,-1},
+	{"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},
+	{"rk_ratmat",Prk_ratmat,7},
 	{0,0,0},
 };
 
+int mpfr_roundmode = MPFR_RNDN;
+
 void Ptodouble(NODE arg,Num *rp)
 {
 	double r,i;
@@ -42,15 +71,13 @@ void Ptodouble(NODE arg,Num *rp)
 	}
 }
 
-void Peval(arg,rp)
-NODE arg;
-Obj *rp;
+void Peval(NODE arg,Obj *rp)
 {
   int prec;
 
 	asir_assert(ARG0(arg),O_R,"eval");
   if ( argc(arg) == 2 ) {
-	  prec = QTOS((Q)ARG1(arg));
+	  prec = QTOS((Q)ARG1(arg))*3.32193;
     if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
     else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
   } else
@@ -58,18 +85,41 @@ Obj *rp;
 	evalr(CO,(Obj)ARG0(arg),prec,rp);
 }
 
-/* bit precision */
+/* set/get decimal precision */
 
 void Psetprec(NODE arg,Obj *rp)
 {
 	int p;
 	Q q;
+	int prec,dprec;
+  
+  prec = mpfr_get_default_prec();
+  /* decimal precision */
+  dprec = prec*0.30103;
+  STOQ(dprec,q); *rp = (Obj)q;
+	if ( arg ) {
+		asir_assert(ARG0(arg),O_N,"setprec");
+		p = QTOS((Q)ARG0(arg))*3.32193;
+		if ( p > 0 )
+			prec = p;
+	}
+  if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
+  else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
+	mpfr_set_default_prec(prec);
+}
 
-	long prec = mpfr_get_default_prec();
+/* set/get bit precision */
 
-    STOQ(prec,q); *rp = (Obj)q;
+void Psetbprec(NODE arg,Obj *rp)
+{
+	int p;
+	Q q;
+	int prec;
+  
+  prec = mpfr_get_default_prec();
+  STOQ(prec,q); *rp = (Obj)q;
 	if ( arg ) {
-		asir_assert(ARG0(arg),O_N,"setprec");
+		asir_assert(ARG0(arg),O_N,"setbprec");
 		p = QTOS((Q)ARG0(arg));
 		if ( p > 0 )
 			prec = p;
@@ -79,23 +129,61 @@ void Psetprec(NODE arg,Obj *rp)
 	mpfr_set_default_prec(prec);
 }
 
+void Psetround(NODE arg,Q *rp)
+{
+  int round;
+
+  STOQ(mpfr_roundmode,*rp);
+  if ( arg ) {
+		asir_assert(ARG0(arg),O_N,"setround");
+    round = QTOS((Q)ARG0(arg)); 
+    switch ( round ) {
+    case 0:
+      mpfr_roundmode = MPFR_RNDN;
+      break;
+    case 1:
+      mpfr_roundmode = MPFR_RNDZ;
+      break;
+    case 2:
+      mpfr_roundmode = MPFR_RNDU;
+      break;
+    case 3:
+      mpfr_roundmode = MPFR_RNDD;
+      break;
+    case 4:
+      mpfr_roundmode = MPFR_RNDA;
+      break;
+    case 5:
+      mpfr_roundmode = MPFR_RNDF;
+      break;
+    case 6:
+      mpfr_roundmode = MPFR_RNDNA;
+      break;
+    default:
+      error("setround : invalid rounding mode");
+      break;
+    }
+  }
+}
+
 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_RNDN);
-    *rp = r; 
+	mpfr_const_pi(r->body,mpfr_roundmode);
+  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;
 
@@ -103,214 +191,534 @@ void mp_e(NODE arg,BF *rp)
 	NEWBF(r);
 	prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
 	mpfr_init(one);
-	mpfr_set_ui(one,1,MPFR_RNDN);
-	mpfr_exp(r->body,one,MPFR_RNDN);
-    *rp = r; 
+	mpfr_set_ui(one,1,mpfr_roundmode);
+	mpfr_exp(r->body,one,mpfr_roundmode);
+  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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_RNDN);
-    *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_ceil(NODE arg,Q *rp)
+{
+	Num a;
+  int prec;
+	BF r;
+	mpz_t t;
+	GZ rz;
+
+  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); 
+}
+
+double **almat_double(int n)
+{
+  int i;
+  double **a;
+
+  a = (double **)MALLOC(n*sizeof(double *));
+  for ( i = 0; i < n; i++ )
+    a[i] = (double *)MALLOC(n*sizeof(double));
+  return a;
+}
+
+/*
+ *  k <- (A(xi)-(sbeta-mn2/xi))f
+ *  A(t) = (num[0]+num[1]t+...+num[d-1]*t^(d-1))/den(t)
+ */
+
+struct jv {
+  int j;
+  double v;
+};
+
+struct smat {
+  int *rlen;
+  struct jv **row;
+};
+
+void eval_pfaffian2(double *k,int n,int d,struct smat *num,P den,double xi,double *f)
+{
+  struct smat ma;
+  struct jv *maj;
+  int i,j,l,s;
+  double t,dn;
+  P r;
+  Real u;
+
+  memset(k,0,n*sizeof(double));
+  for ( i = d-1; i >= 0; i-- ) {
+    ma = num[i];
+    for ( j = 0; j < n; j++ ) {
+      maj = ma.row[j];
+      l = ma.rlen[j];
+      for ( t = 0, s = 0; s < l; s++, maj++ ) t += maj->v*f[maj->j];
+      k[j] = k[j]*xi+t;
+    }
+  }
+  MKReal(xi,u);
+  substp(CO,den,den->v,(P)u,&r); dn = ToReal(r);
+  for ( j = 0; j < n; j++ )
+    k[j] /= dn;
+}
+
+void Prk_ratmat(NODE arg,LIST *rp)
+{
+  VECT mat;
+  P den;
+  int ord;
+  double sbeta,x0,x1,xi,h,mn2,hd;
+  double a2,a3,a4,a5,a6;
+  double b21,b31,b32,b41,b42,b43,b51,b52,b53,b54,b61,b62,b63,b64,b65;
+  double c1,c2,c3,c4,c5,c6,c7;
+  VECT fv;
+  int step,j,i,k,n,d,len,s;
+  struct smat *num;
+  Obj **b;
+  MAT mati;
+  double *f,*w,*k1,*k2,*k3,*k4,*k5,*k6;
+  NODE nd,nd1;
+  Real x,t;
+  LIST l;
+
+  ord = QTOS((Q)ARG0(arg));
+  mat = (VECT)ARG1(arg); den = (P)ARG2(arg);
+  x0 = ToReal((Num)ARG3(arg)); x1 = ToReal((Num)ARG4(arg));
+  step = QTOS((Q)ARG5(arg)); fv = (VECT)ARG6(arg);
+  h = (x1-x0)/step;
+
+  n = fv->len;
+  d = mat->len;
+  num = (struct smat *)MALLOC(d*sizeof(struct smat)); 
+  for ( i = 0; i < d; i++ ) {
+    num[i].rlen = (int *)MALLOC_ATOMIC(n*sizeof(int));
+    num[i].row = (struct jv **)MALLOC(n*sizeof(struct jv *));
+    mati = (MAT)mat->body[i];
+    b = (Obj **)mati->body;
+    for ( j = 0; j < n; j++ ) {
+      for ( len = k = 0; k < n; k++ )
+        if ( b[j][k] ) len++;
+      num[i].rlen[j] = len;
+      if ( !len )
+        num[i].row[j] = 0;
+      else {
+        num[i].row[j] = (struct jv *)MALLOC_ATOMIC((len)*sizeof(struct jv));
+        for ( s = k = 0; k < n; k++ )
+          if ( b[j][k] ) {
+             num[i].row[j][s].j = k;
+             num[i].row[j][s].v = ToReal((Num)b[j][k]);
+             s++;
+          }
+      }
+    }
+  }
+  f = (double *)MALLOC_ATOMIC(n*sizeof(double));
+  for ( j = 0; j < n; j++ )
+    f[j] = ToReal((Num)fv->body[j]);
+  w = (double *)MALLOC_ATOMIC(n*sizeof(double));
+  k1 = (double *)MALLOC_ATOMIC(n*sizeof(double));
+  k2 = (double *)MALLOC_ATOMIC(n*sizeof(double));
+  k3 = (double *)MALLOC_ATOMIC(n*sizeof(double));
+  k4 = (double *)MALLOC_ATOMIC(n*sizeof(double));
+  k5 = (double *)MALLOC_ATOMIC(n*sizeof(double));
+  k6 = (double *)MALLOC(n*sizeof(double));
+  nd = 0;
+  switch ( ord ) {
+  case 4:
+    a2 = 1/2.0*h; b21 = 1/2.0*h;
+    a3 = 1/2.0*h; b31 = 0.0;   b32 = 1/2.0*h;
+    a4 = 1.0*h;   b41 = 0.0;   b42 = 0.0;    b43 = 1.0*h;
+    c1 = 1/6.0*h; c2 = 1/3.0*h;     c3 =  1/3.0*h; c4 = 1/6.0*h;
+    for ( i = 0; i < step; i++ ) {
+      if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
+      xi = x0+i*h;
+      eval_pfaffian2(k1,n,d,num,den,xi,f);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
+      eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
+      eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
+      eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += c1*k1[j]+c2*k2[j]+c3*k3[j]+c4*k4[j];
+      memcpy(f,w,n*sizeof(double));
+      MKReal(f[0],t);
+      MKReal(xi+h,x);
+      nd1 = mknode(2,x,t);
+      MKLIST(l,nd1);
+      MKNODE(nd1,l,nd);
+      nd = nd1;
+      for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
+    }
+    MKLIST(*rp,nd);
+    break;
+  case 5:
+  default:
+    a2 = 1/4.0*h; b21 = 1/4.0*h;
+    a3 = 1/4.0*h; b31 = 1/8.0*h; b32 = 1/8.0*h;
+    a4 = 1/2.0*h; b41 = 0.0;   b42 = 0.0;    b43 = 1/2.0*h;
+    a5 = 3/4.0*h; b51 = 3/16.0*h;b52 = -3/8.0*h; b53 = 3/8.0*h;   b54 = 9/16.0*h;
+    a6 = 1.0*h;   b61 = -3/7.0*h;b62 = 8/7.0*h;  b63 = 6/7.0*h;   b64 = -12/7.0*h; b65 = 8/7.0*h;
+    c1 = 7/90.0*h; c2 = 0.0;     c3 =  16/45.0*h; c4 = 2/15.0*h;   c5 = 16/45.0*h; c6 = 7/90.0*h;
+    for ( i = 0; i < step; i++ ) {
+      if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
+      xi = x0+i*h;
+      eval_pfaffian2(k1,n,d,num,den,xi,f);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
+      eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
+      eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
+      eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b51*k1[j]+b52*k2[j]+b53*k3[j]+b54*k4[j];
+      eval_pfaffian2(k5,n,d,num,den,xi+a5,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b61*k1[j]+b62*k2[j]+b63*k3[j]+b64*k4[j]+b65*k5[j];
+      eval_pfaffian2(k6,n,d,num,den,xi+a6,w);
+      memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += c1*k1[j]+c2*k2[j]+c3*k3[j]+c4*k4[j]+c5*k5[j]+c6*k6[j];
+      memcpy(f,w,n*sizeof(double));
+      MKReal(f[0],t);
+      MKReal(xi+h,x);
+      nd1 = mknode(2,x,t);
+      MKLIST(l,nd1);
+      MKNODE(nd1,l,nd);
+      nd = nd1;
+      for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
+    }
+    MKLIST(*rp,nd);
+    break;
+  }
 }