===================================================================
RCS file: /home/cvs/OpenXM/src/hgm/mh/src/jack-n.c,v
retrieving revision 1.14
retrieving revision 1.46
diff -u -p -r1.14 -r1.46
--- OpenXM/src/hgm/mh/src/jack-n.c	2014/02/24 07:58:05	1.14
+++ OpenXM/src/hgm/mh/src/jack-n.c	2016/02/15 06:02:39	1.46
@@ -4,8 +4,10 @@
 #include <math.h>
 #include <string.h>
 #include "sfile.h"
+
+#define VSTRING "%!version2.0"
 /*
-  $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.13 2013/03/07 05:23:31 takayama Exp $
+  $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.45 2016/02/13 22:56:50 takayama Exp $
   Ref: copied from this11/misc-2011/A1/wishart/Prog
   jack-n.c, translated from mh.rr or tk_jack.rr in the asir-contrib. License: LGPL
   Koev-Edelman for higher order derivatives.
@@ -15,6 +17,13 @@
   2. Use the recurrence to obtain beta().
   3. Easier input data file format for mh-n.c
   Changelog:
+  2016.02.15 log of Ef
+  2016.02.09 unify 2F1 and 1F1. Parser.
+  2016.02.04 Ef_type (exponential or scalar factor type)
+  2016.02.01 ifdef C_2F1 ...
+  2016.01.12 2F1
+  2014.03.15 http://fe.math.kobe-u.ac.jp/Movies/oxvh/2014-03-11-jack-n-c-automatic  see also hgm/doc/ref.html, @s/2014/03/15-my-note-jack-automatic-order-F_A-casta.pdf.
+  2014.03.14, --automatic option. Output estimation data.
   2012.02.21, porting to OpenXM/src/hgm
   2011.12.22, --table option, which is experimental.
   2011.12.19, bug fix.  jjk() should return double. It can become more than max int.
@@ -53,8 +62,14 @@ static double Ef2; 
 #define M_n0 3 /* used for tests. Must be equal to M_n */
 #define M_m_MAX 200
 #define M_nmx  M_m_MAX  /* maximal of M_n */
-#define A_LEN  1 /* (a_1) , (a_1, ..., a_p)*/
-#define B_LEN  1 /* (b_1) */
+
+static int A_LEN=-1;  /* (a_1, ..., a_p), A_LEN=p */
+static int B_LEN=-1;  /* (b_1,..., b_q),  B_LEN=q */
+static double *A_pFq=NULL;
+static double *B_pFq=NULL;
+static int Ef_type=1;  /* 1F1 for Wishart */
+/* Ef_type=2;   2F1, Hashiguchi note (1) */
+
 static int Debug = 0;
 static int Alpha = 2;  /* 2 implies the zonal polynomial */
 static int *Darray = NULL;
@@ -93,6 +108,70 @@ static double Xarray[M_nmx][M_m_MAX];  
 static double *M_qk=NULL;  /* saves pochhammerb */
 static double M_rel_error=0.0; /* relative errors */
 
+/* For automatic degree and X0g setting. */
+/*
+ If automatic == 1, then the series is reevaluated as long as t_success!=1
+ by increasing X0g (evaluation point) and M_m (approx degree);
+ */
+static int M_automatic=1;
+/* Estimated degree bound for series expansion. See mh_t */
+static int M_m_estimated_approx_deg=0;
+/* Let F(i) be the approximation up to degree i. 
+   The i-th series_error is defined
+   by |(F(i)-F(i-1))/F(i-1)|.
+*/
+static double M_series_error;
+/*
+  M_series_error < M_assigend_series_error (A) is required for the 
+  estimated_approx_deg.
+ */
+static double M_assigned_series_error=M_ASSIGNED_SERIES_ERROR_DEFAULT;
+/*
+  Let Ef be the exponential factor ( Ef=(4)/1F1 of [HNTT] )
+  If F(M_m)*Ef < x0value_min (B), the success=0 and X0g is increased. 
+  Note that minimal double is about 2e-308
+ */
+static double M_x0value_min=1e-60;
+/*
+  estimated_X0g is the suggested value of X0g.
+ */
+static double M_estimated_X0g=0.0;
+/*
+ X0g should be less than M_X0g_bound.
+ */
+static double M_X0g_bound = 1e+100;
+/*
+ success is set to 1 when (A) and (B) are satisfied.
+ */
+static int M_mh_t_success=1;
+/*
+  recommended_abserr is the recommended value of the absolute error
+  for the Runge-Kutta method. It is defined as
+  assigend_series_error(standing for significant digits)*Ig[0](initial value)
+ */
+static double M_recommended_abserr;
+/*
+  recommended_relerr is the recommended value of the relative error
+  for the Runge-Kutta method..
+ */
+static double M_recommended_relerr;
+/*
+  max of beta(i)*x/2
+ */
+static double M_beta_i_x_o2_max;
+/*
+  minimum of |beta_i-beta_j|
+ */
+static double M_beta_i_beta_j_min;
+/*
+  Value of matrix hg
+*/
+static double M_mh_t_value;
+/*
+ Show the process of updating degree.
+ */
+int M_show_autosteps=1;
+
 /* prototypes */
 static void *mymalloc(int size);
 static int myfree(void *p);
@@ -104,25 +183,18 @@ static int mysum(int L[]);
 static int plength(int P[]);
 static int plength_t(int P[]);
 static void ptrans(int P[M_nmx],int Pt[]);
-static int test_ptrans();
 static int huk(int K[],int I,int J);
 static int hdk(int K[],int I,int J);
 static double jjk(int K[]);
 static double ppoch(double A,int K[]);
 static double ppoch2(double A,double B,int K[]);
 static double mypower(double x,int n);
-static double qk(int K[],double A[A_LEN],double B[B_LEN]);
+static double qk(int K[],double A[],double B[]);
 static int bb(int N[],int K[],int M[],int I,int J);
 static double beta(int K[],int M[]);
 static int printp(int kappa[]);
-static int printp2(int kappa[]);
-static int test_beta();
 static double q3_10(int K[],int M[],int SK);
-static double q3_5(double A[],double B[],int K[],int I);
-static void mtest4();
-static void mtest4b();
 static int nk(int KK[]);
-static int plength2(int P1[],int P2[]);
 static int myeq(int P1[],int P2[]);
 static int pListPartition(int M,int N);
 static int pListPartition2(int Less,int From,int To, int M);
@@ -138,30 +210,94 @@ static int genDarray2(int M,int N);
 static int isHStrip(int Kap[],int Nu[]);
 static void hsExec_beta(void);
 static int genBeta(int Kap[]);
-static int checkBeta1();
 static int psublen(int Kap[],int Mu[]);
 static int genJack(int M,int N);
-static int checkJack1(int M,int N);
-static int checkJack2(int M,int N);
-static int mtest1b();
 
 static int imypower(int x,int n);
 static int usage();
 static int setParamDefault();
 static int next(struct SFILE *fp,char *s,char *msg);
-static int gopen_file(void);
 static int setParam(char *fname);
 static int showParam(struct SFILE *fp,int fd);
+#ifdef STANDALONE
+static int showParam_v1(struct SFILE *fp,int fd);
+#endif
 static double iv_factor(void);
 static double gammam(double a,int n);
 static double mypower(double a,int n);
 
-double mh_t(double A[A_LEN],double B[B_LEN],int N,int M);
+static double iv_factor_ef_type1(void);
+static double iv_factor_ef_type2(void);
+static double lgammam(double a,int n);
+static double liv_factor_ef_type1(void);
+static double liv_factor_ef_type2(void);
+
+static void setM_x(void);
+static void setM_x_ef_type1(void);
+static void setM_x_ef_type2(void);
+
+#ifdef STANDALONE
+static int test_ptrans();
+static int printp2(int kappa[]);
+static int test_beta();
+static void mtest4();
+static void mtest4b();
+static int plength2(int P1[],int P2[]);
+static int checkBeta1();
+static int checkJack1(int M,int N);
+static int checkJack2(int M,int N);
+static int mtest1b();
+static double q3_5(double A[],double B[],int K[],int I);
+#endif
+
+double mh_t(double A[],double B[],int N,int M);
 double mh_t2(int J);
 struct MH_RESULT *jk_main(int argc,char *argv[]);
+struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree);
 int jk_freeWorkArea();
 int jk_initializeWorkArea();
+static void setA(double a[],int alen);  /* set A_LEN and A_pFq */
+static void setB(double b[],int blen);
 
