===================================================================
RCS file: /home/cvs/OpenXM/src/hgm/mh/src/jack-n.c,v
retrieving revision 1.32
retrieving revision 1.39
diff -u -p -r1.32 -r1.39
--- OpenXM/src/hgm/mh/src/jack-n.c	2015/04/02 01:11:13	1.32
+++ OpenXM/src/hgm/mh/src/jack-n.c	2016/02/04 02:52:19	1.39
@@ -5,7 +5,7 @@
 #include <string.h>
 #include "sfile.h"
 /*
-  $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.31 2015/04/02 00:11:32 takayama Exp $
+  $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.38 2016/02/01 07:05:25 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 +15,8 @@
   2. Use the recurrence to obtain beta().
   3. Easier input data file format for mh-n.c
   Changelog:
+  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
@@ -55,8 +57,24 @@ 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 */
+#ifdef C_2F1
+#define A_LEN  2 /* (a_1) , (a_1, ..., a_p)*/
+#define B_LEN  1 /* (b_1) */
+static int P_pFq=2;
+static int Q_pFq=1;
+#else
 #define A_LEN  1 /* (a_1) , (a_1, ..., a_p)*/
 #define B_LEN  1 /* (b_1) */
+static int P_pFq=1;
+static int Q_pFq=1;
+#endif
+static double A_pFq[A_LEN];
+static double B_pFq[B_LEN];
+#ifndef C_2F1
+static int Orig_1F1=1;
+#else
+static int Orig_1F1=0;
+#endif
 static int Debug = 0;
 static int Alpha = 2;  /* 2 implies the zonal polynomial */
 static int *Darray = NULL;
@@ -446,9 +464,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++) {
@@ -512,11 +530,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);
 }
 
@@ -623,12 +651,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];
@@ -784,20 +812,20 @@ 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) {
@@ -1504,7 +1532,7 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a
   // double ef;
 
   int i,j; // int i,j,rank;
-  double a[1]; double b[1];
+  double a[A_LEN]; double b[B_LEN];
   extern double M_x[];
   extern double *Beta;
   extern int M_2n;
@@ -1572,7 +1600,7 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a
     if (M_n > Mapprox) Mapprox=M_n;
   }
   /* Output by a file=stdout */
-  ofp = mh_fopen("oxstdout","w",JK_byFile);
+  ofp = mh_fopen("stdout","w",JK_byFile);
 
   sprintf(swork,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp);
   sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp);
@@ -1595,9 +1623,19 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a
     }
   }
 
