OpenXM/src/hgm/mh/src/jack-n.c - diff

Return to jack-n.c CVS log

Up to [local] / OpenXM / src / hgm / mh / src

Diff for /OpenXM/src/hgm/mh/src/jack-n.c between version 1.24 and 1.54

version 1.24, 2014/03/16 00:00:46

version 1.54, 2020/02/06 05:18:12

Line 4

#include <math.h>

#include <string.h>

#include "sfile.h"

#define VSTRING "%!version2.0"

$OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.23 2014/03/15 11:23:58 takayama Exp $

$OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.53 2016/10/30 01:10:18 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.

Line 15

Line 17

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

Line 46 static int Sample = Sample_default;

Line 53 static int Sample = Sample_default;

static double *Iv2;

static double Ef2;

static int SAR_warning = 1;

#ifdef NAN

#else

#define NAN 3.40282e+38 /* for old 32 bit machines. Todo, configure */

Line 55 static double Ef2;

Line 64 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;

Line 100 static double M_rel_error=0.0; /* relative errors */

Line 115 static double M_rel_error=0.0; /* relative errors */

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;

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.

Line 112 static double M_series_error;

Line 127 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;

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-30;

static double M_x0value_min=1e-60;

estimated_X0g is the suggested value of X0g.

Line 138 static int M_mh_t_success=1;

Line 153 static int M_mh_t_success=1;

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;

Line 149 static double M_beta_i_beta_j_min;

Line 169 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);

Line 162 static int mysum(int L[]);

Line 185 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);

Line 196 static int genDarray2(int M,int N);

Line 212 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.

Line 234 int jk_freeWorkArea() {

Line 313 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);

}

Line 278 int jk_initializeWorkArea() {

Line 358 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);}

#ifdef STANDALONE2

static int myerror(char *s) { oxprintfe("%s: type in control-C\n",s); getchar(); getchar(); return(0);}

#else

static int myerror(char *s) {

oxprintfe("Error in jack-n.c: %s\n",s);

mh_exit(-1);

return(0);

}

#endif

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);

Line 308 static double jack1(int K) {

Line 395 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);

Line 321 static double jack1diff(int K) {

Line 408 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);}

Line 341 static double xval(int ii,int p) { /* x_i^p */

Line 427 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]);

Line 390 static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m]

Line 476 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};

Line 397 static int test_ptrans() {

Line 484 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)

Line 411 static int huk(int K[],int I,int J) {

Line 498 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);

Line 426 static int hdk(int K[],int I,int J) {

Line 513 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);

Line 436 static int hdk(int K[],int I,int J) {

Line 523 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++) {

Line 446 static double jjk(int K[]) {

Line 533 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);

}

Line 492 static double mypower(double x,int n) {

Line 579 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;

Line 502 static double qk(int K[],double A[A_LEN],double B[B_LE

Line 589 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);

if (P > Q) {

for (I=0; I<Q; I++) V = V/ppoch(B[I],K);

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);

}

Line 520 static int bb(int N[],int K[],int M[],int I,int J) {

Line 617 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(K); oxprintf("K<--, "); printp2(Kp); oxprintf("<--Kp\n");

printp(M); printf("M<--, "); printp2(Mp); printf("<--Mp\n");

printp(M); oxprintf("M<--, "); printp2(Mp); oxprintf("<--Mp\n");

if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J));

Line 542 static double beta(int K[],int M[]) {

Line 639 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); */

}

Line 551 static double beta(int K[],int M[]) {

Line 648 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);*/

}

Line 559 static double beta(int K[],int M[]) {

Line 656 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]);

if (i <M_n-1) oxprintf("%d,",kappa[i]);

else printf("%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]);

if (i <ell+1-1) oxprintf("%d,",kappa[i]);

else printf("%d)",kappa[i]);

else oxprintf("%d)",kappa[i]);

}

return(0);

}

Line 582 static int test_beta() {

Line 680 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); oxprintf(","); printp(mu3); oxprintf(": beta = %lf\n",beta(kappa,mu3));

printp(kappa); printf(","); printp(mu1); printf(": beta = %lf\n",beta(kappa,mu1));

printp(kappa); oxprintf(","); printp(mu1); oxprintf(": beta = %lf\n",beta(kappa,mu1));

printp(kappa); printf(","); printp(mu2); printf(": beta = %lf\n",beta(kappa,mu2));

printp(kappa); oxprintf(","); printp(mu2); oxprintf(": beta = %lf\n",beta(kappa,mu2)); return(0);

}

#endif

/* main() { test_beta(); } */

Line 612 static int test_beta() {

Line 710 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];

Line 638 static double q3_10(int K[],int M[],int SK) {

Line 736 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];