+static void setA(double a[],int alen) {
+  int i;
+  if (alen < 0) {
+	if (A_pFq != NULL) myfree(A_pFq);
+	A_pFq=NULL; A_LEN=-1;
+	return;
+  }
+  if (alen == 0) {
+	A_LEN=0; return;
+  }
+  A_LEN=alen;
+  A_pFq = (double *)mymalloc(A_LEN*sizeof(double));
+  if (a != NULL) {
+	for (i=0; i<alen; i++) A_pFq[i] = a[i];
+  }else{
+	for (i=0; i<alen; i++) A_pFq[i] = 0.0;
+  }
+  return;
+}
+static void setB(double b[],int blen) {
+  int i;
+  if (blen < 0) {
+	if (B_pFq != NULL) myfree(B_pFq);
+	B_pFq=NULL; B_LEN=-1;
+	return;
+  }
+  if (blen == 0) {
+	B_LEN=0; return;
+  }
+  B_LEN=blen;
+  B_pFq = (double *)mymalloc(B_LEN*sizeof(double));
+  if (b != NULL) {
+	for (i=0; i<blen; i++) B_pFq[i] = b[i];
+  }else{
+	for (i=0; i<blen; i++) B_pFq[i] = 0.0;
+  }
+  return;
+}
+
 int jk_freeWorkArea() {
   /* bug, p in the cloneP will not be deallocated.
      Nk in genDarray2 will not be deallocated.
@@ -175,13 +311,14 @@ int jk_freeWorkArea() {
   if (M_beta_1) {myfree(M_beta_1); M_beta_1=NULL;}
   if (M_jack) {
     for (i=0; M_jack[i] != NULL; i++) {
-      if (Debug) printf("Free M_jack[%d]\n",i);
+      if (Debug) oxprintf("Free M_jack[%d]\n",i);
       myfree(M_jack[i]); M_jack[i] = NULL;
     }
     myfree(M_jack); M_jack=NULL;
   }
   if (M_qk) {myfree(M_qk); M_qk=NULL;}
   if (P_pki) {myfree(P_pki); P_pki=NULL;}
+  setA(NULL,-1); setB(NULL,-1);
   JK_deallocate=0;
   return(0);
 }
@@ -219,24 +356,24 @@ int jk_initializeWorkArea() {
 
 static void *mymalloc(int size) {
   void *p;
-  if (Debug) printf("mymalloc(%d)\n",size);
+  if (Debug) oxprintf("mymalloc(%d)\n",size);
   p = (void *)mh_malloc(size);
   if (p == NULL) {
-    fprintf(stderr,"No more memory.\n");
+    oxprintfe("No more memory.\n");
     mh_exit(-1);
   }
   return(p);
 }
 static int myfree(void *p) {
-  if (Debug) printf("myFree at %p\n",p);
+  if (Debug) oxprintf("myFree at %p\n",p);
   return(mh_free(p));
 }
-static int myerror(char *s) { fprintf(stderr,"%s: type in control-C\n",s); getchar(); getchar(); return(0);}
+static int myerror(char *s) { oxprintfe("%s: type in control-C\n",s); getchar(); getchar(); return(0);}
 
 static double jack1(int K) {
   double F;
   extern int Alpha;
-  int I,J,L,II,JJ,N;
+  int J,II,JJ,N;   /*  int I,J,L,II,JJ,N; */
   N = 1;
   if (K == 0) return((double)1);
   F = xval(1,K);
@@ -249,7 +386,7 @@ static double jack1(int K) {
 static double jack1diff(int K) {
   double F;
   extern int Alpha;
-  int I,J,S,L,II,JJ,N;
+  int J,II,JJ,N;  /* int I,J,S,L,II,JJ,N; */
   N = 1;
   if (K == 0) return((double)1);
   F = K*xval(1,K-1);
@@ -262,7 +399,6 @@ static double jack1diff(int K) {
 
 static double xval(int ii,int p) { /* x_i^p */
   extern double M_x[];
-  double F;
   int i,j;
   static int init=0;
   if (JK_deallocate) { init=0; return(0.0);}
@@ -282,7 +418,7 @@ static double xval(int ii,int p) { /* x_i^p */
   if (p > M_m_MAX-2) myerror("xval, p is too large.");
   if (p < 0) {
     myerror("xval, p is negative.");
-    printf("ii=%d, p=%d\n",ii,p);
+    oxprintf("ii=%d, p=%d\n",ii,p);
     mh_exit(-1);
   }
   return(Xarray[ii-1][p]);
@@ -331,6 +467,7 @@ static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m]
   }
 }
 
+#ifdef STANDALONE
 static int test_ptrans() {
   extern int M_m;
   int p[M_n0]={5,3,2};
@@ -338,11 +475,11 @@ static int test_ptrans() {
   int i;
   M_m = 10;
   ptrans(p,pt);
-  if (Debug) {for (i=0; i<10; i++) printf("%d,",pt[i]);  printf("\n");}
+  if (Debug) {for (i=0; i<10; i++) oxprintf("%d,",pt[i]);  oxprintf("\n");}
   return(0);
 }
+#endif
 
-
 /*
   upper hook length
   h_kappa^*(K)
@@ -352,7 +489,7 @@ static int huk(int K[],int I,int J) {
   int Kp[M_m_MAX];
   int A,H;
   A=Alpha;
-  /*printf("h^k(%a,%a,%a,%a)\n",K,I,J,A);*/
+  /*oxprintf("h^k(%a,%a,%a,%a)\n",K,I,J,A);*/
   ptrans(K,Kp);
   H=Kp[J-1]-I+A*(K[I-1]-J+1);
   return(H);
@@ -367,7 +504,7 @@ static int hdk(int K[],int I,int J) {
   int Kp[M_m_MAX];
   int A,H;
   A = Alpha;
-  /*printf("h_k(%a,%a,%a,%a)\n",K,I,J,A);*/
+  /*oxprintf("h_k(%a,%a,%a,%a)\n",K,I,J,A);*/
   ptrans(K,Kp);
   H=Kp[J-1]-I+1+A*(K[I-1]-J);
   return(H);
@@ -377,9 +514,9 @@ static int hdk(int K[],int I,int J) {
 */
 static double jjk(int K[]) {
   extern int Alpha;
-  int A,L,I,J;
+  int L,I,J;
   double V;
-  A=Alpha;
+
   V=1;
   L=plength(K);
   for (I=0; I<L; I++) {
@@ -387,7 +524,7 @@ static double jjk(int K[]) {
       V *= huk(K,I+1,J+1)*hdk(K,I+1,J+1);
     }
   }
-  if (Debug) {printp(K); printf("<--K, jjk=%lg\n",V);}
+  if (Debug) {printp(K); oxprintf("<--K, jjk=%lg\n",V);}
   return(V);
 }
 /*
@@ -433,7 +570,7 @@ static double mypower(double x,int n) {
 }
 /* Q_kappa
  */
-static double qk(int K[],double A[A_LEN],double B[B_LEN]) {
+static double qk(int K[],double A[],double B[]) {
   extern int Alpha;
   int P,Q,I;
   double V;
@@ -443,11 +580,21 @@ static double qk(int K[],double A[A_LEN],double B[B_LE
   /* to reduce numerical errors, temporary. */
   if (P == Q) {
     for (I=0; I<P; I++) V = V*ppoch2(A[I],B[I],K);
+	return(V);
   }
-  return(V);
 
-  for (I=0; I<P; I++) V = V*ppoch(A[I],K);
-  for (I=0; I<Q; I++) V = V/ppoch(B[I],K);
+  if (P > Q) {
+	for (I=0; I<Q; I++) V = V*ppoch2(A[I],B[I],K);
+	for (I=Q; I<P; I++) V = V*ppoch(A[I],K);
+  }else {
+	for (I=0; I<P; I++) V = V*ppoch2(A[I],B[I],K);
+	for (I=P; I<Q; I++) V = V/ppoch(B[I],K);
+  }    
+  /* for debug
+  printf("K="); 
+  for (I=0; I<3; I++) printf("%d, ",K[I]);
+  printf("qk=%lg\n",V); 
+  */
   return(V);
 }
 
@@ -461,8 +608,8 @@ static int bb(int N[],int K[],int M[],int I,int J) {
   ptrans(M,Mp); 
 
   /*
-    printp(K); printf("K<--, "); printp2(Kp); printf("<--Kp\n");
-    printp(M); printf("M<--, "); printp2(Mp); printf("<--Mp\n");
+    printp(K); oxprintf("K<--, "); printp2(Kp); oxprintf("<--Kp\n");
+    printp(M); oxprintf("M<--, "); printp2(Mp); oxprintf("<--Mp\n");
   */
 
   if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J));
@@ -483,7 +630,7 @@ static double beta(int K[],int M[]) {
     for (J=0; J<K[I]; J++) {
       II = I+1; JJ = J+1;
       V *= (double)bb(K,K,M,II,JJ);
-      /* printf("[%d,%d,%lf]\n",I,J,V); */
+      /* oxprintf("[%d,%d,%lf]\n",I,J,V); */
     }
   }
 
@@ -492,7 +639,7 @@ static double beta(int K[],int M[]) {
     for (J=0; J<M[I]; J++) {
       II = I+1; JJ = J+1;
       V /= (double)bb(M,K,M,II,JJ);
-      /* printf("[%d,%d,%lf]\n",I,J,V);*/
+      /* oxprintf("[%d,%d,%lf]\n",I,J,V);*/
     }
   }
 
@@ -500,20 +647,21 @@ static double beta(int K[],int M[]) {
 }
 static int printp(int kappa[]) {
   int i;
-  printf("(");
+  oxprintf("(");
   for (i=0; i<M_n; i++) {
-    if (i <M_n-1) printf("%d,",kappa[i]);
-    else printf("%d)",kappa[i]);
+    if (i <M_n-1) oxprintf("%d,",kappa[i]);
+    else oxprintf("%d)",kappa[i]);
   }
   return(0);
 }
+#ifdef STANDALONE
 static int printp2(int kappa[]) {
   int i,ell;
-  printf("(");
+  oxprintf("(");
   ell = plength_t(kappa);
   for (i=0; i<ell+1; i++) {
-    if (i <ell+1-1) printf("%d,",kappa[i]);
-    else printf("%d)",kappa[i]);
+    if (i <ell+1-1) oxprintf("%d,",kappa[i]);
+    else oxprintf("%d)",kappa[i]);
   }
   return(0);
 }
@@ -523,11 +671,11 @@ static int test_beta() {
   int mu1[M_n0]={1,0,0};
   int mu2[M_n0]={1,1,0};
   int mu3[M_n0]={2,0,0};
-  printp(kappa); printf(","); printp(mu3); printf(": beta = %lf\n",beta(kappa,mu3));
-  printp(kappa); printf(","); printp(mu1); printf(": beta = %lf\n",beta(kappa,mu1));
-  printp(kappa); printf(","); printp(mu2); printf(": beta = %lf\n",beta(kappa,mu2));
+  printp(kappa); oxprintf(","); printp(mu3); oxprintf(": beta = %lf\n",beta(kappa,mu3));
+  printp(kappa); oxprintf(","); printp(mu1); oxprintf(": beta = %lf\n",beta(kappa,mu1));
+  printp(kappa); oxprintf(","); printp(mu2); oxprintf(": beta = %lf\n",beta(kappa,mu2)); return(0);
 }
-
+#endif
 /* main() { test_beta(); } */
 
 
@@ -553,12 +701,12 @@ static int test_beta() {
 static double q3_10(int K[],int M[],int SK) {
   extern int Alpha;
   int Mp[M_m_MAX];
-  int ML[M_nmx];
+  //  int ML[M_nmx];
   int N[M_nmx];
   int i,R;
   double T,Q,V,Ur,Vr,Wr;
   ptrans(M,Mp);
-  for (i=0; i<M_n; i++) {ML[i] = M[i]; N[i] = M[i];}
+  for (i=0; i<M_n; i++) {N[i] = M[i];}
   N[SK-1] = N[SK-1]-1; 
 
   T = SK-Alpha*M[SK-1];
@@ -579,6 +727,7 @@ static double q3_10(int K[],int M[],int SK) {
   return(V);
 }
 
+#ifdef STANDALONE
 static double q3_5(double A[],double B[],int K[],int I) {
   extern int Alpha;
   int Kp[M_m_MAX];
@@ -611,17 +760,19 @@ static double q3_5(double A[],double B[],int K[],int I
   }
   return(V);
 }
-
+#endif
+#ifdef STANDALONE
 static void mtest4() {
-  double A[A_LEN] = {1.5};
-  double B[B_LEN]={6.5};
+  double A[1] = {1.5};
+  double B[1]={6.5};
   int K[M_n0] = {3,2,0};
   int I=2; 
   int Ki[M_n0]={3,1,0};
   double V1,V2;
+  setA(A,1); setB(B,1);
   V1=q3_5(A,B,K,I);
   V2=qk(K,A,B)/qk(Ki,A,B);
-  printf("%lf== %lf?\n",V1,V2);
+  oxprintf("%lf== %lf?\n",V1,V2);
 }
 static void mtest4b() {
   int K[M_n0]={3,2,0}; 
@@ -631,8 +782,9 @@ static void mtest4b() {
   double V1,V2;
   V1=q3_10(K,M,SK);
   V2=beta(K,N)/beta(K,M);
-  printf("%lf== %lf?\n",V1,V2);
+  oxprintf("%lf== %lf?\n",V1,V2);
 }
+#endif
 
 /* main() { mtest4(); mtest4b(); } */
 
@@ -652,6 +804,7 @@ static int nk(int KK[]) {
   /* K = (Kpp,Ki) */
   return(Darray[nk(Kpp)]+Ki-1);
 }
+#ifdef STANDALONE
 static int plength2(int P1[],int P2[]) {
   int S1,S2;
   S1 = plength(P1); S2 = plength(P2);
@@ -664,6 +817,7 @@ static int plength2(int P1[],int P2[]) {
   }
   else return(-1);
 }
+#endif
 static int myeq(int P1[],int P2[]) {
   int I,L1;
   if ((L1=plength(P1)) != plength(P2)) return(0);
@@ -691,7 +845,7 @@ static int pListPartition(int M,int N) {
   int I;
   /* initialize */
   if (M_n != N) {
-    fprintf(stderr,"M_n != N\n"); mh_exit(-1);
+    oxprintfe("M_n != N\n"); mh_exit(-1);
   }
   M_m = M;
   /* M_plist = []; */
@@ -709,26 +863,26 @@ static int pListPartition(int M,int N) {
   Less >= M_kap[From], ..., M_kap[To],  |(M_kap[From],...,M_kap[To])|<=M, 
 */
 static int pListPartition2(int Less,int From,int To, int M) {
-  int I,R;
+  int I;
   mh_check_intr(100);
   if (To < From) {
     (*M_pExec)(); return(0);
   }
   for (I=1; (I<=Less) && (I<=M) ; I++) {
     M_kap[From-1] = I;
-    R = pListPartition2(I,From+1,To,M-I);
+    pListPartition2(I,From+1,To,M-I);
   } 
   return(1);
 }
 
 /*
-  partition に対してやる仕事をこちらへ書く.
+  Commands to do for each partition are given here.
 */
 static void pExec_0() {
   if (Debug) {
-    printf("M_kap=");
+    oxprintf("M_kap=");
     printp(M_kap);
-    printf("\n");
+    oxprintf("\n");
   }
 }
 
@@ -772,8 +926,8 @@ static int pListHS2(int From,int To,int Kap[]) {
 }
 
 static void hsExec_0() {
-  int i;
-  if(Debug) {printf("hsExec: "); printp(HS_mu); printf("\n");}
+  /* int i; */
+  if(Debug) {oxprintf("hsExec: "); printp(HS_mu); oxprintf("\n");}
 }
 
 /*
@@ -811,10 +965,10 @@ static int pmn(int M,int N) {
   }
   P_pmn=S;
   if (Debug) {
-    printf("P_pmn=%d\n",P_pmn);
+    oxprintf("P_pmn=%d\n",P_pmn);
     for (i=0; i<=Min_m_n; i++) {
-      for (j=0; j<=M; j++) printf("%d,",aP_pki(i,j));
-      printf("\n");
+      for (j=0; j<=M; j++) oxprintf("%d,",aP_pki(i,j));
+      oxprintf("\n");
     }
   }
   myfree(P_pki); P_pki=NULL;
@@ -822,7 +976,7 @@ static int pmn(int M,int N) {
 }
 
 /*
-  main() {pmn(4,3); printf("P_pmn=%d\n",P_pmn);}
+  main() {pmn(4,3); oxprintf("P_pmn=%d\n",P_pmn);}
 */
 
 static int *cloneP(int a[]) {
@@ -863,7 +1017,7 @@ static int genDarray2(int M,int N) {
 
   M_m = M;
   Pmn = pmn(M,N)+1;
-  if (Debug) printf("Degree M = %d, N of vars N = %d, Pmn+1=%d\n",M,N,Pmn);
+  if (Debug) oxprintf("Degree M = %d, N of vars N = %d, Pmn+1=%d\n",M,N,Pmn);
   Darray=(int *) mymalloc(sizeof(int)*Pmn);
   for (i=0; i<Pmn; i++) Darray[i] = 0;
   Parray=(int **) mymalloc(sizeof(int *)*Pmn);
@@ -882,7 +1036,7 @@ static int genDarray2(int M,int N) {
     K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */
     Ksize = plength(K);
     if (Ksize >= M_n) {
-      if (Debug) {fprintf(stderr,"Ksize >= M_n\n");}
+      if (Debug) {oxprintfe("Ksize >= M_n\n");}
       continue;
     }
     for (i=0; i<M_n; i++) Kone[i] = 0;
@@ -892,9 +1046,9 @@ static int genDarray2(int M,int N) {
     }
   }
   if (Debug) {
-    printf("Darray=\n");
-    for (i=0; i<Pmn; i++) printf("%d\n",Darray[i]);
-    printf("-----------\n");
+    oxprintf("Darray=\n");
+    for (i=0; i<Pmn; i++) oxprintf("%d\n",Darray[i]);
+    oxprintf("-----------\n");
   }
   return(0);
 }
@@ -923,7 +1077,7 @@ static void hsExec_beta(void) {
   int Nu[M_nmx];
   int rrMax;
   hsExec_0();
-  /* printf("M_beta_pt=%a\n",M_beta_pt); */
+  /* oxprintf("M_beta_pt=%a\n",M_beta_pt); */
   /* Mu = cdr(vtol(HS_mu)); */
   Mu = HS_mu; /* buggy? need cloneP */
   if (M_beta_pt == 0) {
@@ -947,22 +1101,22 @@ static void hsExec_beta(void) {
       if (M_beta_1[J] == Nnu) {
         K=I+1;
         if (Debug) {
-          printf("Found at J=%d, K=%d, q3_10(Kap,Nu,K)=%lf,Nu,Mu= \n",
+          oxprintf("Found at J=%d, K=%d, q3_10(Kap,Nu,K)=%lf,Nu,Mu= \n",
                  J,K,q3_10(M_beta_kap,Nu,K));
-          printp(Nu); printf("\n");
-          printp(Mu); printf("\n");
+          printp(Nu); oxprintf("\n");
+          printp(Mu); oxprintf("\n");
         }
         /* Check other conditions. See Numata's mail on Dec 24, 2011. */
         rrMax = Nu[I]-1;                 
         if ((plength_t(Kapt) < rrMax) || (plength_t(Nut) < rrMax)) {
-          if (Debug) printf(" is not taken (length). \n");
+          if (Debug) oxprintf(" is not taken (length). \n");
           break;
         }
         OK=1; 
         for (RR=0; RR<rrMax; RR++) {
           if (Kapt[RR] != Nut[RR]) { OK=0; break;}
         }
-        if (!OK) { if (Debug) printf(" is not taken.\n"); break; }
+        if (!OK) { if (Debug) oxprintf(" is not taken.\n"); break; }
         /* check done. */
         M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K);
         Done = 1; break;
@@ -971,7 +1125,7 @@ static void hsExec_beta(void) {
     if (Done) break; else Nu[I]--;
   }
   if (!Done) {
-    if (Debug) printf("BUG: not found M_beta_pt=%d.\n",M_beta_pt);
+    if (Debug) oxprintf("BUG: not found M_beta_pt=%d.\n",M_beta_pt);
     /* M_beta_0[M_beta_pt] = NAN;  error("Not found."); */
     M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu); 
   }
@@ -984,8 +1138,8 @@ static int genBeta(int Kap[]) {
   extern int M_beta_pt;
   extern int M_beta_kap[];
   extern int P_pmn;
-  int I,J,N;
-  if (Debug) {printp(Kap); printf("<-Kappa, P_pmn=%d\n",P_pmn);}
+  int I,N;
+  if (Debug) {printp(Kap); oxprintf("<-Kappa, P_pmn=%d\n",P_pmn);}
   /* M_beta = newmat(2,P_pmn+1); */
   M_beta_0 = (double *)mymalloc(sizeof(double)*(P_pmn+1));
   M_beta_1 = (int *)mymalloc(sizeof(int)*(P_pmn+1));
@@ -994,14 +1148,14 @@ static int genBeta(int Kap[]) {
   N = plength(Kap);
   HS_hsExec = hsExec_beta;
   copyP(M_beta_kap,Kap);
-  pListHS(Kap,N);
+  pListHS(Kap,N); return(0);
 }
 /*
   genDarray2(4,3);
   genBeta([2,2,0]);
   genBeta([2,1,1]);
 */
-
+#ifdef STANDALONE
 static int checkBeta1() {
   int Kap[3] = {2,2,0};
   int Kap2[3] = {2,1,0};
@@ -1014,36 +1168,36 @@ static int checkBeta1() {
     Mu = Parray[M_beta_1[I]];
     Beta_km = M_beta_0[I];
     if (Debug) {
-      printp(Kap); printf("<--Kap, ");
-      printp(Mu); printf("<--Mu,");
-      printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap,Mu));
+      printp(Kap); oxprintf("<--Kap, ");
+      printp(Mu); oxprintf("<--Mu,");
+      oxprintf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap,Mu));
     }
   }
-  if (Debug) printf("-------------------------------------\n");
+  if (Debug) oxprintf("-------------------------------------\n");
   genBeta(Kap2);
   for (I=0; I<M_beta_pt; I++) {
     Mu = Parray[M_beta_1[I]];
     Beta_km = M_beta_0[I];
     if (Debug) {
-      printp(Kap2); printf("<--Kap, ");
-      printp(Mu); printf("<--Mu,");
-      printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap2,Mu));
+      printp(Kap2); oxprintf("<--Kap, ");
+      printp(Mu); oxprintf("<--Mu,");
+      oxprintf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap2,Mu));
     }
   }
   return(0);
 }
-
+#endif
 /*
   def checkBeta2() {
   genDarray2(3,3);
   Kap = [2,1,0];
-  printf("Kap=%a\n",Kap);
+  oxprintf("Kap=%a\n",Kap);
   genBeta(Kap);
   for (I=0; I<M_beta_pt; I++) {
   Mu = Parray[M_beta[1][I]];
   Beta_km = M_beta[0][I];
-  printf("Mu=%a,",Mu);
-  printf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));
+  oxprintf("Mu=%a,",Mu);
+  oxprintf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));
   }
   }
 */
@@ -1086,8 +1240,8 @@ static int genJack(int M,int N) {
   int Pmn,I,J,K,L,Nv,H,P;
   int *Kap,*Mu;
   double Jack,Beta_km;
-  int Nk,JJ;
-  if (Debug) printf("genJack(%d,%d)\n",M,N);
+  int Nk,JJ, two_to_I;
+  if (Debug) oxprintf("genJack(%d,%d)\n",M,N);
   M_jack = (double **) mymalloc(sizeof(double *)*(N+2));
   M_2n = imypower(2,N);
   Pmn = pmn(M,N);  /*P_pmn is initializeded. 
@@ -1110,12 +1264,12 @@ static int genJack(int M,int N) {
       for (J=2; J<M_2n; J++) aM_jack(1,J,K) = 0;
     }
   }
-  for (I=1; I<=N; I++) {
+  for (I=1; I<=N; I++) {   two_to_I = imypower(2,I);
     for (K=M+1; K<Pmn+1; K++) {
       aM_jack(I,0,K) = NAN;
       if (M_df) {
         for (J=1; J<M_2n; J++) {
-          if (J >= 2^I) aM_jack(I,J,K) = 0;
+          if (J >= two_to_I) aM_jack(I,J,K) = 0; /* J >= 2^I */
           else aM_jack(I,J,K) = NAN;
         }
       }
@@ -1132,7 +1286,7 @@ static int genJack(int M,int N) {
         for (J=1; J<M_2n; J++) aM_jack(I,J,K) = 0;
       } 
     }
-    if (Debug) {printf("Kappa="); printp(Kap);}
+    if (Debug) {oxprintf("Kappa="); printp(Kap);}
     /* Enumerate horizontal strip of Kappa */
     genBeta(Kap);  /* M_beta_pt stores the number of hs */
     /* Nv is the number of variables */
@@ -1147,11 +1301,11 @@ static int genJack(int M,int N) {
           Beta_km = beta(Kap,Mu);
           /* do not use the M_beta table. It's buggy. UseTable is experimental.*/
         }
-        if (Debug) {printf("Nv(number of variables)=%d, Beta_km=%lf, Mu=",Nv,Beta_km);
-          printp(Mu); printf("\n");}
+        if (Debug) {oxprintf("Nv(number of variables)=%d, Beta_km=%lf, Mu=",Nv,Beta_km);
+          printp(Mu); oxprintf("\n");}
         P = psublen(Kap,Mu);
         Jack += aM_jack(Nv-1,0,Nk)*Beta_km*xval(Nv,P); /* util_v(x,[Nv])^P;*/
-        if (Debug) printf("xval(%d,%d)=%lf\n",Nv,P,xval(Nv,P));
+        if (Debug) oxprintf("xval(%d,%d)=%lf\n",Nv,P,xval(Nv,P));
       }
       aM_jack(Nv,0,K) = Jack;
       if (M_df) {
@@ -1167,8 +1321,8 @@ static int genJack(int M,int N) {
             }else{
               Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */
             }
-            if (Debug) {printf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);
-              printp(Mu); printf("\n"); }
+            if (Debug) {oxprintf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);
+              printp(Mu); oxprintf("\n"); }
             P = psublen(Kap,Mu);
             if (J & (1 << (Nv-1))) {
               JJ = J & ((1 << (Nv-1)) ^ 0xffff);  /* NOTE!! Up to 16 bits. mh-15 */
@@ -1180,14 +1334,14 @@ static int genJack(int M,int N) {
             }
           }
           aM_jack(Nv,J,K) = Jack;
-          if (Debug) printf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);
+          if (Debug) oxprintf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);
         } /* end of J loop */
       }
     }
-  }
+  } return(0);
 }  
 
