===================================================================
RCS file: /home/cvs/OpenXM/src/hgm/mh/src/jack-n.c,v
retrieving revision 1.3
retrieving revision 1.25
diff -u -p -r1.3 -r1.25
--- OpenXM/src/hgm/mh/src/jack-n.c	2013/02/21 07:30:56	1.3
+++ OpenXM/src/hgm/mh/src/jack-n.c	2014/03/16 03:11:07	1.25
@@ -5,25 +5,29 @@
 #include <string.h>
 #include "sfile.h"
 /*
-$OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.2 2013/02/20 07:53:44 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.
-cf. 20-my-note-mh-jack-n.pdf,  /Movies/oxvh/priv-j1/2011-12-14-ke-mh-jack.mov
-Todo:
-1. Cash the transposes of partitions.
-2. Use the recurrence to obtain beta().
-3. Easier input data file format for mh-n.c
-Changelog:
-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.
-2011.12.15, mh.r --> jack-n.c
+  $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.24 2014/03/16 00:00:46 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.
+  cf. 20-my-note-mh-jack-n.pdf,  /Movies/oxvh/priv-j1/2011-12-14-ke-mh-jack.mov
+  Todo:
+  1. Cash the transposes of partitions.
+  2. Use the recurrence to obtain beta().
+  3. Easier input data file format for mh-n.c
+  Changelog:
+  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.
+  2011.12.15, mh.r --> jack-n.c
 */
 
 /****** from mh-n.c *****/
 static int JK_byFile=1;
+static int JK_deallocate=0;
 #define M_n_default 3
+#define Sample_default 1
 static int M_n=0;
 /* global variables. They are set in setParam() */
 static int Mg;  /* n */
@@ -36,7 +40,7 @@ static double Ef;   /* exponential factor at beta*x0 *
 static double Hg;   /* step size of rk defined in rk.c */
 static int Dp;      /* Data sampling period */
 static double Xng=0.0;   /* the last point */
-static int Sample = 1;
+static int Sample = Sample_default;
 
 /* for sample inputs */
 static double *Iv2; 
@@ -54,7 +58,7 @@ static double Ef2; 
 #define A_LEN  1 /* (a_1) , (a_1, ..., a_p)*/
 #define B_LEN  1 /* (b_1) */
 static int Debug = 0;
-static int Alpha = 2;;
+static int Alpha = 2;  /* 2 implies the zonal polynomial */
 static int *Darray = NULL;
 static int **Parray = NULL; /* array of partitions of size M_n */
 static int *ParraySize = NULL; /* length of each partitions */
@@ -91,10 +95,69 @@ 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=0;
+/* 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=0.00001;
+/*
+  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-30;
+/*
+  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;
+/*
+  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 myfree(void *p);
-static myerror(char *s);
+static int myfree(void *p);
+static int myerror(char *s);
 static double jack1(int K);
 static double jack1diff(int k);
 static double xval(int i,int p); /* x_i^p */
@@ -102,7 +165,7 @@ 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 test_ptrans();
+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[]);
@@ -112,9 +175,9 @@ static double mypower(double x,int n);
 static double qk(int K[],double A[A_LEN],double B[B_LEN]);
 static int bb(int N[],int K[],int M[],int I,int J);
 static double beta(int K[],int M[]);
-static printp(int kappa[]);
-static printp2(int kappa[]);
-static test_beta();
+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();
@@ -130,26 +193,26 @@ static int pListHS2(int From,int To,int Kap[]);
 static void hsExec_0();
 static int pmn(int M,int N);
 static int *cloneP(int a[]);
-static copyP(int p[],int a[]);
+static int copyP(int p[],int a[]);
 static void pExec_darray(void);
-static genDarray2(int M,int N);
-static isHStrip(int Kap[],int Nu[]);
+static int genDarray2(int M,int N);
+static int isHStrip(int Kap[],int Nu[]);
 static void hsExec_beta(void);
-static genBeta(int Kap[]);
-static checkBeta1();
+static int genBeta(int Kap[]);
+static int checkBeta1();
 static int psublen(int Kap[],int Mu[]);
-static genJack(int M,int N);
-static checkJack1(int M,int N);
-static checkJack2(int M,int N);
-static mtest1b();
+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 usage();
-static setParamDefault();
-static next(struct SFILE *fp,char *s,char *msg);
+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);
+static int showParam(struct SFILE *fp,int fd);
 static double iv_factor(void);
 static double gammam(double a,int n);
 static double mypower(double a,int n);