Line 670 static double q3_5(double A[],double B[],int K[],int I

Line 769 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 A[1] = {1.5};

double B[B_LEN]={6.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};

Line 690 static void mtest4b() {

Line 791 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(); } */

Line 711 static int nk(int KK[]) {

Line 813 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);

Line 723 static int plength2(int P1[],int P2[]) {

Line 826 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);

Line 750 static int pListPartition(int M,int N) {

Line 854 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 = []; */

Line 768 static int pListPartition(int M,int N) {

Line 872 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 $B$KBP$7$F$d$k;E;v$r$3$A$i$X=q$/(B.

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");

}

Line 831 static int pListHS2(int From,int To,int Kap[]) {

Line 935 static int pListHS2(int From,int To,int Kap[]) {

}

static void hsExec_0() {

int i;

/* int i; */

if(Debug) {printf("hsExec: "); printp(HS_mu); printf("\n");}

if(Debug) {oxprintf("hsExec: "); printp(HS_mu); oxprintf("\n");}

}

Line 870 static int pmn(int M,int N) {

Line 974 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));

for (j=0; j<=M; j++) oxprintf("%d,",aP_pki(i,j));

printf("\n");

oxprintf("\n");

}

myfree(P_pki); P_pki=NULL;

Line 881 static int pmn(int M,int N) {

Line 985 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[]) {

Line 922 static int genDarray2(int M,int N) {

Line 1026 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);

Line 941 static int genDarray2(int M,int N) {

Line 1045 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;

Line 951 static int genDarray2(int M,int N) {

Line 1055 static int genDarray2(int M,int N) {

}

if (Debug) {

printf("Darray=\n");

oxprintf("Darray=\n");

for (i=0; i<Pmn; i++) printf("%d\n",Darray[i]);

for (i=0; i<Pmn; i++) oxprintf("%d\n",Darray[i]);

printf("-----------\n");

oxprintf("-----------\n");

}

return(0);

}

Line 982 static void hsExec_beta(void) {

Line 1086 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) {

Line 1006 static void hsExec_beta(void) {

Line 1110 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(Nu); oxprintf("\n");

printp(Mu); printf("\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;

Line 1030 static void hsExec_beta(void) {

Line 1134 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);

}

Line 1043 static int genBeta(int Kap[]) {

Line 1147 static int genBeta(int Kap[]) {

extern int M_beta_pt;

extern int M_beta_kap[];

extern int P_pmn;

int I,J,N;

int I,N;

if (Debug) {printp(Kap); printf("<-Kappa, P_pmn=%d\n",P_pmn);}

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));

Line 1053 static int genBeta(int Kap[]) {

Line 1157 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};

Line 1073 static int checkBeta1() {

Line 1177 static int checkBeta1() {

Mu = Parray[M_beta_1[I]];

Beta_km = M_beta_0[I];

if (Debug) {

printp(Kap); printf("<--Kap, ");

printp(Kap); oxprintf("<--Kap, ");

printp(Mu); printf("<--Mu,");

printp(Mu); oxprintf("<--Mu,");

printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap,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(Kap2); oxprintf("<--Kap, ");

printp(Mu); printf("<--Mu,");

printp(Mu); oxprintf("<--Mu,");

printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap2,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);

oxprintf("Mu=%a,",Mu);

printf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));

oxprintf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));

}

Line 1145 static int genJack(int M,int N) {

Line 1249 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;

int Nk,JJ, two_to_I;

if (Debug) printf("genJack(%d,%d)\n",M,N);

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.

Line 1169 static int genJack(int M,int N) {

Line 1273 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;

}

Line 1191 static int genJack(int M,int N) {

Line 1295 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 */

Line 1206 static int genJack(int M,int N) {

Line 1310 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);

if (Debug) {oxprintf("Nv(number of variables)=%d, Beta_km=%lf, Mu=",Nv,Beta_km);

printp(Mu); printf("\n");}

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) {

Line 1226 static int genJack(int M,int N) {

Line 1330 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);

if (Debug) {oxprintf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);

printp(Mu); printf("\n"); }

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 */

Line 1239 static int genJack(int M,int N) {

Line 1343 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) {

Line 1263 static int checkJack1(int M,int N) {

Line 1367 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]);

for (I=1; I<=N; I++) oxprintf("%lf, ",M_x[I-1]);

printf("<--x\n");

oxprintf("<--x\n");

return(0);

}

/*main() { checkJack1(3,3); }*/

Line 1289 static int checkJack2(int M,int N) {

Line 1393 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]);

for (I=1; I<=N; I++) oxprintf("%lf, ",M_x[I-1]);

printf("<--x\n");

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));

}

Line 1308 static int checkJack2(int M,int N) {

Line 1412 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;

Line 1338 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

Line 1443 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) oxprintf("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 && (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;

Line 1350 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

Line 1455 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

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) 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;

}

Line 1360 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

Line 1465 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

}

for (K=0; K<=P_pmn; K++) {

printf("Kappa="); for (i=0; i<N; i++) printf("%d ",Parray[K][i]); printf("\n");

oxprintf("Kappa="); for (i=0; i<N; i++) oxprintf("%d ",Parray[K][i]); oxprintf("\n");

printf("ParraySize(%d)=%d (|kappa|), M_m=%d\n",K,ParraySize[K],M_m);

oxprintf("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]);

oxprintf("partial_sum[%d]=%lg\n",i,partial_sum[i]);

}

M_estimated_X0g = X0g;

Line 1377 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

Line 1482 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

M_series_error = serror;

M_recommended_abserr = iv*M_assigned_series_error;

M_recommended_relerr = M_series_error;

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);

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);

}