-
+#ifdef STANDALONE
 /* checkJack1(3,3) 
  */
 static int checkJack1(int M,int N) {
@@ -1204,11 +1358,11 @@ static int checkJack1(int M,int N) {
   for (I=1; I<=N; I++) {
     for (K=0; K<=P_pmn; K++) {
       printp(Parray[K]);
-      printf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
+      oxprintf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
     }
   }
-  for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);
-  printf("<--x\n");
+  for (I=1; I<=N; I++) oxprintf("%lf, ",M_x[I-1]);
+  oxprintf("<--x\n");
   return(0);
 }
 /*main() {  checkJack1(3,3);  }*/
@@ -1230,17 +1384,17 @@ static int checkJack2(int M,int N) {
   for (I=1; I<=N; I++) {
     for (K=0; K<=P_pmn; K++) {
       printp(Parray[K]);
-      printf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
+      oxprintf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
     }
   }
-  for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);
-  printf("<--x\n");
+  for (I=1; I<=N; I++) oxprintf("%lf, ",M_x[I-1]);
+  oxprintf("<--x\n");
 
   for (I=1; I<=N; I++) {
     for (K=0; K<=P_pmn; K++) {
       for (J=0; J<M_2n; J++) {
         printp(Parray[K]);
-        printf("<--Kap, Nv=%d,J(diff)=%d, D^J Jack=%lf\n",
+        oxprintf("<--Kap, Nv=%d,J(diff)=%d, D^J Jack=%lf\n",
                I,J,aM_jack(I,J,K));
       }
     }
@@ -1249,17 +1403,25 @@ static int checkJack2(int M,int N) {
 }
 
 /* main() { checkJack2(3,3); } */
+#endif
 
-double mh_t(double A[A_LEN],double B[B_LEN],int N,int M) {
+double mh_t(double A[],double B[],int N,int M) {
   double F,F2;
   extern int M_df;
   extern int P_pmn;
   extern double *M_qk;
   extern double M_rel_error;
+  extern int M_m;
+  extern int M_m_estimated_approx_deg;
+  extern double M_assigned_series_error;
   int Pmn;
   int K;
   int *Kap;
   int size;
+  int i;
+  double partial_sum[M_m_MAX+1];
+  double iv;
+  double serror;
   F = 0; F2=0;
   M_df=1;
   genJack(M,N);
@@ -1272,10 +1434,57 @@ double mh_t(double A[A_LEN],double B[B_LEN],int N,int 
     M_qk[K] = qk(Kap,A,B);
     F += M_qk[K]*aM_jack(N,0,K);
     if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K);
-    if (Debug) printf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size);
-    if (Debug && (ParraySize[K] == size)) printf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K));
+    if (Debug) oxprintf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size);
+    if (Debug && (ParraySize[K] == size)) oxprintf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K));
   }
   M_rel_error = F-F2;