@@ -157,30 +220,38 @@ static double mypower(double a,int n);
 double mh_t(double A[A_LEN],double B[B_LEN],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();
 
 int jk_freeWorkArea() {
   /* bug, p in the cloneP will not be deallocated.
-	      Nk in genDarray2 will not be deallocated.
-   */
+     Nk in genDarray2 will not be deallocated.
+  */
   int i;
+  JK_deallocate=1;
   if (Darray) {myfree(Darray); Darray=NULL;}
   if (Parray) {myfree(Parray); Parray=NULL;}
   if (ParraySize) {myfree(ParraySize); ParraySize=NULL;}
   if (M_beta_0) {myfree(M_beta_0); M_beta_0=NULL;}
   if (M_beta_1) {myfree(M_beta_1); M_beta_1=NULL;}
   if (M_jack) {
-	for (i=0; M_jack[i] != NULL; i++) {
-	  myfree(M_jack[i]); M_jack[i] = NULL;
-	}
-	myfree(M_jack); M_jack=NULL;
+    for (i=0; M_jack[i] != NULL; i++) {
+      if (Debug) printf("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;}
+  JK_deallocate=0;
+  return(0);
 }
 int jk_initializeWorkArea() {
   int i,j;
+  JK_deallocate=1;
+  xval(0,0);
+  JK_deallocate=0;
   Darray=NULL;
   Parray=NULL;
   ParraySize=NULL;
@@ -191,6 +262,7 @@ int jk_initializeWorkArea() {
   for (i=0; i<M_nmx; i++) M_kap[i]=HS_mu[i]=0;
   for (i=0; i<M_nmx; i++) M_x[i]=0;
   for (i=0; i<M_nmx; i++) for (j=0; j<M_m_MAX; j++) Xarray[i][j]=0;
+  for (i=0; i<M_nmx; i++) M_beta_kap[i]=0;
   M_m=M_m_MAX-2; 
   Alpha = 2;
   HS_n=M_nmx; 
@@ -201,20 +273,27 @@ int jk_initializeWorkArea() {
   M_df=1;
   M_2n=0;
   M_rel_error=0.0;
+  Sample = Sample_default;
+  Xng=0.0;
+  M_n=0;
+  return(0);
 }
 
 static void *mymalloc(int size) {
   void *p;
   if (Debug) printf("mymalloc(%d)\n",size);
-  p = (void *)malloc(size);
+  p = (void *)mh_malloc(size);
   if (p == NULL) {
     fprintf(stderr,"No more memory.\n");
     mh_exit(-1);
   }
   return(p);
 }
-static myfree(void *p) { free(p); }
-static myerror(char *s) { fprintf(stderr,"%s: type in control-C\n",s); getchar(); getchar();}
+static int myfree(void *p) {
+  if (Debug) printf("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 double jack1(int K) {
   double F;
@@ -247,15 +326,16 @@ static double xval(int ii,int p) { /* x_i^p */
   extern double M_x[];
   double F;
   int i,j;
-  static init=0;
+  static int init=0;
+  if (JK_deallocate) { init=0; return(0.0);}
   if (!init) {
     for (i=1; i<=M_n; i++) {
       for (j=0; j<M_m_MAX; j++) {
-	if (j != 0) {
-	  Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1];
-	}else{
-	  Xarray[i-1][0] = 1;
-	}
+        if (j != 0) {
+          Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1];
+        }else{
+          Xarray[i-1][0] = 1;
+        }
       }
     }
     init = 1;
@@ -279,7 +359,7 @@ static int mysum(int L[]) {
 }
 
 /*
- (3,2,2,0,0) --> 3
+  (3,2,2,0,0) --> 3
 */
 static int plength(int P[]) {
   int I;
@@ -298,7 +378,7 @@ static int plength_t(int P[]) {  
 }
 
 /*  
- ptrans(P)  returns Pt
+    ptrans(P)  returns Pt
 */
 static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] */
   extern int M_m;
@@ -313,7 +393,7 @@ static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m]
   }
 }
 
-static test_ptrans() {
+static int test_ptrans() {
   extern int M_m;
   int p[M_n0]={5,3,2};
   int pt[10];
@@ -321,6 +401,7 @@ static test_ptrans() {
   M_m = 10;
   ptrans(p,pt);
   if (Debug) {for (i=0; i<10; i++) printf("%d,",pt[i]);  printf("\n");}
+  return(0);
 }
 
 
@@ -413,7 +494,7 @@ static double mypower(double x,int n) {
   return(v);
 }
 /* Q_kappa
-*/
+ */
 static double qk(int K[],double A[A_LEN],double B[B_LEN]) {
   extern int Alpha;
   int P,Q,I;
@@ -433,8 +514,8 @@ static double qk(int K[],double A[A_LEN],double B[B_LE
 }
 
 /*
- B^nu_{kappa,mu}(i,j)
- bb(N,K,M,I,J)
+  B^nu_{kappa,mu}(i,j)
+  bb(N,K,M,I,J)
 */
 static int bb(int N[],int K[],int M[],int I,int J) {
   int Kp[M_m_MAX]; int Mp[M_m_MAX];
@@ -442,8 +523,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); printf("K<--, "); printp2(Kp); printf("<--Kp\n");
+    printp(M); printf("M<--, "); printp2(Mp); printf("<--Mp\n");
   */
 
   if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J));
@@ -479,15 +560,16 @@ static double beta(int K[],int M[]) {
 
   return(V);
 }
-static printp(int kappa[]) {
+static int printp(int kappa[]) {
   int i;
   printf("(");
   for (i=0; i<M_n; i++) {
     if (i <M_n-1) printf("%d,",kappa[i]);
     else printf("%d)",kappa[i]);
   }
+  return(0);
 }
-static printp2(int kappa[]) {
+static int printp2(int kappa[]) {
   int i,ell;
   printf("(");
   ell = plength_t(kappa);
@@ -495,9 +577,10 @@ static printp2(int kappa[]) {
     if (i <ell+1-1) printf("%d,",kappa[i]);
     else printf("%d)",kappa[i]);
   }
+  return(0);
 }
 
-static test_beta() {
+static int test_beta() {
   int kappa[M_n0]={2,1,0};
   int mu1[M_n0]={1,0,0};
   int mu2[M_n0]={1,1,0};
@@ -511,20 +594,20 @@ static test_beta() {
 
 
 /*
- cf. w1m.rr 
+  cf. w1m.rr 
   matrix hypergeometric by jack
   N variables, up to degree M.
 */
 /* todo
-def mhgj(A,B,N,M) {
-  F = 0;
-  P = partition_a(N,M);
-  for (I=0; I<length(P); I++) {
-     K = P[I];
-     F += qk(K,A,B)*jack(K,N);
-  }
-  return(F);
-}
+   def mhgj(A,B,N,M) {
+   F = 0;
+   P = partition_a(N,M);
+   for (I=0; I<length(P); I++) {
+   K = P[I];
+   F += qk(K,A,B)*jack(K,N);
+   }
+   return(F);
+   }
 */
 
 
@@ -571,22 +654,22 @@ static double q3_5(double A[],double B[],int K[],int I
   V=1; 
 
   for (J=1; J<=P; J++)  {
-     V *= (A[J-1]+C);
+    V *= (A[J-1]+C);
   }
   for (J=1; J<=Q; J++) {
-     V /= (B[J-1]+C);
+    V /= (B[J-1]+C);
   }
 
   for (J=1; J<=K[I-1]-1; J++) {
-     Ej = D-J*Alpha+Kp[J-1];
-     Gj = Ej+1;
-     V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha));
+    Ej = D-J*Alpha+Kp[J-1];
+    Gj = Ej+1;
+    V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha));
   } 
   for (J=1; J<=I-1; J++) {
-     Fj=K[J-1]*Alpha-J-D;
-     Hj=Fj+Alpha;
-     Lj=Hj*Fj;
-     V *= (Lj-Fj)/(Lj+Hj);
+    Fj=K[J-1]*Alpha-J-D;
+    Hj=Fj+Alpha;
+    Lj=Hj*Fj;
+    V *= (Lj-Fj)/(Lj+Hj);
   }
   return(V);
 }
@@ -616,7 +699,7 @@ static void mtest4b() {
 /* main() { mtest4(); mtest4b(); } */
 
 /* nk in (4.1),
-*/
+ */
 static int nk(int KK[]) {
   extern int *Darray;
   int N,I,Ki;
@@ -689,8 +772,9 @@ static int pListPartition(int M,int N) {
 */
 static int pListPartition2(int Less,int From,int To, int M) {
   int I,R;
+  mh_check_intr(100);
   if (To < From) {
-     (*M_pExec)(); return(0);
+    (*M_pExec)(); return(0);
   }
   for (I=1; (I<=Less) && (I<=M) ; I++) {
     M_kap[From-1] = I;
@@ -704,23 +788,23 @@ static int pListPartition2(int Less,int From,int To, i
 */
 static void pExec_0() {
   if (Debug) {
-	printf("M_kap=");
-	printp(M_kap);
-	printf("\n");
+    printf("M_kap=");
+    printp(M_kap);
+    printf("\n");
   }
 }
 
 /* Test.
-  Compare pListPartition(4,3);  genDarray(4,3);
-  Compare pListPartition(5,3);  genDarray(5,3);
+   Compare pListPartition(4,3);  genDarray(4,3);
+   Compare pListPartition(5,3);  genDarray(5,3);
 
 */
 
 /*
-main() {
+  main() {
   M_pExec = pExec_0;
   pListPartition(5,3);
-}
+  }
 */
 
 
@@ -735,7 +819,7 @@ static int pListHS(int Kap[],int N) {
   HS_n = N;
   /* Clear HS_mu. Do not forget when N < M_n  */
   for (i=0; i<M_n; i++) HS_mu[i] = 0;
-  pListHS2(1,N,Kap);
+  return(pListHS2(1,N,Kap));
 }
 
 static int pListHS2(int From,int To,int Kap[]) {
@@ -758,15 +842,15 @@ static void hsExec_0() {
   pListHS([0,4,2,1],3);
 */
 /*
-main() {
+  main() {
   int Kap[3]={4,2,1};
   HS_hsExec = hsExec_0;
   pListHS(Kap,3);
-}
+  }
 */
 
 /* The number of partitions <= M, with N parts.
-  (0,0,...,0) is excluded.
+   (0,0,...,0) is excluded.
 */
 #define aP_pki(i,j) P_pki[(i)*(M+1)+(j)]
 static int pmn(int M,int N) {
@@ -800,7 +884,7 @@ static int pmn(int M,int N) {
 }
 
 /*
-main() {pmn(4,3); printf("P_pmn=%d\n",P_pmn);}
+  main() {pmn(4,3); printf("P_pmn=%d\n",P_pmn);}
 */
 
 static int *cloneP(int a[]) {
@@ -810,9 +894,10 @@ static int *cloneP(int a[]) {
   for (i=0; i<M_n; i++) p[i] = a[i];
   return(p);
 }
-static copyP(int p[],int a[]) {
+static int copyP(int p[],int a[]) {
   int i;
   for (i=0; i<M_n; i++) p[i] = a[i];
+  return(0);
 }
 
 static void pExec_darray(void) {
@@ -827,7 +912,7 @@ static void pExec_darray(void) {
   ParraySize[DR_parray] = mysum(K);
   DR_parray++;
 }
-static genDarray2(int M,int N) {
+static int genDarray2(int M,int N) {
   extern int *Darray;
   extern int **Parray;
   extern int DR_parray;
@@ -855,30 +940,32 @@ static genDarray2(int M,int N) {
   Nk = (int *) mymalloc(sizeof(int)*(Pmn+1));
   for (I=0; I<Pmn; I++) Nk[I] = I;
   for (I=0; I<Pmn; I++) {
-     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");}
-	   continue;
-     }
-     for (i=0; i<M_n; i++) Kone[i] = 0;
-     for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1;
-     for (J=0; J<Pmn; J++) {
-        if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */
-     }
+    mh_check_intr(100);
+    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");}
+      continue;
+    }
+    for (i=0; i<M_n; i++) Kone[i] = 0;
+    for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1;
+    for (J=0; J<Pmn; J++) {
+      if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */
+    }
   }
   if (Debug) {
-	printf("Darray=\n");
-	for (i=0; i<Pmn; i++) printf("%d\n",Darray[i]);
-	printf("-----------\n");
+    printf("Darray=\n");
+    for (i=0; i<Pmn; i++) printf("%d\n",Darray[i]);
+    printf("-----------\n");
   }
+  return(0);
 }
 
 
 /* main() {  genDarray2(4,3);}  */
 
 /* M_beta_0[*] value of beta_{kappa,mu}, M_beta_1[*] N_mu */
-static isHStrip(int Kap[],int Nu[]) {
+static int isHStrip(int Kap[],int Nu[]) {
   int N1,N2,I,P;
   N1 = plength(Kap); N2 = plength(Nu);
   if (N2 > N1) return(0);
@@ -920,40 +1007,40 @@ static void hsExec_beta(void) {
     Done=0;
     for (J=M_beta_pt-1; J>=0; J--) {
       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",
-				  J,K,q3_10(M_beta_kap,Nu,K));
-		   printp(Nu); printf("\n");
-		   printp(Mu); printf("\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");
-			 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; }
-         /* check done. */
-         M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K);
-         Done = 1; break;
+        K=I+1;
+        if (Debug) {
+          printf("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");
+        }
+        /* 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");
+          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; }
+        /* check done. */
+        M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K);
+        Done = 1; break;
       } 
     }
     if (Done) break; else Nu[I]--;
   }
   if (!Done) {
     if (Debug) printf("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] = NAN;  error("Not found."); */
     M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu); 
   }
   /* Fix the bug of mh.rr */
   M_beta_pt++;
 }
-static genBeta(int Kap[]) {
+static int genBeta(int Kap[]) {
   extern double *M_beta_0;
   extern int *M_beta_1;
   extern int M_beta_pt;
@@ -963,7 +1050,7 @@ static genBeta(int Kap[]) {
   if (Debug) {printp(Kap); printf("<-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(double)*(P_pmn+1));
+  M_beta_1 = (int *)mymalloc(sizeof(int)*(P_pmn+1));
   M_beta_pt = 0;
   for (I=0; I<=P_pmn; I++) {M_beta_0[I] = NAN; M_beta_1[I] = -1;}
   N = plength(Kap);
@@ -977,7 +1064,7 @@ static genBeta(int Kap[]) {
   genBeta([2,1,1]);
 */
 
-static checkBeta1() {
+static int checkBeta1() {
   int Kap[3] = {2,2,0};
   int Kap2[3] = {2,1,0};
   int I;
@@ -988,38 +1075,39 @@ static checkBeta1() {
   for (I=0; I<M_beta_pt; I++) {
     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));
-	}
+    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));
+    }
   }
   if (Debug) printf("-------------------------------------\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));
-	}
+    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));
+    }
   }
+  return(0);
 }
 
 /*
-def checkBeta2() {
+  def checkBeta2() {
   genDarray2(3,3);
   Kap = [2,1,0];
   printf("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));
+  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));
   }
-}
+  }
 */
 
 /* main() { checkBeta1(); } */
@@ -1040,19 +1128,19 @@ static int psublen(int Kap[],int Mu[]) {
 
 
 /* Table of Jack polynomials
-  Jack[1]* one variable.
-  Jack[2]* two variables.
+   Jack[1]* one variable.
+   Jack[2]* two variables.
    ...
-  Jack[M_n]* n variables.
-  Jack[P][J]*  
-       D^J(P variables jack of p variables). Example. J=001 d_1, 010 d_2, 100 d_3
-       0<=J<=2^{M_n}-1
-  Jack[P][J][nk(Kappa)]  Jack_Kappa, Kappa is a partition.
-  0<=nk(Kappa)<=pmn(M_m,M_n)  
+   Jack[M_n]* n variables.
+   Jack[P][J]*  
+   D^J(P variables jack of p variables). Example. J=001 d_1, 010 d_2, 100 d_3
+   0<=J<=2^{M_n}-1
+   Jack[P][J][nk(Kappa)]  Jack_Kappa, Kappa is a partition.
+   0<=nk(Kappa)<=pmn(M_m,M_n)  
 */
 
 #define aM_jack(i,j,k) ((M_jack[i])[(j)*(Pmn+1)+(k)])
-static genJack(int M,int N) {
+static int genJack(int M,int N) {
   extern double **M_jack;
   extern int M_2n;
   extern int P_pmn;
@@ -1061,6 +1149,7 @@ static genJack(int M,int N) {
   int *Kap,*Mu;
   double Jack,Beta_km;
   int Nk,JJ;
+  if (Debug) printf("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. 
@@ -1089,7 +1178,7 @@ static genJack(int M,int N) {
       if (M_df) {
         for (J=1; J<M_2n; J++) {
           if (J >= 2^I) aM_jack(I,J,K) = 0;
-	  else aM_jack(I,J,K) = NAN;
+          else aM_jack(I,J,K) = NAN;
         }
       }
     }
@@ -1114,43 +1203,46 @@ static genJack(int M,int N) {
       for (H=0; H<M_beta_pt; H++) {
         Nk = M_beta_1[H];
         Mu = Parray[Nk];
-	if (UseTable) {
-	  Beta_km = M_beta_0[H];
-	}else{
-	  Beta_km = beta(Kap,Mu);
-  /* do not use the M_beta table. It's buggy. UseTable is experimental.*/
-	}
+        if (UseTable) {
+          Beta_km = M_beta_0[H];
+        }else{
+          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");}
+          printp(Mu); printf("\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));
       }
       aM_jack(Nv,0,K) = Jack;
       if (M_df) {
         /* The case of M_df > 0. */ 
         for (J=1; J<M_2n; J++) {
-      Jack = 0;
-      for (H=0; H<M_beta_pt; H++) {
-        Nk = M_beta_1[H];
-        Mu = Parray[Nk];
-	if (UseTable) {
-	  Beta_km = M_beta_0[H];
-	}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"); }
-        P = psublen(Kap,Mu);
-        if (J & (1 << (Nv-1))) {
-          JJ = J & ((1 << (Nv-1)) ^ 0xffff);  /* NOTE!! Up to 16 bits. mh-15 */
-	  if (P != 0) {
-	    Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1);
-	  }
-        }else{
-          Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P);
-        }
-      }
-      aM_jack(Nv,J,K) = Jack;
+          mh_check_intr(100);
+          Jack = 0;
+          for (H=0; H<M_beta_pt; H++) {
+            Nk = M_beta_1[H];
+            Mu = Parray[Nk];
+            if (UseTable) {
+              Beta_km = M_beta_0[H];
+            }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"); }
+            P = psublen(Kap,Mu);
+            if (J & (1 << (Nv-1))) {
+              JJ = J & ((1 << (Nv-1)) ^ 0xffff);  /* NOTE!! Up to 16 bits. mh-15 */
+              if (P != 0) {
+                Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1);
+              }
+            }else{
+              Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P);
+            }
+          }
+          aM_jack(Nv,J,K) = Jack;
+          if (Debug) printf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);
         } /* end of J loop */
       }
     }
@@ -1159,8 +1251,8 @@ static genJack(int M,int N) {
 
 
 /* checkJack1(3,3) 
-*/
-static checkJack1(int M,int N) {
+ */
+static int checkJack1(int M,int N) {
   int I,K;
   extern int P_pmn;
   extern double M_x[];
@@ -1179,15 +1271,16 @@ static checkJack1(int M,int N) {
   }
   for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);
   printf("<--x\n");
+  return(0);
 }
 /*main() {  checkJack1(3,3);  }*/
 
 
-static checkJack2(int M,int N) {
+static int checkJack2(int M,int N) {
   int I,K,J;
   extern int P_pmn;
   extern double M_x[];
-  extern M_df;
+  extern int M_df;
   int Pmn; /* used in aM_jack */
   M_df=1;
   /* initialize x vars. */
@@ -1208,12 +1301,13 @@ static checkJack2(int M,int 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",
-			   I,J,aM_jack(I,J,K));
-	  }
-	}
+        printp(Parray[K]);
+        printf("<--Kap, Nv=%d,J(diff)=%d, D^J Jack=%lf\n",
+               I,J,aM_jack(I,J,K));
+      }
+    }
   }
+  return(0);
 }
 
 /* main() { checkJack2(3,3); } */
@@ -1224,25 +1318,79 @@ double mh_t(double A[A_LEN],double B[B_LEN],int N,int 
   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);
-  M_qk = (double *)mymalloc(sizeof(double)*P_pmn);
+  M_qk = (double *)mymalloc(sizeof(double)*(P_pmn+1)); /* found a bug. */
   Pmn = P_pmn;
   size = ParraySize[P_pmn];
   for (K=0; K<=P_pmn; K++) {
-	Kap = Parray[K];
-	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));
+    mh_check_intr(100);
+    Kap = Parray[K];
+    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));
   }
   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];
+    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++) {
+    printf("Kappa="); for (i=0; i<N; i++) printf("%d ",Parray[K][i]); printf("\n"); 
+    printf("ParraySize(%d)=%d (|kappa|),   M_m=%d\n",K,ParraySize[K],M_m);
+  }
+  for (i=0; i<=M_m; i++) {
+    printf("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;
+
+  if (M_show_autosteps) {
+    printf("%%%%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);
+    printf("%%%%x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
+    printf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
+    fprintf(stderr,"%%%%(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);
+    fprintf(stderr,"%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
+    fprintf(stderr,"%%%%(stderr) F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
+  }
+
+  M_mh_t_value=F;
   return(F);
 }
 double mh_t2(int J) {
@@ -1255,12 +1403,12 @@ double mh_t2(int J) {
   F = 0; 
   Pmn = P_pmn;
   for (K=0; K<P_pmn; K++) {
-     F += M_qk[K]*aM_jack(M_n,J,K);
+    F += M_qk[K]*aM_jack(M_n,J,K);
   }
   return(F);
 }
 
-static mtest1b() {
+static int mtest1b() {
   double A[1] = {1.5};
   double B[1] = {1.5+5};
   int I,N,M,J;
@@ -1271,8 +1419,8 @@ static 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);
+    F=mh_t2(J);
+    printf("J=%d, D^J mh_t=%lf\n",J,F);
   }
 }
 
@@ -1300,16 +1448,40 @@ static int imypower(int x,int n) {
   return(v);
 }
 
-#ifdef STANDALONE
+#ifdef STANDALONE2
 main(int argc,char *argv[]) {
   mh_exit(MH_RESET_EXIT);
-/*  jk_main(argc,argv);
-	printf("second run.\n"); */
+  /*  jk_main(argc,argv);
+      printf("second run.\n"); */
   jk_main(argc,argv);
+  return(0);
 }
 #endif
 
 struct MH_RESULT *jk_main(int argc,char *argv[]) {
+  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;
@@ -1318,46 +1490,64 @@ struct MH_RESULT *jk_main(int argc,char *argv[]) {
   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;
   int idata=0;
   JK_byFile = 1;
   jk_initializeWorkArea();
-  Debug = 0;
   UseTable = 1;
   Mapprox=6;
   for (i=1; i<argc; i++) {
-	if (strcmp(argv[i],"--idata")==0) {
-	  inci(i);
-	  setParam(argv[i]); idata=1;
-	}else if (strcmp(argv[i],"--degree")==0) {
-	  inci(i);
-	  sscanf(argv[i],"%d",&Mapprox);
-	}else if (strcmp(argv[i],"--x0")==0) {
-	  inci(i);
-	  sscanf(argv[i],"%lg",&X0g);
-	}else if (strcmp(argv[i],"--notable")==0) {
-	  UseTable = 0;
-	}else if (strcmp(argv[i],"--debug")==0) {
-	  Debug = 1;
-	}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);}
-	  JK_byFile = 0;
-	}else {
-	  fprintf(stderr,"Unknown option %s\n",argv[i]);
-	  usage();
-	  return(NULL);
-	}
+    if (strcmp(argv[i],"--idata")==0) {
+      inci(i);
+      setParam(argv[i]); idata=1;
+    }else if (strcmp(argv[i],"--degree")==0) {
+      inci(i);
+      sscanf(argv[i],"%d",&Mapprox);
+    }else if (strcmp(argv[i],"--x0")==0) {
+      inci(i);
+      sscanf(argv[i],"%lg",&X0g);
+    }else if (strcmp(argv[i],"--notable")==0) {
+      UseTable = 0;
+    }else if (strcmp(argv[i],"--debug")==0) {
+      Debug = 1;
+    }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);}
+      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]);
+      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;
   /* Output by a file=stdout */
   ofp = mh_fopen("stdout","w",JK_byFile);
@@ -1365,69 +1555,101 @@ struct MH_RESULT *jk_main(int argc,char *argv[]) {
   sprintf(swork,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp);
   sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp);
   if (M_n != Mg) {
-	myerror("Mg must be equal to M_n\n"); mh_exit(-1);
+    myerror("Mg must be equal to M_n\n"); mh_exit(-1);
   }
-  if (Debug) showParam(ofp);
+  if (Debug) showParam(NULL,1);
   for (i=0; i<M_n; i++) {
-	M_x[i] = Beta[i]*X0g;
+    M_x[i] = Beta[i]*X0g;
   } 
+
+  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++) {
+    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]);
+    }
+  }
+
   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 ((!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);
+    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);
   }
   sprintf(swork,"%%%%M_rel_error=%lg\n",M_rel_error); mh_fputs(swork,ofp);
   for (j=0; j<M_2n; j++) {
-	Iv[j] = mh_t2(j);
+    Iv[j] = mh_t2(j);
   }
   Ef = iv_factor();
-  showParam(ofp);
+  showParam(ofp,0);
 
   /* return the result */
-  if (!JK_byFile) {
-	ans->x = X0g;
-	ans->rank = imypower(2,Mg);
-	ans->y = (double *)mymalloc(sizeof(double)*(ans->rank));
-	for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i];
-	ans->size=1;
-	ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size));
-	(ans->sfpp)[0] = ofp;
+  if (!JK_byFile) { 
+    ans->x = X0g;
+    ans->rank = imypower(2,Mg);
+    ans->y = (double *)mymalloc(sizeof(double)*(ans->rank));
+    for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i];
+    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");
   jk_freeWorkArea();
+  if (Debug) printf("jk_freeWorkArea() has finished.\n");
   return(ans);
 }
 
-static usage() {
+static int usage() {
   fprintf(stderr,"Usages:\n");
-  fprintf(stderr,"mh-m [--idata input_data_file --x0 x0 --degree approxm]\n");
-  fprintf(stderr,"\nThe command mh-m [options] generates an input for w-m, 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 mh-m uses the Koev-Edelman algorithm to evalute the matrix hypergeometric function.\n");
+  fprintf(stderr,"hgm_jack-n [--idata input_data_file --x0 x0 --degree approxm]\n");
+  fprintf(stderr,"           [--automatic n --assigned_series_error e --x0value_min e2]\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," Iv: initial values at X0g*Beta (see our paper how to order them), are evaluated in this program. Give zeros. \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 w-m, X0g: starting value of x(when --x0 option is used, this value is used), Xng: terminating value of x which is for w-m.\n");
-  fprintf(stderr," Dp: output data is stored in every Dp steps when output_data_file is specified, which is for w-m.\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,"Optional parameters automatic, ... are interpreted by a parser. See setParam() in jack-n.c and Testdata/tmp-idata2.txt\n");
+  fprintf(stderr,"Parameters are redefined when they appear more than once in the idata file and the command line options.\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 mh-m with no option.\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,"By --automatic option, X0g and degree are automatically searched. The current strategy is described in mh_t in jack-n.c\n");
+  fprintf(stderr,"Default values for the papameters of the automatic mode: assigned_series_error=%lg, x0value_min=%lg\n",M_assigned_series_error,M_x0value_min);
+  fprintf(stderr,"Todo: automatic mode throws away the table of Jack polynomials of the previous degrees and reevaluate them. They should be kept.\n");
   fprintf(stderr,"\nExamples:\n");
-  fprintf(stderr,"[1] ./mh-2 \n");
-  fprintf(stderr,"[2] ./mh-2 --x0 0.3 \n");
-  fprintf(stderr,"[3] ./mh-3 --x0 0.1 \n");
-  fprintf(stderr,"[4] ./mh-3 --x0 0.1 --degree 15 \n");
-  fprintf(stderr,"[5] ./mh-3 --idata test.txt --degree 15 \n");
-  fprintf(stderr,"[6] ./mh-3 --degree 15 >test2.txt\n");
-  fprintf(stderr,"    ./w-3 --idata test2.txt --gnuplotf test-g\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");
+  fprintf(stderr,"[6] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --automatic 1\n");
+  return(0);
 }
 
-static setParamDefault() {
+static int setParamDefault() {
   int rank;
   int i;
   Mg = M_n_default ;
@@ -1450,29 +1672,41 @@ static setParamDefault() {
   Hg = 0.001;
   Dp = 1;
   Xng = 10.0;
+  return(0);
 }
 
-static next(struct SFILE *sfp,char *s,char *msg) {
+static int next(struct SFILE *sfp,char *s,char *msg) {
+  static int check=1;
+  char *ng="%Ng=";
+  int i;
   s[0] = '%';
   while (s[0] == '%') {
-	if (!mh_fgets(s,SMAX,sfp)) {
-	  fprintf(stderr,"Data format error at %s\n",msg);
-	  mh_exit(-1);
-	}
-	if (s[0] != '%') return(0);
+    if (!mh_fgets(s,SMAX,sfp)) {
+      fprintf(stderr,"Data format error at %s\n",msg);
+      mh_exit(-1);
+    }
+    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);
+      }
+      check=0;
+    }
+    if (s[0] != '%') return(0);
   }
+  return(0);
 }
-static setParam(char *fname) {
+static int setParam(char *fname) {
   int rank;
   char s[SMAX];
   struct SFILE *fp;
   int i;
+  struct mh_token tk;
   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);
-	mh_exit(-1);
+    if (JK_byFile) fprintf(stderr,"File %s is not found.\n",fp->s);
+    mh_exit(-1);
   }
   next(fp,s,"Mg(m)");
   sscanf(s,"%d",&Mg);
@@ -1481,24 +1715,29 @@ static setParam(char *fname) {
   Beta = (double *)mymalloc(sizeof(double)*Mg);
   for (i=0; i<Mg; i++) {
     next(fp,s,"Beta");
-	sscanf(s,"%lf",&(Beta[i]));
+    sscanf(s,"%lf",&(Beta[i]));
   }
 
   Ng = (double *)mymalloc(sizeof(double));
-  next(fp,s,"Ng(freedom parameter n)");
+  next(fp,s,"%Ng= (freedom parameter n)");
   sscanf(s,"%lf",Ng);
   
   next(fp,s,"X0g(initial point)");
   sscanf(s,"%lf",&X0g);
-  
+
   Iv = (double *)mymalloc(sizeof(double)*rank);
   for (i=0; i<rank; i++) {
-	next(fp,s,"Iv(initial values)");
-	sscanf(s,"%lg",&(Iv[i]));
+    next(fp,s,"Iv(initial values)");
+	if (strncmp(s,"*",1)==0) {
+	  for (i=0; i<rank; i++) Iv[i] = 0.0;
+	  break;
+	}
+    sscanf(s,"%lg",&(Iv[i]));
   }
 
   next(fp,s,"Ef(exponential factor)");
-  sscanf(s,"%lg",&Ef);
+  if (strncmp(s,"*",1)==0) Ef=0.0;
+  else sscanf(s,"%lg",&Ef);
 
   next(fp,s,"Hg (step size of rk)");
   sscanf(s,"%lf",&Hg);
@@ -1508,29 +1747,114 @@ static setParam(char *fname) {
 
   next(fp,s,"Xng (the last point, cf. --xmax)");
   sscanf(s,"%lf",&Xng);
+
+  /* Reading the optional parameters */
+  while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) {
+    /* expect ID */
+    if (tk.type != MH_TOKEN_ID) {
+      fprintf(stderr,"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) {
+        fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        fprintf(stderr,"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) {
+        fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
+        fprintf(stderr,"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) {
+        fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
+        fprintf(stderr,"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) {
+        fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        fprintf(stderr,"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) {
+        fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
+        fprintf(stderr,"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) {
+        fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
+      }
+      if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
+        fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
+      }
+      M_show_autosteps = tk.ival;
+      continue;
+    }
+    fprintf(stderr,"Unknown ID at %s\n",s); mh_exit(-1);
+  }
   mh_fclose(fp);
+  return(0);
 }
 
-static showParam(struct SFILE *fp) {
+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,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp);
   for (i=0; i<Mg; i++) {
-	sprintf(swork,"%%Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp);
+    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);
   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);
-	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,"%%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);
+  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);
+  return(0);
 }
 
 static double gammam(double a,int n) {
@@ -1562,7 +1886,7 @@ static double iv_factor(void) {
   detSigma = 1.0/(b*mypower(2.0,M_n));
 
   c = gammam(((double)(M_n+1))/2.0, M_n)/
-     ( mypower(sqrt(2), M_n*n)*mypower(sqrt(detSigma),n)*gammam(((double)(M_n+n+1))/2.0,M_n) );
+    ( mypower(sqrt(2), M_n*n)*mypower(sqrt(detSigma),n)*gammam(((double)(M_n+n+1))/2.0,M_n) );
   return( c*v1);
 }