if (isnan(F)) myerror("F is nan. Make q0 smaller (or X0g smaller on the standalone system).\n");

if (!isnormal(iv)) myerror("F*Ef (initial value) is zero. Make q0 larger (or X0g larger on the standalone system).\n");

M_mh_t_value=F;

return(F);

}

Line 1403 double mh_t2(int J) {

Line 1513 double mh_t2(int J) {

return(F);

}

#ifdef STANDALONE

static int mtest1b() {

double A[1] = {1.5};

double B[1] = {1.5+5};

Line 1415 static int mtest1b() {

Line 1526 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

Line 1432 static int mtest1b() {

Line 1544 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;

Line 1444 static int imypower(int x,int n) {

Line 1556 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);

}

Line 1477 struct MH_RESULT *jk_main(int argc,char *argv[]) {

Line 1589 struct MH_RESULT *jk_main(int argc,char *argv[]) {

}

struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree) {

double *y0;

// double *y0;

double x0,xn;

// double x0,xn;

double ef;

// double ef;

int i,j,rank;

double a[1]; double b[1];

int i,j; // int i,j,rank;

extern double M_x[];

extern double *Beta;

extern int M_2n;

Line 1490 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

Line 1602 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

struct MH_RESULT *ans=NULL;

struct SFILE *ofp = NULL;

int idata=0;

JK_byFile = 1;

JK_byFile = 1; for (i=0; i<1024; i++) swork[i]=0;

jk_initializeWorkArea();

UseTable = 1;

Mapprox=6;

Line 1511 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

Line 1623 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

}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 */

Line 1522 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

Line 1634 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

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);

}

Line 1544 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

Line 1656 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

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);

Line 1553 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

Line 1669 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

myerror("Mg must be equal to M_n\n"); mh_exit(-1);

}

if (Debug) showParam(NULL,1);

for (i=0; i<M_n; i++) {

setM_x();

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]);

Line 1568 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

Line 1682 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

}

a[0] = ((double)Mg+1.0)/2.0;

mh_t(A_pFq,B_pFq,M_n,Mapprox);

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;

Line 1601 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

Line 1712 struct MH_RESULT *jk_main2(int argc,char *argv[],int a

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");

oxprintfe("Usages:\n");

fprintf(stderr,"hgm_jack-n [--idata input_data_file --x0 x0 --degree approxm]\n");

#include "usage-jack-n.h"

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," (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,"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 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] ./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);

}

Line 1667 static int setParamDefault() {

Line 1747 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=";

char *ng="%%Ng=";

int i;

// 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) {

if (strncmp(s,ng,5) != 0) {

fprintf(stderr,"Warning, there is no %%Ng= at the border of Beta's and Ng, s=%s\n",s);

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);

}

Line 1696 static int setParam(char *fname) {

Line 1785 static int setParam(char *fname) {

struct SFILE *fp;

int i;

struct mh_token tk;

int version;

if (fname == NULL) return(setParamDefault());

Sample = 0;

Sample = 0; for (i=0; i<SMAX; i++) s[i]=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);

Line 1713 static int setParam(char *fname) {

Line 1815 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);

Line 1744 static int setParam(char *fname) {

Line 1845 static int setParam(char *fname) {

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) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

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) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("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);

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) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("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);

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) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("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);

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) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("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);

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) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("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);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

M_X0g_bound = tk.dval;

continue;

}

fprintf(stderr,"Unknown ID at %s\n",s); mh_exit(-1);

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);

}

#ifdef STANDALONE

static int showParam(struct SFILE *fp,int fd) {

/* may remove */

static int showParam_v1(struct SFILE *fp,int fd) {

int rank,i;

char swork[1024];

if (fd) {

Line 1816 static int showParam(struct SFILE *fp,int fd) {

Line 2114 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);

Line 1835 static int showParam(struct SFILE *fp,int fd) {

Line 2135 static int showParam(struct SFILE *fp,int fd) {

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;

Line 1875 static double iv_factor(void) {

Line 2282 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) {setM_x_ef_type1(); return;}

else if (Ef_type==2) {setM_x_ef_type2(); return;}

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;

if (myabs(M_x[i]) > SERIES_ADMISSIBLE_RADIUS_TYPE1) {

if (SAR_warning) oxprintfe("Warning: evaluation point %lf for %d-th variable of the series 1F1 might be far from 0. Decrease q0 (or X0g for the standalone) if necessary.\n",M_x[i],i);

SAR_warning=0;

}

static void setM_x_ef_type2(void) {

int i;

for (i=0; i<M_n; i++) {

M_x[i] = X0g/(Beta[i]+X0g);

if (myabs(M_x[i]) > SERIES_ADMISSIBLE_RADIUS_TYPE2) {

if (SAR_warning) oxprintfe("Warning: evaluation point %lf for %d-th point of the series 2F1 is near 1. Decrease q0 (or X0g for the standalone).\n",M_x[i],i);

SAR_warning=0;

}

int reset_SAR_warning(int n) {

SAR_warning = n;

return(n);

}

/* 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);

}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>