+
+  M_m_estimated_approx_deg = -1; serror=1;
+  for (i=0; i<=M_m; i++) {
+    partial_sum[i] = 0.0; partial_sum[i+1] = 0.0;
+    for (K=0; K<=P_pmn; K++) {
+      if (ParraySize[K] == i) partial_sum[i] += M_qk[K]*aM_jack(N,0,K);
+    }
+    if (i>0) partial_sum[i] += partial_sum[i-1];
+    if (i>0) serror = myabs((partial_sum[i]-partial_sum[i-1])/partial_sum[i-1]);
+    if ((i>0)&&(M_m_estimated_approx_deg < 0)&&(serror<M_assigned_series_error)) {
+      M_m_estimated_approx_deg = i; break;
+    }
+  }
+  if (M_m_estimated_approx_deg < 0) {
+    M_m_estimated_approx_deg = M_m+mymin(5,mymax(1,(int)log(serror/M_assigned_series_error))); /* Heuristic */
+  }
+  /*
+  for (K=0; K<=P_pmn; K++) {
+    oxprintf("Kappa="); for (i=0; i<N; i++) oxprintf("%d ",Parray[K][i]); oxprintf("\n"); 
+    oxprintf("ParraySize(%d)=%d (|kappa|),   M_m=%d\n",K,ParraySize[K],M_m);
+  }
+  for (i=0; i<=M_m; i++) {
+    oxprintf("partial_sum[%d]=%lg\n",i,partial_sum[i]);
+  }
+  */
+  M_estimated_X0g = X0g;
+  iv=myabs(F*iv_factor());
+  if (iv < M_x0value_min) M_estimated_X0g = X0g*mymax(2,log(log(1/iv)));   /* This is heuristic */
+  M_estimated_X0g = mymin(M_estimated_X0g,M_X0g_bound);
+  M_mh_t_success = 1;
+  if (M_estimated_X0g != X0g) M_mh_t_success=0;
+  if (M_m_estimated_approx_deg > M_m) M_mh_t_success=0;
+
+  M_series_error = serror;
+  M_recommended_abserr = iv*M_assigned_series_error;
+  M_recommended_relerr = M_series_error;
+
+  if (M_show_autosteps) {
+    oxprintf("%%%%serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m);
+    oxprintf("%%%%x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
+    oxprintf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
+    oxprintfe("%%%%(stderr) serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m);
+    oxprintfe("%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
+    oxprintfe("%%%%(stderr) F=%lg,Ef=%lg,iv=|F*Ef|=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),iv,M_estimated_X0g,X0g);
+  }
+
+  M_mh_t_value=F;
   return(F);
 }
 double mh_t2(int J) {
@@ -1293,6 +1502,7 @@ double mh_t2(int J) {
   return(F);
 }
 
+#ifdef STANDALONE
 static int mtest1b() {
   double A[1] = {1.5};
   double B[1] = {1.5+5};
@@ -1305,15 +1515,16 @@ static int mtest1b() {
   mh_t(A,B,N,M);
   for (J=0; J<M_2n; J++) {
     F=mh_t2(J);
-    printf("J=%d, D^J mh_t=%lf\n",J,F);
+    oxprintf("J=%d, D^J mh_t=%lf\n",J,F);
   }
+  return(0);
 }
 
 /* main() { mtest1b(); }*/
+#endif
 
 
 
-
 #define TEST 1
 #ifndef TEST
 
@@ -1322,7 +1533,7 @@ static int mtest1b() {
 /****** from mh-n.c *****/ 
 
 #define SMAX 4096
-#define inci(i) { i++; if (i >= argc) { fprintf(stderr,"Option argument is not given.\n"); return(NULL); }}
+#define inci(i) { i++; if (i >= argc) { oxprintfe("Option argument is not given.\n"); return(NULL); }}
 
 static int imypower(int x,int n) {
   int i;
@@ -1334,23 +1545,48 @@ static int imypower(int x,int n) {
 }
 
 #ifdef STANDALONE2
-main(int argc,char *argv[]) {
+int main(int argc,char *argv[]) {
   mh_exit(MH_RESET_EXIT);
   /*  jk_main(argc,argv);
-      printf("second run.\n"); */
+      oxprintf("second run.\n"); */
   jk_main(argc,argv);
+  return(0);
 }
 #endif
 
 struct MH_RESULT *jk_main(int argc,char *argv[]) {
-  double *y0;
-  double x0,xn;
-  double ef;
-  int i,j,rank;
-  double a[1]; double b[1];
+  int i;
+  struct MH_RESULT *ans;
+  extern int M_automatic;
+  extern int M_mh_t_success;
+  extern double M_estimated_X0g;
+  extern int M_m_estimated_approx_deg;
+  for (i=1; i<argc; i++) {
+    if (strcmp(argv[i],"--automatic")==0) {
+      inci(i);
+      sscanf(argv[i],"%d",&M_automatic);
+      break;
+    }
+  }
+  ans=jk_main2(argc,argv,0,0.0,0);
+  if (!M_automatic) return(ans);
+  if (M_mh_t_success) return(ans);
+  while (!M_mh_t_success) {
+    ans=jk_main2(argc,argv,1,M_estimated_X0g,M_m_estimated_approx_deg);
+  }
+  return(ans);
+}
+
+struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree) {
+  // double *y0;
+  //  double x0,xn;  
+  // double ef;
+
+  int i,j; // int i,j,rank;
   extern double M_x[];
   extern double *Beta;
   extern int M_2n;
+  extern int M_mh_t_success;
   char swork[1024];
   struct MH_RESULT *ans=NULL;
   struct SFILE *ofp = NULL;
@@ -1376,21 +1612,43 @@ struct MH_RESULT *jk_main(int argc,char *argv[]) {
     }else if (strcmp(argv[i],"--help")==0) {
       usage(); return(0);
     }else if (strcmp(argv[i],"--bystring")==0) {
-      if (idata) {fprintf(stderr,"--bystring must come before --idata option.\n"); mh_exit(-1);}
+      if (idata) {oxprintfe("--bystring must come before --idata option.\n"); mh_exit(-1);}
       JK_byFile = 0;
+    }else if (strcmp(argv[i],"--automatic")==0) {
+      inci(i); /* ignore, in this function */
+    }else if (strcmp(argv[i],"--assigned_series_error")==0) {
+      inci(i);
+      sscanf(argv[i],"%lg",&M_assigned_series_error);
+    }else if (strcmp(argv[i],"--x0value_min")==0) {
+      inci(i);
+      sscanf(argv[i],"%lg",&M_x0value_min);
     }else {
-      fprintf(stderr,"Unknown option %s\n",argv[i]);
+      oxprintfe("Unknown option %s\n",argv[i]);
       usage();
       return(NULL);
     }
   }
-  if (!idata) setParam(NULL);  
+  if (!idata) setParam(NULL);
+  if (automode) {
+    Mapprox = newDegree;
+    X0g = newX0g;
+  }
 
   /* Initialize global variables */
   M_n = Mg;
   HS_n=M_n;
-  if (!JK_byFile) ans = (struct MH_RESULT *)mymalloc(sizeof(struct MH_RESULT));
+  if (!JK_byFile) {
+    ans = (struct MH_RESULT *)mymalloc(sizeof(struct MH_RESULT));
+    ans->message = NULL;
+    ans->t_success = 0;
+    ans->series_error = 1.0e+10;
+    ans->recommended_abserr = 1.0e-10;
+  }
   else ans = NULL;
+  if (M_automatic) {
+    /* Differentiation can be M_m in the bit pattern in the M_n variable case.*/
+    if (M_n > Mapprox) Mapprox=M_n;
+  }
   /* Output by a file=stdout */
   ofp = mh_fopen("stdout","w",JK_byFile);
 
@@ -1400,13 +1658,24 @@ struct MH_RESULT *jk_main(int argc,char *argv[]) {
     myerror("Mg must be equal to M_n\n"); mh_exit(-1);
   }
   if (Debug) showParam(NULL,1);
+  setM_x();
+
+  M_beta_i_x_o2_max=myabs(M_x[0]/2);
+  if (M_n <= 1) M_beta_i_beta_j_min = myabs(Beta[0]);
+  else M_beta_i_beta_j_min = myabs(Beta[1]-Beta[0]);
   for (i=0; i<M_n; i++) {
-    M_x[i] = Beta[i]*X0g;
-  } 
-  a[0] = ((double)Mg+1.0)/2.0;
-  b[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */
-  if (Debug) printf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",a[0],b[0],M_n,Mapprox);
-  mh_t(a,b,M_n,Mapprox);
+    if (myabs(M_x[i]/2) > M_beta_i_x_o2_max) M_beta_i_x_o2_max = myabs(M_x[i]/2);
+    for (j=i+1; j<M_n; j++) {
+      if (myabs(Beta[i]-Beta[j]) < M_beta_i_beta_j_min) 
+        M_beta_i_beta_j_min = myabs(Beta[i]-Beta[j]);
+    }
+  }
+
+  mh_t(A_pFq,B_pFq,M_n,Mapprox);
+  if ((!M_mh_t_success) && M_automatic) {
+    jk_freeWorkArea();
+    return NULL;
+  }
   if (imypower(2,M_n) != M_2n) {
     sprintf(swork,"M_n=%d,M_2n=%d\n",M_n,M_2n); mh_fputs(swork,ofp);
     myerror("2^M_n != M_2n\n"); mh_exit(-1);
@@ -1419,7 +1688,7 @@ struct MH_RESULT *jk_main(int argc,char *argv[]) {
   showParam(ofp,0);
 
   /* return the result */
-  if (!JK_byFile) {
+  if (!JK_byFile) { 
     ans->x = X0g;
     ans->rank = imypower(2,Mg);
     ans->y = (double *)mymalloc(sizeof(double)*(ans->rank));
@@ -1427,40 +1696,20 @@ struct MH_RESULT *jk_main(int argc,char *argv[]) {
     ans->size=1;
     ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size));
     (ans->sfpp)[0] = ofp;
+
+    ans->t_success = M_mh_t_success;
+    ans->series_error = M_series_error;
+    ans->recommended_abserr = M_recommended_abserr;
   }
-  if (Debug) printf("jk_freeWorkArea() starts\n");
+  if (Debug) oxprintf("jk_freeWorkArea() starts\n");
   jk_freeWorkArea();
-  if (Debug) printf("jk_freeWorkArea() has finished.\n");
+  if (Debug) oxprintf("jk_freeWorkArea() has finished.\n");
   return(ans);
 }
 
 static int usage() {
-  fprintf(stderr,"Usages:\n");
-  fprintf(stderr,"hgm_jack-n [--idata input_data_file --x0 x0 --degree approxm]\n");
-  fprintf(stderr,"\nThe command hgm_jack-n [options] generates an input for hgm_w-n, Pr({y | y<xmax}), which is the cumulative distribution function of the largest root of the m by m Wishart matrix with n degrees of freedom and the covariantce matrix sigma.\n");
-  fprintf(stderr,"The hgm_jack-n uses the Koev-Edelman algorithm to evalute the matrix hypergeometric function.\n");
-  fprintf(stderr,"The degree of the approximation (Mapprox) is given by the --degree option.\n");
-  fprintf(stderr,"Parameters are specified by the input_data_file. Otherwise, default values are used.\n");
-  fprintf(stderr,"The format of the input_data_file. The orders of the input numbers must be kept.\n");
-  fprintf(stderr," Mg: m(the number of variables), Beta: beta=sigma^(-1)/2 (diagonized), Ng: n,\n");
-  fprintf(stderr," (Add a comment line %Ng= before the data Ng to check the number of beta.)\n");
-  fprintf(stderr," X0g: starting value of x(when --x0 option is used, this value is used)\n");
-  fprintf(stderr," Iv: initial values at X0g*Beta (see our paper how to order them), are evaluated in this program. Give zeros or the symbol * to skip rank many inputs.\n");
-  fprintf(stderr," Ef: a scalar factor to the initial value. It is calculated by this program. Give the zero.\n");
-  fprintf(stderr," Hg: h (step size) which is for hgm_w-n, \n"); 
-  fprintf(stderr," Dp: output data is stored in every Dp steps when output_data_file is specified, which is for hgm_w-n.\n");
-  fprintf(stderr," Xng: terminating value of x which is for hgm_w-n.\n");
-  fprintf(stderr," The line started with %% is a comment line.\n");
-  fprintf(stderr," With the --notable option, it does not use the Lemma 3.2 of Koev-Edelman (there is a typo: kappa'_r = mu'_r for 1<=r<=mu_k).\n");
-  fprintf(stderr," An example format of the input_data_file can be obtained by executing hgm_jack-n with no option.\n");
-  fprintf(stderr,"\nExamples:\n");
-  fprintf(stderr,"[1] ./hgm_jack-n \n");
-  fprintf(stderr,"[2] ./hgm_jack-n --x0 0.1 \n");
-  fprintf(stderr,"[3] ./hgm_jack-n --x0 0.1 --degree 15 \n");
-  fprintf(stderr,"[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --degree 15 \n");
-  fprintf(stderr,"[5] ./hgm_jack-n --degree 15 >test2.txt\n");
-  fprintf(stderr,"    ./hgm_w-n --idata test2.txt --gnuplotf test-g\n");
-  fprintf(stderr,"    gnuplot -persist <test-g-gp.txt\n");
+  oxprintfe("Usages:\n");
+#include "usage-jack-n.h"
   return(0);
 }
 
@@ -1487,24 +1736,33 @@ static int setParamDefault() {
   Hg = 0.001;
   Dp = 1;
   Xng = 10.0;
+  setA(NULL,1); setB(NULL,1);
+  A_pFq[0] = ((double)Mg+1.0)/2.0;
+  B_pFq[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; 
   return(0);
 }
 
 static int next(struct SFILE *sfp,char *s,char *msg) {
   static int check=1;
-  char *ng="%Ng=";
-  int i;
+  char *ng="%%Ng=";
+  // int i;
   s[0] = '%';
-  while (s[0] == '%') {
+  while ((s[0] == '%') || (s[0] == '#')) {
     if (!mh_fgets(s,SMAX,sfp)) {
-      fprintf(stderr,"Data format error at %s\n",msg);
+      oxprintfe("Data format error at %s\n",msg);
+      oxprintfe("Is it version 2.0 format? If so, add\n%s\nat the top.\n",VSTRING);
       mh_exit(-1);
     }
+	if ((s[0] == '%') && (s[1] == '%')) continue;
+    if (s[0] == '#') continue;
+    if (strncmp(s,VSTRING,strlen(VSTRING)) == 0) {
+	  return(2);
+	}
     if (check && (strncmp(msg,ng,4)==0)) {
-      if (strncmp(s,ng,4) != 0) {
-        fprintf(stderr,"Warning, there is no %%Ng= at the border of Beta's and Ng, s=%s\n",s);
+      if (strncmp(s,ng,5) != 0) {
+        oxprintfe("Warning, there is no %%Ng= at the border of Beta's and Ng, s=%s\n",s);
       }
-      check=0;
+      /* check=0; */
     }
     if (s[0] != '%') return(0);
   }
@@ -1515,14 +1773,28 @@ static int setParam(char *fname) {
   char s[SMAX];
   struct SFILE *fp;
   int i;
+  struct mh_token tk;
+  int version;
   if (fname == NULL) return(setParamDefault());
 
   Sample = 0;
   if ((fp=mh_fopen(fname,"r",JK_byFile)) == NULL) {
-    if (JK_byFile) fprintf(stderr,"File %s is not found.\n",fp->s);
+    if (JK_byFile) oxprintfe("File %s is not found.\n",fname);
     mh_exit(-1);
   }
-  next(fp,s,"Mg(m)");
+  /* set default initial values */
+  Mg=-1;  /* number of variables */
+  Ng=(double *) mymalloc(sizeof(double)); *Ng=-1; /* *Ng is the degree of freedom 1F1 */
+  X0g=0.1;   /* evaluation point */
+  Ef=1.0;    /* exponential factor */
+  Ef_type=1;
+  Hg=0.01;   /* step size for RK */
+  Dp = 1;   /* sampling rate */
+  Xng = 10.0; /* terminal point for RK */
+
+  /* Parser for the old style (version <2.0) input */
+  version=next(fp,s,"Mg(m)");
+  if (version == 2) goto myparse;
   sscanf(s,"%d",&Mg);
   rank = imypower(2,Mg);
 
@@ -1532,7 +1804,6 @@ static int setParam(char *fname) {
     sscanf(s,"%lf",&(Beta[i]));
   }
 
-  Ng = (double *)mymalloc(sizeof(double));
   next(fp,s,"%Ng= (freedom parameter n)");
   sscanf(s,"%lf",Ng);
   
@@ -1561,11 +1832,267 @@ static int setParam(char *fname) {
 
   next(fp,s,"Xng (the last point, cf. --xmax)");
   sscanf(s,"%lf",&Xng);
+
+  /* Reading the optional parameters */
+ myparse: 
+  while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) {
+    /* expect ID */
+    if (tk.type != MH_TOKEN_ID) {
+      oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+    }
+    if (strcmp(s,"automatic")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      M_automatic = tk.ival;
+      continue;
+    }
+    if (strcmp(s,"assigned_series_error")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      M_assigned_series_error = tk.dval;
+      continue;
+    }
+    if (strcmp(s,"x0value_min")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      M_x0value_min = tk.dval;
+      continue;
+    }
+    if ((strcmp(s,"Mapprox")==0) || (strcmp(s,"degree")==0)) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      Mapprox = tk.ival;
+      continue;
+    }
+    if (strcmp(s,"X0g_bound")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      M_X0g_bound = tk.dval;
+      continue;
+    }
+    if (strcmp(s,"show_autosteps")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      M_show_autosteps = tk.ival;
+      continue;
+    }
+    // Format: #p_pFq=2  1.5  3.2
+    if (strcmp(s,"p_pFq")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      setA(NULL,tk.ival);
+      for (i=0; i<A_LEN; i++) { 
+	if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+	  A_pFq[i] = tk.dval;
+	}else if (tk.type == MH_TOKEN_INT) {
+	  A_pFq[i] = tk.ival;
+	}else{
+	  oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+	}
+      }
+      continue;
+    }
+    if (strcmp(s,"q_pFq")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      setB(NULL,tk.ival);
+      for (i=0; i<B_LEN; i++) { 
+	if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+	  B_pFq[i] = tk.dval;
+	}else if (tk.type == MH_TOKEN_INT) {
+	  B_pFq[i] = tk.ival;
+	}else{
+	  oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+	}
+      }
+      continue;
+    }
+    if (strcmp(s,"ef_type")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      Ef_type = tk.ival;
+      continue;
+    }
+
+    if (strcmp(s,"Mg")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      Mg = tk.ival;
+      rank = imypower(2,Mg);
+      continue;
+    }
+    if (strcmp(s,"Beta")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+	  if (Mg <= 0) {
+        oxprintfe("Mg should be set before reading Beta.\n"); mh_exit(-1);
+      }
+      Beta = (double *)mymalloc(sizeof(double)*Mg);
+	  for (i=0; i<Mg; i++) {
+        if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+           Beta[i] = tk.dval;
+        }else if (tk.type == MH_TOKEN_INT) {
+           Beta[i] = tk.ival;
+        }else {
+          oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+        }
+      }
+      Iv = (double *)mymalloc(sizeof(double)*rank);
+	  for (i=0; i<rank; i++) {
+		Iv[i] = 0.0;
+	  }
+      continue;
+    }
+    if (strcmp(s,"Ng")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type== MH_TOKEN_DOUBLE) {
+		*Ng = tk.dval;
+      }else if (tk.type == MH_TOKEN_INT) {
+        *Ng = tk.ival;
+      }else{
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      continue;
+    }
+    if (strcmp(s,"X0g")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+		X0g = tk.dval;
+	  }else if (tk.type == MH_TOKEN_INT) {
+		X0g = tk.ival;
+      }else{
+		oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+	  }			
+      continue;
+    }
+	if (strcmp(s,"Iv")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+	  for (i=0; i<rank; i++) {
+        if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+		  Iv[i] = tk.dval;
+		}else if (tk.type == MH_TOKEN_INT) {
+		  Iv[i] = tk.ival;
+        }else{
+          oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+		}
+      }
+      continue;
+    }
+    if (strcmp(s,"Ef")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+		Ef = tk.dval;
+	  }else if (tk.type == MH_TOKEN_INT) {
+		Ef = tk.ival;
+      }else{
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+	  }
+      continue;
+    }
+    if (strcmp(s,"Hg")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+		Hg = tk.dval;
+	  }else if (tk.type == MH_TOKEN_INT) {
+		Hg = tk.ival;
+      }else{
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      continue;
+    }
+    if (strcmp(s,"Dp")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      Dp = tk.dval;
+      continue;
+    }
+    if (strcmp(s,"Xng")==0) {
+      if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
+		Xng = tk.dval;
+	  }else if (tk.type == MH_TOKEN_INT) {
+		Xng = tk.ival;
+      }else{
+        oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
+	  }
+      continue;
+    }
+
+    oxprintfe("Unknown ID at %s\n",s); mh_exit(-1);
+  }
+  
+  /* 1F1, original wishart case. */  
+  if ((A_LEN <= 1) && (B_LEN <= 1) && (Ef_type==1) && (*Ng >= 0)) {
+    if (A_LEN <1) setA(NULL,1);
+	if (B_LEN <1) setB(NULL,1);
+    A_pFq[0] = ((double)Mg+1.0)/2.0;
+    B_pFq[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */
+    if (Debug) oxprintf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",A_pFq[0],B_pFq[0],M_n,Mapprox);
+  }
   mh_fclose(fp);
   return(0);
 }
-
-static int showParam(struct SFILE *fp,int fd) {
+#ifdef STANDALONE
+/* may remove */
+static int showParam_v1(struct SFILE *fp,int fd) {
   int rank,i;
   char swork[1024];
   if (fd) {
@@ -1576,7 +2103,9 @@ static int showParam(struct SFILE *fp,int fd) {
   for (i=0; i<Mg; i++) {
     sprintf(swork,"%%Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp);
   }
-  sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp);
+  if (*Ng >= 0) {
+    sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp);
+  }
   sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp);
   for (i=0; i<rank; i++) {
     sprintf(swork,"%%Iv[%d]=\n%lg\n",i,Iv[i]); mh_fputs(swork,fp);
@@ -1588,22 +2117,140 @@ static int showParam(struct SFILE *fp,int fd) {
   sprintf(swork,"%%Hg=\n%lf\n",Hg); mh_fputs(swork,fp);
   sprintf(swork,"%%Dp=\n%d\n",Dp);  mh_fputs(swork,fp);
   sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp);
+  
+  sprintf(swork,"%%%% Optional paramters\n"); mh_fputs(swork,fp);
+  sprintf(swork,"#success=%d\n",M_mh_t_success); mh_fputs(swork,fp);
+  sprintf(swork,"#automatic=%d\n",M_automatic); mh_fputs(swork,fp);
+  sprintf(swork,"#series_error=%lg\n",M_series_error); mh_fputs(swork,fp);
+  sprintf(swork,"#recommended_abserr\n"); mh_fputs(swork,fp);
+  sprintf(swork,"#abserror=%lg\n",M_recommended_abserr); mh_fputs(swork,fp);
+  if (M_recommended_relerr < MH_RELERR_DEFAULT) {
+    sprintf(swork,"%%relerror=%lg\n",M_recommended_relerr); mh_fputs(swork,fp);
+  }
+  sprintf(swork,"#mh_t_value=%lg # Value of matrix hg at X0g.\n",M_mh_t_value); mh_fputs(swork,fp);
+  sprintf(swork,"# M_m=%d  # Approximation degree of matrix hg.\n",M_m); mh_fputs(swork,fp);
+  sprintf(swork,"#beta_i_x_o2_max=%lg #max(|beta[i]*x|/2)\n",M_beta_i_x_o2_max); mh_fputs(swork,fp);
+  sprintf(swork,"#beta_i_beta_j_min=%lg #min(|beta[i]-beta[j]|)\n",M_beta_i_beta_j_min); mh_fputs(swork,fp);
+  sprintf(swork,"# change # to %% to read as an optional parameter.\n"); mh_fputs(swork,fp);
+  sprintf(swork,"%%p_pFq=%d, ",A_LEN); mh_fputs(swork,fp);
+  for (i=0; i<A_LEN; i++) {
+    if (i != A_LEN-1) sprintf(swork," %lg,",A_pFq[i]); 
+    else sprintf(swork," %lg\n",A_pFq[i]); 
+    mh_fputs(swork,fp);
+  }
+  sprintf(swork,"%%q_pFq=%d, ",B_LEN); mh_fputs(swork,fp);
+  for (i=0; i<B_LEN; i++) {
+    if (i != B_LEN-1) sprintf(swork," %lg,",B_pFq[i]); 
+    else sprintf(swork," %lg\n",B_pFq[i]); 
+    mh_fputs(swork,fp);
+  }
+  sprintf(swork,"%%ef_type=%d\n",Ef_type); mh_fputs(swork,fp);
   return(0);
 }
+#endif
+/* version2.0 format */
+static int showParam(struct SFILE *fp,int fd) {
+  int rank,i;
+  char swork[1024];
+  if (fd) {
+    fp = mh_fopen("stdout","w",1);
+  }
+  rank = imypower(2,Mg);
+  sprintf(swork,"%s\n",VSTRING); mh_fputs(swork,fp);
+  sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp);
+  sprintf(swork,"%%p_pFq=%d, ",A_LEN); mh_fputs(swork,fp);
+  for (i=0; i<A_LEN; i++) {
+    if (i != A_LEN-1) sprintf(swork," %lg,",A_pFq[i]); 
+    else sprintf(swork," %lg\n",A_pFq[i]); 
+    mh_fputs(swork,fp);
+  }
+  sprintf(swork,"%%q_pFq=%d, ",B_LEN); mh_fputs(swork,fp);
+  for (i=0; i<B_LEN; i++) {
+    if (i != B_LEN-1) sprintf(swork," %lg,",B_pFq[i]); 
+    else sprintf(swork," %lg\n",B_pFq[i]); 
+    mh_fputs(swork,fp);
+  }
+  sprintf(swork,"%%ef_type=%d\n",Ef_type); mh_fputs(swork,fp);
+  sprintf(swork,"%%Beta=\n"); mh_fputs(swork,fp);
+  for (i=0; i<Mg; i++) {
+    sprintf(swork,"#Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp);
+  }
+  if (*Ng >= 0) {
+    sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp);
+  }
+  sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp);
+  sprintf(swork,"%%Iv=\n"); mh_fputs(swork,fp);
+  for (i=0; i<rank; i++) {
+    sprintf(swork,"#Iv[%d]=\n%lg\n",i,Iv[i]); mh_fputs(swork,fp);
+    if (Sample && (M_n == 2) && (X0g == 0.3)) {
+      sprintf(swork,"%%Iv[%d]-Iv2[%d]=%lg\n",i,i,Iv[i]-Iv2[i]); mh_fputs(swork,fp);
+    }
+  }
+  sprintf(swork,"%%Ef=\n%lg\n",Ef); mh_fputs(swork,fp);
+  sprintf(swork,"%%Hg=\n%lf\n",Hg); mh_fputs(swork,fp);
+  sprintf(swork,"%%Dp=\n%d\n",Dp);  mh_fputs(swork,fp);
+  sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp);
+  
+  sprintf(swork,"%%%% Optional paramters\n"); mh_fputs(swork,fp);
+  sprintf(swork,"#success=%d\n",M_mh_t_success); mh_fputs(swork,fp);
+  sprintf(swork,"#automatic=%d\n",M_automatic); mh_fputs(swork,fp);
+  sprintf(swork,"#series_error=%lg\n",M_series_error); mh_fputs(swork,fp);
+  sprintf(swork,"#recommended_abserr\n"); mh_fputs(swork,fp);
+  sprintf(swork,"#abserror=%lg\n",M_recommended_abserr); mh_fputs(swork,fp);
+  if (M_recommended_relerr < MH_RELERR_DEFAULT) {
+    sprintf(swork,"%%relerror=%lg\n",M_recommended_relerr); mh_fputs(swork,fp);
+  }
+  sprintf(swork,"#mh_t_value=%lg # Value of matrix hg at X0g.\n",M_mh_t_value); mh_fputs(swork,fp);
+  sprintf(swork,"# M_m=%d  # Approximation degree of matrix hg.\n",M_m); mh_fputs(swork,fp);
+  sprintf(swork,"#beta_i_x_o2_max=%lg #max(|beta[i]*x|/2)\n",M_beta_i_x_o2_max); mh_fputs(swork,fp);
+  sprintf(swork,"#beta_i_beta_j_min=%lg #min(|beta[i]-beta[j]|)\n",M_beta_i_beta_j_min); mh_fputs(swork,fp);
+  sprintf(swork,"# change # to %% to read as an optional parameter.\n"); mh_fputs(swork,fp);
+  return(0);
+}
 
 static double gammam(double a,int n) {
   double v,v2;
   int i;
-  v=mypower(sqrt(M_PI),(n*(n-1))/2);
+  v=mypower(sqrt(M_PI),(n*(n-1))/2); /* pi^(n*(n-1)/2) */
   v2=0;
   for (i=1; i<=n; i++) {
     v2 += lgamma(a-((double)(i-1))/2.0); /* not for big n */
   }
-  if (Debug) printf("gammam(%lg,%d)=%lg\n",a,n,v*exp(v2));
+  if (Debug) oxprintf("gammam(%lg,%d)=%lg\n",a,n,v*exp(v2));
   return(v*exp(v2));
 }
 
 static double iv_factor(void) {
+  double ef;
+  double lef;
+  double r;
+  ef=1; lef=0;
+  if (Ef_type < 1) return(1.0);
+  if (Ef_type == 1) {
+	ef=iv_factor_ef_type1();
+	lef=liv_factor_ef_type1();
+  }else if (Ef_type==2) {
+	ef = iv_factor_ef_type2();
+	lef = liv_factor_ef_type2();
+  }else{
+	return(1.0);
+  }
+  if (isnan(ef) || isinf(ef)) {
+	if (Debug) oxprintfe("Exponential factor (Ef) seems to be large or ill-conditioned.\n");
+	return(exp(lef));
+  }else {
+	r = ef/exp(lef);
+	if ((0.9 < r) && (r < 1.1)) return(ef);
+	else {
+	  oxprintfe("Warning: There seems to be a numerical error to get Ef. We use a log value of Ef");
+      oxprintfe(" Ef=%lg, exp(lef)=%lg\n",ef,exp(lef));
+	  return(exp(lef));
+	}
+  }
+  return(exp(lef));
+}
+
+static double iv_factor_ef_type1(void) {
   double v1;
   double t;
   double b;
@@ -1624,5 +2271,108 @@ static double iv_factor(void) {
   return( c*v1);
 }
 
+static double iv_factor_ef_type2(void) {
+  double ef;
+  int i,m;
+  double a,b,c;
+  double t;
+  a = A_pFq[0]; b = A_pFq[1]; c= B_pFq[0];
+  m = M_n;
+  ef = 1.0;
 
+  /*fprintf(stderr,"iv_factor_ef_type2: a=%lf, b=%lf, c=%lf, m=%d\n",a,b,c,m);*/
+  /* Ref: note 2016.02.04 */
+  if (X0g<0){ myerror("Negative X0g\n"); mh_exit(-1);}
+  t = 0;
+  for (i=0; i<m; i++) if (Beta[i]<0){ myerror("Negative beta\n"); mh_exit(-1);}
+  for (i=0; i<m; i++) t += log(Beta[i]);
+  ef *= exp((a+b-c)*t);
 
+  t = 0;
+  for (i=0; i<m; i++) t += log(Beta[i]+X0g);
+  ef *= exp(-b*t);
+
+  ef *= exp((c-a)*(2*a-1)*log(X0g));
+
+  ef *= gammam(b,m)/gammam(a+b-c,m);
+  ef *= gammam(a,m)/gammam(c,m);
+  return(ef);  
+}
+static void setM_x(void) {
+  if (Ef_type <= 1)    return(setM_x_ef_type1());
+  else if (Ef_type==2) return(setM_x_ef_type2());
+  setM_x_ef_type1();
+}
+static void setM_x_ef_type1(void) {
+  int i;
+  for (i=0; i<M_n; i++) {
+	M_x[i] = Beta[i]*X0g;
+  }
+}
+static void setM_x_ef_type2(void) {
+  int i;
+  for (i=0; i<M_n; i++) {
+	M_x[i] = X0g/(Beta[i]+X0g);
+  }
+}
+
+/* log of gammam */
+static double lgammam(double a,int n) {
+  double v,v2;
+  int i;
+  v=log(M_PI)*n*(n-1)/2.0; /* log pi^(n*(n-1)/2) */
+  v2=0;
+  for (i=1; i<=n; i++) {
+    v2 += lgamma(a-((double)(i-1))/2.0); /* not for big n */
+  }
+  return(v+v2);
+}
+
+/* log of iv_factor_ef_type1() */
+static double liv_factor_ef_type1(void) {
+  double v1;
+  double t;
+  double b;
+  double detSigma;
+  double c;
+  int i,n;
+  if (X0g<0){ myerror("Negative X0g\n"); mh_exit(-1);}
+  n = (int) (*Ng);
+  v1= log(X0g)*n*M_n/2.0;
+  t = 0.0;
+  for (i=0; i<M_n; i++)  t += -X0g*Beta[i];
+  v1 += t;
+
+  b = 1; for (i=0; i<M_n; i++) b *= Beta[i];
+  detSigma = -log(b)-M_n*log(2);
+
+  c = lgammam(((double)(M_n+1))/2.0, M_n)-log(2)*M_n*n/2.0
+	-detSigma*n/2.0-lgammam(((double)(M_n+n+1))/2.0,M_n);
+  return(c+v1);
+}
+
+static double liv_factor_ef_type2(void) {
+  double ef;
+  int i,m;
+  double a,b,c;
+  double t;
+  a = A_pFq[0]; b = A_pFq[1]; c= B_pFq[0];
+  m = M_n;
+  ef = 0.0;
+
+  if (X0g<0){ myerror("Negative X0g\n"); mh_exit(-1);}
+  t = 0;
+  for (i=0; i<m; i++) if (Beta[i]<0){ myerror("Negative beta\n"); mh_exit(-1);}
+  for (i=0; i<m; i++) t += log(Beta[i]);
+  ef += (a+b-c)*t;
+
+  t = 0;
+  for (i=0; i<m; i++) t += log(Beta[i]+X0g);
+  ef += -b*t;
+
+  ef += (c-a)*(2*a-1)*log(X0g);
+
+  ef += lgammam(b,m)-lgammam(a+b-c,m);
+  ef += lgammam(a,m)-lgammam(c,m);
+  return(ef);  
+}