-  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) oxprintf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",a[0],b[0],M_n,Mapprox);
+  if ((P_pFq != A_LEN) || (Q_pFq != B_LEN)) {
+    oxprintfe("It must be P_pFq == A_LEN and Q_pFq == B_LEN in this version. %s\n","");
+    mh_exit(-1);
+  }
+  oxprintfe("%%%%(stderr) Orig_1F1=%d\n",Orig_1F1);
+  if ((P_pFq == 1) && (Q_pFq == 1) && (Orig_1F1)) {
+    A_pFq[0] = a[0] = ((double)Mg+1.0)/2.0;
+    B_pFq[0] = b[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[0],b[0],M_n,Mapprox);
+  }else{
+    for (i=0; i<P_pFq; i++) a[i] = A_pFq[i];
+    for (i=0; i<Q_pFq; i++) b[i] = B_pFq[i];
+  }
   mh_t(a,b,M_n,Mapprox);
   if ((!M_mh_t_success) && M_automatic) {
     jk_freeWorkArea();
@@ -1636,7 +1674,12 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a
 
 static int usage() {
   oxprintfe("Usages:\n");
-  oxprintfe("hgm_jack-n [--idata input_data_file --x0 x0 --degree approxm]\n");
+#ifdef C_2F1
+  oxprintfe("hgm_jack-n-2f1");
+#else
+  oxprintfe("hgm_jack-n    ");
+#endif
+  oxprintfe(" [--idata input_data_file --x0 x0 --degree approxm]\n");
   oxprintfe("           [--automatic n --assigned_series_error e --x0value_min e2]\n");
   oxprintfe("\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 matrices with n degrees of freedom and the covariantce matrix sigma.\n");
   oxprintfe("The hgm_jack-n uses the Koev-Edelman algorithm to evalute the matrix hypergeometric function.\n");
@@ -1657,6 +1700,9 @@ static int usage() {
   oxprintfe("An example format of the input_data_file can be obtained by executing hgm_jack-n with no option. When there is no --idata file, all options are ignored.\n");
   oxprintfe("By --automatic option, X0g and degree are automatically determined from assigend_series_error. The current strategy is described in mh_t in jack-n.c\n");
   oxprintfe("Default values for the papameters of the automatic mode: automatic=%d, assigned_series_error=%lg, x0value_min=%lg\n",M_automatic,M_assigned_series_error,M_x0value_min);
+#ifdef C_2F1
+  oxprintfe("The parameters a,b,c of 2F1 are given by %%p_pFq=2, a,b  and  %%q_pFq=1, c\nNg is ignored.\n");
+#endif
   oxprintfe("Todo: automatic mode throws away the table of Jack polynomials of the previous degrees and reevaluate them. They should be kept.\n");
   oxprintfe("\nExamples:\n");
   oxprintfe("[1] ./hgm_jack-n \n");
@@ -1666,6 +1712,9 @@ static int usage() {
   oxprintfe("    gnuplot -persist <test-g-gp.txt\n");
   oxprintfe("[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --automatic 1 --assigned_series_error=1e-12\n");
   oxprintfe("[5] ./hgm_jack-n --idata Testdata/tmp-idata4.txt\n");
+#ifdef C_2F1
+  oxprintfe("Todo for 2F1: example for hgm_jack-n-2f1 has not been written.\niv_factor? Ef?");
+#endif
   return(0);
 }
 
@@ -1691,7 +1740,13 @@ static int setParamDefault() {
   X0g = (Beta[0]/Beta[Mg-1])*0.5;
   Hg = 0.001;
   Dp = 1;
-  Xng = 10.0;
+  if ((P_pFq == 1) && (Q_pFq == 1)) {
+    Xng = 10.0;
+  }else {
+	Xng=0.25;
+	for (i=0; i<A_LEN; i++) A_pFq[i] = (i+1)/5.0;
+	for (i=0; i<B_LEN; i++) B_pFq[i] = (A_LEN+i+1)/7.0;
+  }
   return(0);
 }
 
@@ -1834,6 +1889,46 @@ static int setParam(char *fname) {
       M_show_autosteps = tk.ival;
       continue;
     }
+    // Format: #p_pFq=2  1.5  3.2
+    if (strcmp(s,"p_pFq")==0) {
+      Orig_1F1=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);
+      }
+      P_pFq = tk.ival;
+      for (i=0; i<P_pFq; 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);
+      }
+      Q_pFq = tk.ival;
+      for (i=0; i<Q_pFq; 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;
+    }
     oxprintfe("Unknown ID at %s\n",s); mh_exit(-1);
   }
   mh_fclose(fp);
@@ -1844,7 +1939,7 @@ static int showParam(struct SFILE *fp,int fd) {
   int rank,i;
   char swork[1024];
   if (fd) {
-    fp = mh_fopen("oxstdout","w",1);
+    fp = mh_fopen("stdout","w",1);
   }
   rank = imypower(2,Mg);
   sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp);
@@ -1869,7 +1964,7 @@ static int showParam(struct SFILE *fp,int fd) {
   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,"#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);
   }
@@ -1877,13 +1972,26 @@ static int showParam(struct SFILE *fp,int fd) {
   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, ",P_pFq); mh_fputs(swork,fp);
+  for (i=0; i<P_pFq; i++) {
+    if (i != P_pFq-1) sprintf(swork," %lg,",A_pFq[i]); 
+    else sprintf(swork," %lg\n",A_pFq[i]); 
+    mh_fputs(swork,fp);
+  }
+  sprintf(swork,"%%q_pFq=%d, ",Q_pFq); mh_fputs(swork,fp);
+  for (i=0; i<Q_pFq; i++) {
+    if (i != Q_pFq-1) sprintf(swork," %lg,",B_pFq[i]); 
+    else sprintf(swork," %lg\n",B_pFq[i]); 
+    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 */
@@ -1899,6 +2007,7 @@ static double iv_factor(void) {
   double detSigma;
   double c;
   int i,n;
+  if ((P_pFq != 1) || (Q_pFq != 1)) return(1.0);
   n = (int) (*Ng);
   v1= mypower(sqrt(X0g),n*M_n);
   t = 0.0;