version 1.3, 2013/02/21 07:30:56 |
version 1.48, 2016/05/30 00:38:18 |
|
|
#include <math.h> |
#include <math.h> |
#include <string.h> |
#include <string.h> |
#include "sfile.h" |
#include "sfile.h" |
|
|
|
#define VSTRING "%!version2.0" |
/* |
/* |
$OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.2 2013/02/20 07:53:44 takayama Exp $ |
$OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.47 2016/03/06 23:51:36 takayama Exp $ |
Ref: copied from this11/misc-2011/A1/wishart/Prog |
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 |
jack-n.c, translated from mh.rr or tk_jack.rr in the asir-contrib. License: LGPL |
Koev-Edelman for higher order derivatives. |
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 |
cf. 20-my-note-mh-jack-n.pdf, /Movies/oxvh/priv-j1/2011-12-14-ke-mh-jack.mov |
Todo: |
Todo: |
1. Cash the transposes of partitions. |
1. Cash the transposes of partitions. |
2. Use the recurrence to obtain beta(). |
2. Use the recurrence to obtain beta(). |
3. Easier input data file format for mh-n.c |
3. Easier input data file format for mh-n.c |
Changelog: |
Changelog: |
2012.02.21, porting to OpenXM/src/hgm |
2016.02.15 log of Ef |
2011.12.22, --table option, which is experimental. |
2016.02.09 unify 2F1 and 1F1. Parser. |
2011.12.19, bug fix. jjk() should return double. It can become more than max int. |
2016.02.04 Ef_type (exponential or scalar factor type) |
2011.12.15, mh.r --> jack-n.c |
2016.02.01 ifdef C_2F1 ... |
|
2016.01.12 2F1 |
|
2014.03.15 http://fe.math.kobe-u.ac.jp/Movies/oxvh/2014-03-11-jack-n-c-automatic see also hgm/doc/ref.html, @s/2014/03/15-my-note-jack-automatic-order-F_A-casta.pdf. |
|
2014.03.14, --automatic option. Output estimation data. |
|
2012.02.21, porting to OpenXM/src/hgm |
|
2011.12.22, --table option, which is experimental. |
|
2011.12.19, bug fix. jjk() should return double. It can become more than max int. |
|
2011.12.15, mh.r --> jack-n.c |
*/ |
*/ |
|
|
/****** from mh-n.c *****/ |
/****** from mh-n.c *****/ |
static int JK_byFile=1; |
static int JK_byFile=1; |
|
static int JK_deallocate=0; |
#define M_n_default 3 |
#define M_n_default 3 |
|
#define Sample_default 1 |
static int M_n=0; |
static int M_n=0; |
/* global variables. They are set in setParam() */ |
/* global variables. They are set in setParam() */ |
static int Mg; /* n */ |
static int Mg; /* n */ |
Line 36 static double Ef; /* exponential factor at beta*x0 * |
|
Line 47 static double Ef; /* exponential factor at beta*x0 * |
|
static double Hg; /* step size of rk defined in rk.c */ |
static double Hg; /* step size of rk defined in rk.c */ |
static int Dp; /* Data sampling period */ |
static int Dp; /* Data sampling period */ |
static double Xng=0.0; /* the last point */ |
static double Xng=0.0; /* the last point */ |
static int Sample = 1; |
static int Sample = Sample_default; |
|
|
/* for sample inputs */ |
/* for sample inputs */ |
static double *Iv2; |
static double *Iv2; |
static double Ef2; |
static double Ef2; |
|
|
|
static int SAR_warning = 1; |
|
|
#ifdef NAN |
#ifdef NAN |
#else |
#else |
#define NAN 3.40282e+38 /* for old 32 bit machines. Todo, configure */ |
#define NAN 3.40282e+38 /* for old 32 bit machines. Todo, configure */ |
Line 51 static double Ef2; |
|
Line 64 static double Ef2; |
|
#define M_n0 3 /* used for tests. Must be equal to M_n */ |
#define M_n0 3 /* used for tests. Must be equal to M_n */ |
#define M_m_MAX 200 |
#define M_m_MAX 200 |
#define M_nmx M_m_MAX /* maximal of M_n */ |
#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 Debug = 0; |
static int Alpha = 2;; |
static int Alpha = 2; /* 2 implies the zonal polynomial */ |
static int *Darray = NULL; |
static int *Darray = NULL; |
static int **Parray = NULL; /* array of partitions of size M_n */ |
static int **Parray = NULL; /* array of partitions of size M_n */ |
static int *ParraySize = NULL; /* length of each partitions */ |
static int *ParraySize = NULL; /* length of each partitions */ |
Line 91 static double Xarray[M_nmx][M_m_MAX]; |
|
Line 110 static double Xarray[M_nmx][M_m_MAX]; |
|
static double *M_qk=NULL; /* saves pochhammerb */ |
static double *M_qk=NULL; /* saves pochhammerb */ |
static double M_rel_error=0.0; /* relative errors */ |
static double M_rel_error=0.0; /* relative errors */ |
|
|
|
/* For automatic degree and X0g setting. */ |
|
/* |
|
If automatic == 1, then the series is reevaluated as long as t_success!=1 |
|
by increasing X0g (evaluation point) and M_m (approx degree); |
|
*/ |
|
static int M_automatic=1; |
|
/* Estimated degree bound for series expansion. See mh_t */ |
|
static int M_m_estimated_approx_deg=0; |
|
/* Let F(i) be the approximation up to degree i. |
|
The i-th series_error is defined |
|
by |(F(i)-F(i-1))/F(i-1)|. |
|
*/ |
|
static double M_series_error; |
|
/* |
|
M_series_error < M_assigend_series_error (A) is required for the |
|
estimated_approx_deg. |
|
*/ |
|
static double M_assigned_series_error=M_ASSIGNED_SERIES_ERROR_DEFAULT; |
|
/* |
|
Let Ef be the exponential factor ( Ef=(4)/1F1 of [HNTT] ) |
|
If F(M_m)*Ef < x0value_min (B), the success=0 and X0g is increased. |
|
Note that minimal double is about 2e-308 |
|
*/ |
|
static double M_x0value_min=1e-60; |
|
/* |
|
estimated_X0g is the suggested value of X0g. |
|
*/ |
|
static double M_estimated_X0g=0.0; |
|
/* |
|
X0g should be less than M_X0g_bound. |
|
*/ |
|
static double M_X0g_bound = 1e+100; |
|
/* |
|
success is set to 1 when (A) and (B) are satisfied. |
|
*/ |
|
static int M_mh_t_success=1; |
|
/* |
|
recommended_abserr is the recommended value of the absolute error |
|
for the Runge-Kutta method. It is defined as |
|
assigend_series_error(standing for significant digits)*Ig[0](initial value) |
|
*/ |
|
static double M_recommended_abserr; |
|
/* |
|
recommended_relerr is the recommended value of the relative error |
|
for the Runge-Kutta method.. |
|
*/ |
|
static double M_recommended_relerr; |
|
/* |
|
max of beta(i)*x/2 |
|
*/ |
|
static double M_beta_i_x_o2_max; |
|
/* |
|
minimum of |beta_i-beta_j| |
|
*/ |
|
static double M_beta_i_beta_j_min; |
|
/* |
|
Value of matrix hg |
|
*/ |
|
static double M_mh_t_value; |
|
/* |
|
Show the process of updating degree. |
|
*/ |
|
int M_show_autosteps=1; |
|
|
/* prototypes */ |
/* prototypes */ |
static void *mymalloc(int size); |
static void *mymalloc(int size); |
static myfree(void *p); |
static int myfree(void *p); |
static myerror(char *s); |
static int myerror(char *s); |
static double jack1(int K); |
static double jack1(int K); |
static double jack1diff(int k); |
static double jack1diff(int k); |
static double xval(int i,int p); /* x_i^p */ |
static double xval(int i,int p); /* x_i^p */ |
Line 102 static int mysum(int L[]); |
|
Line 185 static int mysum(int L[]); |
|
static int plength(int P[]); |
static int plength(int P[]); |
static int plength_t(int P[]); |
static int plength_t(int P[]); |
static void ptrans(int P[M_nmx],int Pt[]); |
static void ptrans(int P[M_nmx],int Pt[]); |
static test_ptrans(); |
|
static int huk(int K[],int I,int J); |
static int huk(int K[],int I,int J); |
static int hdk(int K[],int I,int J); |
static int hdk(int K[],int I,int J); |
static double jjk(int K[]); |
static double jjk(int K[]); |
static double ppoch(double A,int K[]); |
static double ppoch(double A,int K[]); |
static double ppoch2(double A,double B,int K[]); |
static double ppoch2(double A,double B,int K[]); |
static double mypower(double x,int n); |
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 int bb(int N[],int K[],int M[],int I,int J); |
static double beta(int K[],int M[]); |
static double beta(int K[],int M[]); |
static printp(int kappa[]); |
static int printp(int kappa[]); |
static printp2(int kappa[]); |
|
static test_beta(); |
|
static double q3_10(int K[],int M[],int SK); |
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 nk(int KK[]); |
static int plength2(int P1[],int P2[]); |
|
static int myeq(int P1[],int P2[]); |
static int myeq(int P1[],int P2[]); |
static int pListPartition(int M,int N); |
static int pListPartition(int M,int N); |
static int pListPartition2(int Less,int From,int To, int M); |
static int pListPartition2(int Less,int From,int To, int M); |
Line 130 static int pListHS2(int From,int To,int Kap[]); |
|
Line 206 static int pListHS2(int From,int To,int Kap[]); |
|
static void hsExec_0(); |
static void hsExec_0(); |
static int pmn(int M,int N); |
static int pmn(int M,int N); |
static int *cloneP(int a[]); |
static int *cloneP(int a[]); |
static copyP(int p[],int a[]); |
static int copyP(int p[],int a[]); |
static void pExec_darray(void); |
static void pExec_darray(void); |
static genDarray2(int M,int N); |
static int genDarray2(int M,int N); |
static isHStrip(int Kap[],int Nu[]); |
static int isHStrip(int Kap[],int Nu[]); |
static void hsExec_beta(void); |
static void hsExec_beta(void); |
static genBeta(int Kap[]); |
static int genBeta(int Kap[]); |
static checkBeta1(); |
|
static int psublen(int Kap[],int Mu[]); |
static int psublen(int Kap[],int Mu[]); |
static genJack(int M,int N); |
static int genJack(int M,int N); |
static checkJack1(int M,int N); |
|
static checkJack2(int M,int N); |
|
static mtest1b(); |
|
|
|
static int imypower(int x,int n); |
static int imypower(int x,int n); |
static usage(); |
static int usage(); |
static setParamDefault(); |
static int setParamDefault(); |
static next(struct SFILE *fp,char *s,char *msg); |
static int next(struct SFILE *fp,char *s,char *msg); |
static int gopen_file(void); |
|
static int setParam(char *fname); |
static int setParam(char *fname); |
static int showParam(struct SFILE *fp); |
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 iv_factor(void); |
static double gammam(double a,int n); |
static double gammam(double a,int n); |
static double mypower(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); |
double mh_t2(int J); |
struct MH_RESULT *jk_main(int argc,char *argv[]); |
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_freeWorkArea(); |
int jk_initializeWorkArea(); |
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() { |
int jk_freeWorkArea() { |
/* bug, p in the cloneP will not be deallocated. |
/* 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; |
int i; |
|
JK_deallocate=1; |
if (Darray) {myfree(Darray); Darray=NULL;} |
if (Darray) {myfree(Darray); Darray=NULL;} |
if (Parray) {myfree(Parray); Parray=NULL;} |
if (Parray) {myfree(Parray); Parray=NULL;} |
if (ParraySize) {myfree(ParraySize); ParraySize=NULL;} |
if (ParraySize) {myfree(ParraySize); ParraySize=NULL;} |
if (M_beta_0) {myfree(M_beta_0); M_beta_0=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_beta_1) {myfree(M_beta_1); M_beta_1=NULL;} |
if (M_jack) { |
if (M_jack) { |
for (i=0; M_jack[i] != NULL; i++) { |
for (i=0; M_jack[i] != NULL; i++) { |
myfree(M_jack[i]); M_jack[i] = NULL; |
if (Debug) oxprintf("Free M_jack[%d]\n",i); |
} |
myfree(M_jack[i]); M_jack[i] = NULL; |
myfree(M_jack); M_jack=NULL; |
} |
|
myfree(M_jack); M_jack=NULL; |
} |
} |
if (M_qk) {myfree(M_qk); M_qk=NULL;} |
if (M_qk) {myfree(M_qk); M_qk=NULL;} |
if (P_pki) {myfree(P_pki); P_pki=NULL;} |
if (P_pki) {myfree(P_pki); P_pki=NULL;} |
|
setA(NULL,-1); setB(NULL,-1); |
|
JK_deallocate=0; |
|
return(0); |
} |
} |
int jk_initializeWorkArea() { |
int jk_initializeWorkArea() { |
int i,j; |
int i,j; |
|
JK_deallocate=1; |
|
xval(0,0); |
|
JK_deallocate=0; |
Darray=NULL; |
Darray=NULL; |
Parray=NULL; |
Parray=NULL; |
ParraySize=NULL; |
ParraySize=NULL; |
Line 191 int jk_initializeWorkArea() { |
|
Line 339 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_kap[i]=HS_mu[i]=0; |
for (i=0; i<M_nmx; i++) M_x[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++) 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; |
M_m=M_m_MAX-2; |
Alpha = 2; |
Alpha = 2; |
HS_n=M_nmx; |
HS_n=M_nmx; |
Line 201 int jk_initializeWorkArea() { |
|
Line 350 int jk_initializeWorkArea() { |
|
M_df=1; |
M_df=1; |
M_2n=0; |
M_2n=0; |
M_rel_error=0.0; |
M_rel_error=0.0; |
|
Sample = Sample_default; |
|
Xng=0.0; |
|
M_n=0; |
|
return(0); |
} |
} |
|
|
static void *mymalloc(int size) { |
static void *mymalloc(int size) { |
void *p; |
void *p; |
if (Debug) printf("mymalloc(%d)\n",size); |
if (Debug) oxprintf("mymalloc(%d)\n",size); |
p = (void *)malloc(size); |
p = (void *)mh_malloc(size); |
if (p == NULL) { |
if (p == NULL) { |
fprintf(stderr,"No more memory.\n"); |
oxprintfe("No more memory.\n"); |
mh_exit(-1); |
mh_exit(-1); |
} |
} |
return(p); |
return(p); |
} |
} |
static myfree(void *p) { free(p); } |
static int myfree(void *p) { |
static myerror(char *s) { fprintf(stderr,"%s: type in control-C\n",s); getchar(); getchar();} |
if (Debug) oxprintf("myFree at %p\n",p); |
|
return(mh_free(p)); |
|
} |
|
static int myerror(char *s) { oxprintfe("%s: type in control-C\n",s); getchar(); getchar(); return(0);} |
|
|
static double jack1(int K) { |
static double jack1(int K) { |
double F; |
double F; |
extern int Alpha; |
extern int Alpha; |
int I,J,L,II,JJ,N; |
int J,II,JJ,N; /* int I,J,L,II,JJ,N; */ |
N = 1; |
N = 1; |
if (K == 0) return((double)1); |
if (K == 0) return((double)1); |
F = xval(1,K); |
F = xval(1,K); |
Line 232 static double jack1(int K) { |
|
Line 388 static double jack1(int K) { |
|
static double jack1diff(int K) { |
static double jack1diff(int K) { |
double F; |
double F; |
extern int Alpha; |
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; |
N = 1; |
if (K == 0) return((double)1); |
if (K == 0) return((double)1); |
F = K*xval(1,K-1); |
F = K*xval(1,K-1); |
Line 245 static double jack1diff(int K) { |
|
Line 401 static double jack1diff(int K) { |
|
|
|
static double xval(int ii,int p) { /* x_i^p */ |
static double xval(int ii,int p) { /* x_i^p */ |
extern double M_x[]; |
extern double M_x[]; |
double F; |
|
int i,j; |
int i,j; |
static init=0; |
static int init=0; |
|
if (JK_deallocate) { init=0; return(0.0);} |
if (!init) { |
if (!init) { |
for (i=1; i<=M_n; i++) { |
for (i=1; i<=M_n; i++) { |
for (j=0; j<M_m_MAX; j++) { |
for (j=0; j<M_m_MAX; j++) { |
if (j != 0) { |
if (j != 0) { |
Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1]; |
Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1]; |
}else{ |
}else{ |
Xarray[i-1][0] = 1; |
Xarray[i-1][0] = 1; |
} |
} |
} |
} |
} |
} |
init = 1; |
init = 1; |
Line 264 static double xval(int ii,int p) { /* x_i^p */ |
|
Line 420 static double xval(int ii,int p) { /* x_i^p */ |
|
if (p > M_m_MAX-2) myerror("xval, p is too large."); |
if (p > M_m_MAX-2) myerror("xval, p is too large."); |
if (p < 0) { |
if (p < 0) { |
myerror("xval, p is negative."); |
myerror("xval, p is negative."); |
printf("ii=%d, p=%d\n",ii,p); |
oxprintf("ii=%d, p=%d\n",ii,p); |
mh_exit(-1); |
mh_exit(-1); |
} |
} |
return(Xarray[ii-1][p]); |
return(Xarray[ii-1][p]); |
Line 279 static int mysum(int L[]) { |
|
Line 435 static int mysum(int L[]) { |
|
} |
} |
|
|
/* |
/* |
(3,2,2,0,0) --> 3 |
(3,2,2,0,0) --> 3 |
*/ |
*/ |
static int plength(int P[]) { |
static int plength(int P[]) { |
int I; |
int I; |
Line 298 static int plength_t(int P[]) { |
|
Line 454 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] */ |
static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] */ |
extern int M_m; |
extern int M_m; |
Line 313 static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] |
|
Line 469 static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] |
|
} |
} |
} |
} |
|
|
static test_ptrans() { |
#ifdef STANDALONE |
|
static int test_ptrans() { |
extern int M_m; |
extern int M_m; |
int p[M_n0]={5,3,2}; |
int p[M_n0]={5,3,2}; |
int pt[10]; |
int pt[10]; |
int i; |
int i; |
M_m = 10; |
M_m = 10; |
ptrans(p,pt); |
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 |
upper hook length |
h_kappa^*(K) |
h_kappa^*(K) |
Line 333 static int huk(int K[],int I,int J) { |
|
Line 491 static int huk(int K[],int I,int J) { |
|
int Kp[M_m_MAX]; |
int Kp[M_m_MAX]; |
int A,H; |
int A,H; |
A=Alpha; |
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); |
ptrans(K,Kp); |
H=Kp[J-1]-I+A*(K[I-1]-J+1); |
H=Kp[J-1]-I+A*(K[I-1]-J+1); |
return(H); |
return(H); |
Line 348 static int hdk(int K[],int I,int J) { |
|
Line 506 static int hdk(int K[],int I,int J) { |
|
int Kp[M_m_MAX]; |
int Kp[M_m_MAX]; |
int A,H; |
int A,H; |
A = Alpha; |
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); |
ptrans(K,Kp); |
H=Kp[J-1]-I+1+A*(K[I-1]-J); |
H=Kp[J-1]-I+1+A*(K[I-1]-J); |
return(H); |
return(H); |
Line 358 static int hdk(int K[],int I,int J) { |
|
Line 516 static int hdk(int K[],int I,int J) { |
|
*/ |
*/ |
static double jjk(int K[]) { |
static double jjk(int K[]) { |
extern int Alpha; |
extern int Alpha; |
int A,L,I,J; |
int L,I,J; |
double V; |
double V; |
A=Alpha; |
|
V=1; |
V=1; |
L=plength(K); |
L=plength(K); |
for (I=0; I<L; I++) { |
for (I=0; I<L; I++) { |
Line 368 static double jjk(int K[]) { |
|
Line 526 static double jjk(int K[]) { |
|
V *= huk(K,I+1,J+1)*hdk(K,I+1,J+1); |
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); |
return(V); |
} |
} |
/* |
/* |
Line 413 static double mypower(double x,int n) { |
|
Line 571 static double mypower(double x,int n) { |
|
return(v); |
return(v); |
} |
} |
/* Q_kappa |
/* 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; |
extern int Alpha; |
int P,Q,I; |
int P,Q,I; |
double V; |
double V; |
Line 424 static double qk(int K[],double A[A_LEN],double B[B_LE |
|
Line 582 static double qk(int K[],double A[A_LEN],double B[B_LE |
|
/* to reduce numerical errors, temporary. */ |
/* to reduce numerical errors, temporary. */ |
if (P == Q) { |
if (P == Q) { |
for (I=0; I<P; I++) V = V*ppoch2(A[I],B[I],K); |
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); |
return(V); |
} |
} |
|
|
/* |
/* |
B^nu_{kappa,mu}(i,j) |
B^nu_{kappa,mu}(i,j) |
bb(N,K,M,I,J) |
bb(N,K,M,I,J) |
*/ |
*/ |
static int bb(int N[],int K[],int M[],int I,int J) { |
static int bb(int N[],int K[],int M[],int I,int J) { |
int Kp[M_m_MAX]; int Mp[M_m_MAX]; |
int Kp[M_m_MAX]; int Mp[M_m_MAX]; |
Line 442 static int bb(int N[],int K[],int M[],int I,int J) { |
|
Line 610 static int bb(int N[],int K[],int M[],int I,int J) { |
|
ptrans(M,Mp); |
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)); |
if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J)); |
Line 464 static double beta(int K[],int M[]) { |
|
Line 632 static double beta(int K[],int M[]) { |
|
for (J=0; J<K[I]; J++) { |
for (J=0; J<K[I]; J++) { |
II = I+1; JJ = J+1; |
II = I+1; JJ = J+1; |
V *= (double)bb(K,K,M,II,JJ); |
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 473 static double beta(int K[],int M[]) { |
|
Line 641 static double beta(int K[],int M[]) { |
|
for (J=0; J<M[I]; J++) { |
for (J=0; J<M[I]; J++) { |
II = I+1; JJ = J+1; |
II = I+1; JJ = J+1; |
V /= (double)bb(M,K,M,II,JJ); |
V /= (double)bb(M,K,M,II,JJ); |
/* printf("[%d,%d,%lf]\n",I,J,V);*/ |
/* oxprintf("[%d,%d,%lf]\n",I,J,V);*/ |
} |
} |
} |
} |
|
|
return(V); |
return(V); |
} |
} |
static printp(int kappa[]) { |
static int printp(int kappa[]) { |
int i; |
int i; |
printf("("); |
oxprintf("("); |
for (i=0; i<M_n; i++) { |
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); |
} |
} |
static printp2(int kappa[]) { |
#ifdef STANDALONE |
|
static int printp2(int kappa[]) { |
int i,ell; |
int i,ell; |
printf("("); |
oxprintf("("); |
ell = plength_t(kappa); |
ell = plength_t(kappa); |
for (i=0; i<ell+1; i++) { |
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); |
} |
} |
|
|
static test_beta() { |
static int test_beta() { |
int kappa[M_n0]={2,1,0}; |
int kappa[M_n0]={2,1,0}; |
int mu1[M_n0]={1,0,0}; |
int mu1[M_n0]={1,0,0}; |
int mu2[M_n0]={1,1,0}; |
int mu2[M_n0]={1,1,0}; |
int mu3[M_n0]={2,0,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(); } */ |
/* main() { test_beta(); } */ |
|
|
|
|
/* |
/* |
cf. w1m.rr |
cf. w1m.rr |
matrix hypergeometric by jack |
matrix hypergeometric by jack |
N variables, up to degree M. |
N variables, up to degree M. |
*/ |
*/ |
/* todo |
/* todo |
def mhgj(A,B,N,M) { |
def mhgj(A,B,N,M) { |
F = 0; |
F = 0; |
P = partition_a(N,M); |
P = partition_a(N,M); |
for (I=0; I<length(P); I++) { |
for (I=0; I<length(P); I++) { |
K = P[I]; |
K = P[I]; |
F += qk(K,A,B)*jack(K,N); |
F += qk(K,A,B)*jack(K,N); |
} |
} |
return(F); |
return(F); |
} |
} |
*/ |
*/ |
|
|
|
|
Line 532 def mhgj(A,B,N,M) { |
|
Line 703 def mhgj(A,B,N,M) { |
|
static double q3_10(int K[],int M[],int SK) { |
static double q3_10(int K[],int M[],int SK) { |
extern int Alpha; |
extern int Alpha; |
int Mp[M_m_MAX]; |
int Mp[M_m_MAX]; |
int ML[M_nmx]; |
// int ML[M_nmx]; |
int N[M_nmx]; |
int N[M_nmx]; |
int i,R; |
int i,R; |
double T,Q,V,Ur,Vr,Wr; |
double T,Q,V,Ur,Vr,Wr; |
ptrans(M,Mp); |
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; |
N[SK-1] = N[SK-1]-1; |
|
|
T = SK-Alpha*M[SK-1]; |
T = SK-Alpha*M[SK-1]; |
Line 558 static double q3_10(int K[],int M[],int SK) { |
|
Line 729 static double q3_10(int K[],int M[],int SK) { |
|
return(V); |
return(V); |
} |
} |
|
|
|
#ifdef STANDALONE |
static double q3_5(double A[],double B[],int K[],int I) { |
static double q3_5(double A[],double B[],int K[],int I) { |
extern int Alpha; |
extern int Alpha; |
int Kp[M_m_MAX]; |
int Kp[M_m_MAX]; |
Line 571 static double q3_5(double A[],double B[],int K[],int I |
|
Line 743 static double q3_5(double A[],double B[],int K[],int I |
|
V=1; |
V=1; |
|
|
for (J=1; J<=P; J++) { |
for (J=1; J<=P; J++) { |
V *= (A[J-1]+C); |
V *= (A[J-1]+C); |
} |
} |
for (J=1; J<=Q; J++) { |
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++) { |
for (J=1; J<=K[I-1]-1; J++) { |
Ej = D-J*Alpha+Kp[J-1]; |
Ej = D-J*Alpha+Kp[J-1]; |
Gj = Ej+1; |
Gj = Ej+1; |
V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha)); |
V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha)); |
} |
} |
for (J=1; J<=I-1; J++) { |
for (J=1; J<=I-1; J++) { |
Fj=K[J-1]*Alpha-J-D; |
Fj=K[J-1]*Alpha-J-D; |
Hj=Fj+Alpha; |
Hj=Fj+Alpha; |
Lj=Hj*Fj; |
Lj=Hj*Fj; |
V *= (Lj-Fj)/(Lj+Hj); |
V *= (Lj-Fj)/(Lj+Hj); |
} |
} |
return(V); |
return(V); |
} |
} |
|
#endif |
|
#ifdef STANDALONE |
static void mtest4() { |
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 K[M_n0] = {3,2,0}; |
int I=2; |
int I=2; |
int Ki[M_n0]={3,1,0}; |
int Ki[M_n0]={3,1,0}; |
double V1,V2; |
double V1,V2; |
|
setA(A,1); setB(B,1); |
V1=q3_5(A,B,K,I); |
V1=q3_5(A,B,K,I); |
V2=qk(K,A,B)/qk(Ki,A,B); |
V2=qk(K,A,B)/qk(Ki,A,B); |
printf("%lf== %lf?\n",V1,V2); |
oxprintf("%lf== %lf?\n",V1,V2); |
} |
} |
static void mtest4b() { |
static void mtest4b() { |
int K[M_n0]={3,2,0}; |
int K[M_n0]={3,2,0}; |
Line 610 static void mtest4b() { |
|
Line 784 static void mtest4b() { |
|
double V1,V2; |
double V1,V2; |
V1=q3_10(K,M,SK); |
V1=q3_10(K,M,SK); |
V2=beta(K,N)/beta(K,M); |
V2=beta(K,N)/beta(K,M); |
printf("%lf== %lf?\n",V1,V2); |
oxprintf("%lf== %lf?\n",V1,V2); |
} |
} |
|
#endif |
|
|
/* main() { mtest4(); mtest4b(); } */ |
/* main() { mtest4(); mtest4b(); } */ |
|
|
/* nk in (4.1), |
/* nk in (4.1), |
*/ |
*/ |
static int nk(int KK[]) { |
static int nk(int KK[]) { |
extern int *Darray; |
extern int *Darray; |
int N,I,Ki; |
int N,I,Ki; |
Line 631 static int nk(int KK[]) { |
|
Line 806 static int nk(int KK[]) { |
|
/* K = (Kpp,Ki) */ |
/* K = (Kpp,Ki) */ |
return(Darray[nk(Kpp)]+Ki-1); |
return(Darray[nk(Kpp)]+Ki-1); |
} |
} |
|
#ifdef STANDALONE |
static int plength2(int P1[],int P2[]) { |
static int plength2(int P1[],int P2[]) { |
int S1,S2; |
int S1,S2; |
S1 = plength(P1); S2 = plength(P2); |
S1 = plength(P1); S2 = plength(P2); |
Line 643 static int plength2(int P1[],int P2[]) { |
|
Line 819 static int plength2(int P1[],int P2[]) { |
|
} |
} |
else return(-1); |
else return(-1); |
} |
} |
|
#endif |
static int myeq(int P1[],int P2[]) { |
static int myeq(int P1[],int P2[]) { |
int I,L1; |
int I,L1; |
if ((L1=plength(P1)) != plength(P2)) return(0); |
if ((L1=plength(P1)) != plength(P2)) return(0); |
Line 670 static int pListPartition(int M,int N) { |
|
Line 847 static int pListPartition(int M,int N) { |
|
int I; |
int I; |
/* initialize */ |
/* initialize */ |
if (M_n != N) { |
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_m = M; |
/* M_plist = []; */ |
/* M_plist = []; */ |
Line 688 static int pListPartition(int M,int N) { |
|
Line 865 static int pListPartition(int M,int N) { |
|
Less >= M_kap[From], ..., M_kap[To], |(M_kap[From],...,M_kap[To])|<=M, |
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) { |
static int pListPartition2(int Less,int From,int To, int M) { |
int I,R; |
int I; |
|
mh_check_intr(100); |
if (To < From) { |
if (To < From) { |
(*M_pExec)(); return(0); |
(*M_pExec)(); return(0); |
} |
} |
for (I=1; (I<=Less) && (I<=M) ; I++) { |
for (I=1; (I<=Less) && (I<=M) ; I++) { |
M_kap[From-1] = I; |
M_kap[From-1] = I; |
R = pListPartition2(I,From+1,To,M-I); |
pListPartition2(I,From+1,To,M-I); |
} |
} |
return(1); |
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() { |
static void pExec_0() { |
if (Debug) { |
if (Debug) { |
printf("M_kap="); |
oxprintf("M_kap="); |
printp(M_kap); |
printp(M_kap); |
printf("\n"); |
oxprintf("\n"); |
} |
} |
} |
} |
|
|
/* Test. |
/* Test. |
Compare pListPartition(4,3); genDarray(4,3); |
Compare pListPartition(4,3); genDarray(4,3); |
Compare pListPartition(5,3); genDarray(5,3); |
Compare pListPartition(5,3); genDarray(5,3); |
|
|
*/ |
*/ |
|
|
/* |
/* |
main() { |
main() { |
M_pExec = pExec_0; |
M_pExec = pExec_0; |
pListPartition(5,3); |
pListPartition(5,3); |
} |
} |
*/ |
*/ |
|
|
|
|
Line 735 static int pListHS(int Kap[],int N) { |
|
Line 913 static int pListHS(int Kap[],int N) { |
|
HS_n = N; |
HS_n = N; |
/* Clear HS_mu. Do not forget when N < M_n */ |
/* Clear HS_mu. Do not forget when N < M_n */ |
for (i=0; i<M_n; i++) HS_mu[i] = 0; |
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[]) { |
static int pListHS2(int From,int To,int Kap[]) { |
Line 750 static int pListHS2(int From,int To,int Kap[]) { |
|
Line 928 static int pListHS2(int From,int To,int Kap[]) { |
|
} |
} |
|
|
static void hsExec_0() { |
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");} |
} |
} |
|
|
/* |
/* |
pListHS([0,4,2,1],3); |
pListHS([0,4,2,1],3); |
*/ |
*/ |
/* |
/* |
main() { |
main() { |
int Kap[3]={4,2,1}; |
int Kap[3]={4,2,1}; |
HS_hsExec = hsExec_0; |
HS_hsExec = hsExec_0; |
pListHS(Kap,3); |
pListHS(Kap,3); |
} |
} |
*/ |
*/ |
|
|
/* The number of partitions <= M, with N parts. |
/* 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)] |
#define aP_pki(i,j) P_pki[(i)*(M+1)+(j)] |
static int pmn(int M,int N) { |
static int pmn(int M,int N) { |
Line 789 static int pmn(int M,int N) { |
|
Line 967 static int pmn(int M,int N) { |
|
} |
} |
P_pmn=S; |
P_pmn=S; |
if (Debug) { |
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 (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; |
myfree(P_pki); P_pki=NULL; |
Line 800 static int pmn(int M,int N) { |
|
Line 978 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[]) { |
static int *cloneP(int a[]) { |
Line 810 static int *cloneP(int a[]) { |
|
Line 988 static int *cloneP(int a[]) { |
|
for (i=0; i<M_n; i++) p[i] = a[i]; |
for (i=0; i<M_n; i++) p[i] = a[i]; |
return(p); |
return(p); |
} |
} |
static copyP(int p[],int a[]) { |
static int copyP(int p[],int a[]) { |
int i; |
int i; |
for (i=0; i<M_n; i++) p[i] = a[i]; |
for (i=0; i<M_n; i++) p[i] = a[i]; |
|
return(0); |
} |
} |
|
|
static void pExec_darray(void) { |
static void pExec_darray(void) { |
Line 827 static void pExec_darray(void) { |
|
Line 1006 static void pExec_darray(void) { |
|
ParraySize[DR_parray] = mysum(K); |
ParraySize[DR_parray] = mysum(K); |
DR_parray++; |
DR_parray++; |
} |
} |
static genDarray2(int M,int N) { |
static int genDarray2(int M,int N) { |
extern int *Darray; |
extern int *Darray; |
extern int **Parray; |
extern int **Parray; |
extern int DR_parray; |
extern int DR_parray; |
Line 840 static genDarray2(int M,int N) { |
|
Line 1019 static genDarray2(int M,int N) { |
|
|
|
M_m = M; |
M_m = M; |
Pmn = pmn(M,N)+1; |
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); |
Darray=(int *) mymalloc(sizeof(int)*Pmn); |
for (i=0; i<Pmn; i++) Darray[i] = 0; |
for (i=0; i<Pmn; i++) Darray[i] = 0; |
Parray=(int **) mymalloc(sizeof(int *)*Pmn); |
Parray=(int **) mymalloc(sizeof(int *)*Pmn); |
Line 855 static genDarray2(int M,int N) { |
|
Line 1034 static genDarray2(int M,int N) { |
|
Nk = (int *) mymalloc(sizeof(int)*(Pmn+1)); |
Nk = (int *) mymalloc(sizeof(int)*(Pmn+1)); |
for (I=0; I<Pmn; I++) Nk[I] = I; |
for (I=0; I<Pmn; I++) Nk[I] = I; |
for (I=0; I<Pmn; I++) { |
for (I=0; I<Pmn; I++) { |
K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */ |
mh_check_intr(100); |
Ksize = plength(K); |
K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */ |
if (Ksize >= M_n) { |
Ksize = plength(K); |
if (Debug) {fprintf(stderr,"Ksize >= M_n\n");} |
if (Ksize >= M_n) { |
continue; |
if (Debug) {oxprintfe("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 (i=0; i<M_n; i++) Kone[i] = 0; |
for (J=0; J<Pmn; J++) { |
for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1; |
if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */ |
for (J=0; J<Pmn; J++) { |
} |
if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */ |
|
} |
} |
} |
if (Debug) { |
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); |
} |
} |
|
|
|
|
/* main() { genDarray2(4,3);} */ |
/* main() { genDarray2(4,3);} */ |
|
|
/* M_beta_0[*] value of beta_{kappa,mu}, M_beta_1[*] N_mu */ |
/* 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; |
int N1,N2,I,P; |
N1 = plength(Kap); N2 = plength(Nu); |
N1 = plength(Kap); N2 = plength(Nu); |
if (N2 > N1) return(0); |
if (N2 > N1) return(0); |
Line 898 static void hsExec_beta(void) { |
|
Line 1079 static void hsExec_beta(void) { |
|
int Nu[M_nmx]; |
int Nu[M_nmx]; |
int rrMax; |
int rrMax; |
hsExec_0(); |
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 = cdr(vtol(HS_mu)); */ |
Mu = HS_mu; /* buggy? need cloneP */ |
Mu = HS_mu; /* buggy? need cloneP */ |
if (M_beta_pt == 0) { |
if (M_beta_pt == 0) { |
Line 920 static void hsExec_beta(void) { |
|
Line 1101 static void hsExec_beta(void) { |
|
Done=0; |
Done=0; |
for (J=M_beta_pt-1; J>=0; J--) { |
for (J=M_beta_pt-1; J>=0; J--) { |
if (M_beta_1[J] == Nnu) { |
if (M_beta_1[J] == Nnu) { |
K=I+1; |
K=I+1; |
if (Debug) { |
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)); |
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. */ |
/* Check other conditions. See Numata's mail on Dec 24, 2011. */ |
rrMax = Nu[I]-1; |
rrMax = Nu[I]-1; |
if ((plength_t(Kapt) < rrMax) || (plength_t(Nut) < rrMax)) { |
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; |
break; |
} |
} |
OK=1; |
OK=1; |
for (RR=0; RR<rrMax; RR++) { |
for (RR=0; RR<rrMax; RR++) { |
if (Kapt[RR] != Nut[RR]) { OK=0; break;} |
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. */ |
/* check done. */ |
M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K); |
M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K); |
Done = 1; break; |
Done = 1; break; |
} |
} |
} |
} |
if (Done) break; else Nu[I]--; |
if (Done) break; else Nu[I]--; |
} |
} |
if (!Done) { |
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] = NAN; error("Not found."); */ |
M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu); |
M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu); |
} |
} |
/* Fix the bug of mh.rr */ |
/* Fix the bug of mh.rr */ |
M_beta_pt++; |
M_beta_pt++; |
} |
} |
static genBeta(int Kap[]) { |
static int genBeta(int Kap[]) { |
extern double *M_beta_0; |
extern double *M_beta_0; |
extern int *M_beta_1; |
extern int *M_beta_1; |
extern int M_beta_pt; |
extern int M_beta_pt; |
extern int M_beta_kap[]; |
extern int M_beta_kap[]; |
extern int P_pmn; |
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 = newmat(2,P_pmn+1); */ |
M_beta_0 = (double *)mymalloc(sizeof(double)*(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; |
M_beta_pt = 0; |
for (I=0; I<=P_pmn; I++) {M_beta_0[I] = NAN; M_beta_1[I] = -1;} |
for (I=0; I<=P_pmn; I++) {M_beta_0[I] = NAN; M_beta_1[I] = -1;} |
N = plength(Kap); |
N = plength(Kap); |
HS_hsExec = hsExec_beta; |
HS_hsExec = hsExec_beta; |
copyP(M_beta_kap,Kap); |
copyP(M_beta_kap,Kap); |
pListHS(Kap,N); |
pListHS(Kap,N); return(0); |
} |
} |
/* |
/* |
genDarray2(4,3); |
genDarray2(4,3); |
genBeta([2,2,0]); |
genBeta([2,2,0]); |
genBeta([2,1,1]); |
genBeta([2,1,1]); |
*/ |
*/ |
|
#ifdef STANDALONE |
static checkBeta1() { |
static int checkBeta1() { |
int Kap[3] = {2,2,0}; |
int Kap[3] = {2,2,0}; |
int Kap2[3] = {2,1,0}; |
int Kap2[3] = {2,1,0}; |
int I; |
int I; |
Line 988 static checkBeta1() { |
|
Line 1169 static checkBeta1() { |
|
for (I=0; I<M_beta_pt; I++) { |
for (I=0; I<M_beta_pt; I++) { |
Mu = Parray[M_beta_1[I]]; |
Mu = Parray[M_beta_1[I]]; |
Beta_km = M_beta_0[I]; |
Beta_km = M_beta_0[I]; |
if (Debug) { |
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); |
genBeta(Kap2); |
for (I=0; I<M_beta_pt; I++) { |
for (I=0; I<M_beta_pt; I++) { |
Mu = Parray[M_beta_1[I]]; |
Mu = Parray[M_beta_1[I]]; |
Beta_km = M_beta_0[I]; |
Beta_km = M_beta_0[I]; |
if (Debug) { |
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() { |
def checkBeta2() { |
genDarray2(3,3); |
genDarray2(3,3); |
Kap = [2,1,0]; |
Kap = [2,1,0]; |
printf("Kap=%a\n",Kap); |
oxprintf("Kap=%a\n",Kap); |
genBeta(Kap); |
genBeta(Kap); |
for (I=0; I<M_beta_pt; I++) { |
for (I=0; I<M_beta_pt; I++) { |
Mu = Parray[M_beta[1][I]]; |
Mu = Parray[M_beta[1][I]]; |
Beta_km = M_beta[0][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)); |
} |
} |
} |
} |
*/ |
*/ |
|
|
/* main() { checkBeta1(); } */ |
/* main() { checkBeta1(); } */ |
Line 1040 static int psublen(int Kap[],int Mu[]) { |
|
Line 1222 static int psublen(int Kap[],int Mu[]) { |
|
|
|
|
|
/* Table of Jack polynomials |
/* Table of Jack polynomials |
Jack[1]* one variable. |
Jack[1]* one variable. |
Jack[2]* two variables. |
Jack[2]* two variables. |
... |
... |
Jack[M_n]* n variables. |
Jack[M_n]* n variables. |
Jack[P][J]* |
Jack[P][J]* |
D^J(P variables jack of p variables). Example. J=001 d_1, 010 d_2, 100 d_3 |
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 |
0<=J<=2^{M_n}-1 |
Jack[P][J][nk(Kappa)] Jack_Kappa, Kappa is a partition. |
Jack[P][J][nk(Kappa)] Jack_Kappa, Kappa is a partition. |
0<=nk(Kappa)<=pmn(M_m,M_n) |
0<=nk(Kappa)<=pmn(M_m,M_n) |
*/ |
*/ |
|
|
#define aM_jack(i,j,k) ((M_jack[i])[(j)*(Pmn+1)+(k)]) |
#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 double **M_jack; |
extern int M_2n; |
extern int M_2n; |
extern int P_pmn; |
extern int P_pmn; |
Line 1060 static genJack(int M,int N) { |
|
Line 1242 static genJack(int M,int N) { |
|
int Pmn,I,J,K,L,Nv,H,P; |
int Pmn,I,J,K,L,Nv,H,P; |
int *Kap,*Mu; |
int *Kap,*Mu; |
double Jack,Beta_km; |
double Jack,Beta_km; |
int Nk,JJ; |
int Nk,JJ, two_to_I; |
|
if (Debug) oxprintf("genJack(%d,%d)\n",M,N); |
M_jack = (double **) mymalloc(sizeof(double *)*(N+2)); |
M_jack = (double **) mymalloc(sizeof(double *)*(N+2)); |
M_2n = imypower(2,N); |
M_2n = imypower(2,N); |
Pmn = pmn(M,N); /*P_pmn is initializeded. |
Pmn = pmn(M,N); /*P_pmn is initializeded. |
Line 1083 static genJack(int M,int N) { |
|
Line 1266 static genJack(int M,int N) { |
|
for (J=2; J<M_2n; J++) aM_jack(1,J,K) = 0; |
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++) { |
for (K=M+1; K<Pmn+1; K++) { |
aM_jack(I,0,K) = NAN; |
aM_jack(I,0,K) = NAN; |
if (M_df) { |
if (M_df) { |
for (J=1; J<M_2n; J++) { |
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; |
else aM_jack(I,J,K) = NAN; |
} |
} |
} |
} |
} |
} |
Line 1105 static genJack(int M,int N) { |
|
Line 1288 static genJack(int M,int N) { |
|
for (J=1; J<M_2n; J++) aM_jack(I,J,K) = 0; |
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 */ |
/* Enumerate horizontal strip of Kappa */ |
genBeta(Kap); /* M_beta_pt stores the number of hs */ |
genBeta(Kap); /* M_beta_pt stores the number of hs */ |
/* Nv is the number of variables */ |
/* Nv is the number of variables */ |
Line 1114 static genJack(int M,int N) { |
|
Line 1297 static genJack(int M,int N) { |
|
for (H=0; H<M_beta_pt; H++) { |
for (H=0; H<M_beta_pt; H++) { |
Nk = M_beta_1[H]; |
Nk = M_beta_1[H]; |
Mu = Parray[Nk]; |
Mu = Parray[Nk]; |
if (UseTable) { |
if (UseTable) { |
Beta_km = M_beta_0[H]; |
Beta_km = M_beta_0[H]; |
}else{ |
}else{ |
Beta_km = beta(Kap,Mu); |
Beta_km = beta(Kap,Mu); |
/* do not use the M_beta table. It's buggy. UseTable is experimental.*/ |
/* 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); |
P = psublen(Kap,Mu); |
Jack += aM_jack(Nv-1,0,Nk)*Beta_km*xval(Nv,P); /* util_v(x,[Nv])^P;*/ |
Jack += aM_jack(Nv-1,0,Nk)*Beta_km*xval(Nv,P); /* util_v(x,[Nv])^P;*/ |
|
if (Debug) oxprintf("xval(%d,%d)=%lf\n",Nv,P,xval(Nv,P)); |
} |
} |
aM_jack(Nv,0,K) = Jack; |
aM_jack(Nv,0,K) = Jack; |
if (M_df) { |
if (M_df) { |
/* The case of M_df > 0. */ |
/* The case of M_df > 0. */ |
for (J=1; J<M_2n; J++) { |
for (J=1; J<M_2n; J++) { |
Jack = 0; |
mh_check_intr(100); |
for (H=0; H<M_beta_pt; H++) { |
Jack = 0; |
Nk = M_beta_1[H]; |
for (H=0; H<M_beta_pt; H++) { |
Mu = Parray[Nk]; |
Nk = M_beta_1[H]; |
if (UseTable) { |
Mu = Parray[Nk]; |
Beta_km = M_beta_0[H]; |
if (UseTable) { |
}else{ |
Beta_km = M_beta_0[H]; |
Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */ |
}else{ |
} |
Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */ |
if (Debug) {printf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km); |
} |
printp(Mu); printf("\n"); } |
if (Debug) {oxprintf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km); |
P = psublen(Kap,Mu); |
printp(Mu); oxprintf("\n"); } |
if (J & (1 << (Nv-1))) { |
P = psublen(Kap,Mu); |
JJ = J & ((1 << (Nv-1)) ^ 0xffff); /* NOTE!! Up to 16 bits. mh-15 */ |
if (J & (1 << (Nv-1))) { |
if (P != 0) { |
JJ = J & ((1 << (Nv-1)) ^ 0xffff); /* NOTE!! Up to 16 bits. mh-15 */ |
Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1); |
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); |
}else{ |
} |
Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P); |
} |
} |
aM_jack(Nv,J,K) = Jack; |
} |
|
aM_jack(Nv,J,K) = Jack; |
|
if (Debug) oxprintf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack); |
} /* end of J loop */ |
} /* end of J loop */ |
} |
} |
} |
} |
} |
} return(0); |
} |
} |
|
|
|
#ifdef STANDALONE |
/* checkJack1(3,3) |
/* checkJack1(3,3) |
*/ |
*/ |
static checkJack1(int M,int N) { |
static int checkJack1(int M,int N) { |
int I,K; |
int I,K; |
extern int P_pmn; |
extern int P_pmn; |
extern double M_x[]; |
extern double M_x[]; |
Line 1174 static checkJack1(int M,int N) { |
|
Line 1360 static checkJack1(int M,int N) { |
|
for (I=1; I<=N; I++) { |
for (I=1; I<=N; I++) { |
for (K=0; K<=P_pmn; K++) { |
for (K=0; K<=P_pmn; K++) { |
printp(Parray[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); }*/ |
/*main() { checkJack1(3,3); }*/ |
|
|
|
|
static checkJack2(int M,int N) { |
static int checkJack2(int M,int N) { |
int I,K,J; |
int I,K,J; |
extern int P_pmn; |
extern int P_pmn; |
extern double M_x[]; |
extern double M_x[]; |
extern M_df; |
extern int M_df; |
int Pmn; /* used in aM_jack */ |
int Pmn; /* used in aM_jack */ |
M_df=1; |
M_df=1; |
/* initialize x vars. */ |
/* initialize x vars. */ |
Line 1199 static checkJack2(int M,int N) { |
|
Line 1386 static checkJack2(int M,int N) { |
|
for (I=1; I<=N; I++) { |
for (I=1; I<=N; I++) { |
for (K=0; K<=P_pmn; K++) { |
for (K=0; K<=P_pmn; K++) { |
printp(Parray[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 (I=1; I<=N; I++) { |
for (K=0; K<=P_pmn; K++) { |
for (K=0; K<=P_pmn; K++) { |
for (J=0; J<M_2n; J++) { |
for (J=0; J<M_2n; J++) { |
printp(Parray[K]); |
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)); |
I,J,aM_jack(I,J,K)); |
} |
} |
} |
} |
} |
} |
|
return(0); |
} |
} |
|
|
/* main() { checkJack2(3,3); } */ |
/* 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; |
double F,F2; |
extern int M_df; |
extern int M_df; |
extern int P_pmn; |
extern int P_pmn; |
extern double *M_qk; |
extern double *M_qk; |
extern double M_rel_error; |
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 Pmn; |
int K; |
int K; |
int *Kap; |
int *Kap; |
int size; |
int size; |
|
int i; |
|
double partial_sum[M_m_MAX+1]; |
|
double iv; |
|
double serror; |
F = 0; F2=0; |
F = 0; F2=0; |
M_df=1; |
M_df=1; |
genJack(M,N); |
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; |
Pmn = P_pmn; |
size = ParraySize[P_pmn]; |
size = ParraySize[P_pmn]; |
for (K=0; K<=P_pmn; K++) { |
for (K=0; K<=P_pmn; K++) { |
Kap = Parray[K]; |
mh_check_intr(100); |
M_qk[K] = qk(Kap,A,B); |
Kap = Parray[K]; |
F += M_qk[K]*aM_jack(N,0,K); |
M_qk[K] = qk(Kap,A,B); |
if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K); |
F += M_qk[K]*aM_jack(N,0,K); |
if (Debug) printf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size); |
if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K); |
if (Debug && (ParraySize[K] == size)) printf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K)); |
if (Debug) oxprintf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size); |
|
if (Debug && (ParraySize[K] == size)) oxprintf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K)); |
} |
} |
M_rel_error = F-F2; |
M_rel_error = F-F2; |
|
|
|
M_m_estimated_approx_deg = -1; serror=1; |
|
for (i=0; i<=M_m; i++) { |
|
partial_sum[i] = 0.0; partial_sum[i+1] = 0.0; |
|
for (K=0; K<=P_pmn; K++) { |
|
if (ParraySize[K] == i) partial_sum[i] += M_qk[K]*aM_jack(N,0,K); |
|
} |
|
if (i>0) partial_sum[i] += partial_sum[i-1]; |
|
if (i>0) serror = myabs((partial_sum[i]-partial_sum[i-1])/partial_sum[i-1]); |
|
if ((i>0)&&(M_m_estimated_approx_deg < 0)&&(serror<M_assigned_series_error)) { |
|
M_m_estimated_approx_deg = i; break; |
|
} |
|
} |
|
if (M_m_estimated_approx_deg < 0) { |
|
M_m_estimated_approx_deg = M_m+mymin(5,mymax(1,(int)log(serror/M_assigned_series_error))); /* Heuristic */ |
|
} |
|
/* |
|
for (K=0; K<=P_pmn; K++) { |
|
oxprintf("Kappa="); for (i=0; i<N; i++) oxprintf("%d ",Parray[K][i]); oxprintf("\n"); |
|
oxprintf("ParraySize(%d)=%d (|kappa|), M_m=%d\n",K,ParraySize[K],M_m); |
|
} |
|
for (i=0; i<=M_m; i++) { |
|
oxprintf("partial_sum[%d]=%lg\n",i,partial_sum[i]); |
|
} |
|
*/ |
|
M_estimated_X0g = X0g; |
|
iv=myabs(F*iv_factor()); |
|
if (iv < M_x0value_min) M_estimated_X0g = X0g*mymax(2,log(log(1/iv))); /* This is heuristic */ |
|
M_estimated_X0g = mymin(M_estimated_X0g,M_X0g_bound); |
|
M_mh_t_success = 1; |
|
if (M_estimated_X0g != X0g) M_mh_t_success=0; |
|
if (M_m_estimated_approx_deg > M_m) M_mh_t_success=0; |
|
|
|
M_series_error = serror; |
|
M_recommended_abserr = iv*M_assigned_series_error; |
|
M_recommended_relerr = M_series_error; |
|
|
|
if (M_show_autosteps) { |
|
oxprintf("%%%%serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m); |
|
oxprintf("%%%%x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound); |
|
oxprintf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g); |
|
oxprintfe("%%%%(stderr) serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m); |
|
oxprintfe("%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound); |
|
oxprintfe("%%%%(stderr) F=%lg,Ef=%lg,iv=|F*Ef|=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),iv,M_estimated_X0g,X0g); |
|
} |
|
|
|
M_mh_t_value=F; |
return(F); |
return(F); |
} |
} |
double mh_t2(int J) { |
double mh_t2(int J) { |
Line 1255 double mh_t2(int J) { |
|
Line 1499 double mh_t2(int J) { |
|
F = 0; |
F = 0; |
Pmn = P_pmn; |
Pmn = P_pmn; |
for (K=0; K<P_pmn; K++) { |
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); |
return(F); |
} |
} |
|
|
static mtest1b() { |
#ifdef STANDALONE |
|
static int mtest1b() { |
double A[1] = {1.5}; |
double A[1] = {1.5}; |
double B[1] = {1.5+5}; |
double B[1] = {1.5+5}; |
int I,N,M,J; |
int I,N,M,J; |
Line 1271 static mtest1b() { |
|
Line 1516 static mtest1b() { |
|
} |
} |
mh_t(A,B,N,M); |
mh_t(A,B,N,M); |
for (J=0; J<M_2n; J++) { |
for (J=0; J<M_2n; J++) { |
F=mh_t2(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(); }*/ |
/* main() { mtest1b(); }*/ |
|
#endif |
|
|
|
|
|
|
|
|
#define TEST 1 |
#define TEST 1 |
#ifndef TEST |
#ifndef TEST |
|
|
Line 1289 static mtest1b() { |
|
Line 1535 static mtest1b() { |
|
/****** from mh-n.c *****/ |
/****** from mh-n.c *****/ |
|
|
#define SMAX 4096 |
#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) { |
static int imypower(int x,int n) { |
int i; |
int i; |
Line 1300 static int imypower(int x,int n) { |
|
Line 1546 static int imypower(int x,int n) { |
|
return(v); |
return(v); |
} |
} |
|
|
#ifdef STANDALONE |
#ifdef STANDALONE2 |
main(int argc,char *argv[]) { |
int main(int argc,char *argv[]) { |
mh_exit(MH_RESET_EXIT); |
mh_exit(MH_RESET_EXIT); |
/* jk_main(argc,argv); |
/* jk_main(argc,argv); |
printf("second run.\n"); */ |
oxprintf("second run.\n"); */ |
jk_main(argc,argv); |
jk_main(argc,argv); |
|
return(0); |
} |
} |
#endif |
#endif |
|
|
struct MH_RESULT *jk_main(int argc,char *argv[]) { |
struct MH_RESULT *jk_main(int argc,char *argv[]) { |
double *y0; |
int i; |
double x0,xn; |
struct MH_RESULT *ans; |
double ef; |
extern int M_automatic; |
int i,j,rank; |
extern int M_mh_t_success; |
double a[1]; double b[1]; |
extern double M_estimated_X0g; |
|
extern int M_m_estimated_approx_deg; |
|
for (i=1; i<argc; i++) { |
|
if (strcmp(argv[i],"--automatic")==0) { |
|
inci(i); |
|
sscanf(argv[i],"%d",&M_automatic); |
|
break; |
|
} |
|
} |
|
ans=jk_main2(argc,argv,0,0.0,0); |
|
if (!M_automatic) return(ans); |
|
if (M_mh_t_success) return(ans); |
|
while (!M_mh_t_success) { |
|
ans=jk_main2(argc,argv,1,M_estimated_X0g,M_m_estimated_approx_deg); |
|
} |
|
return(ans); |
|
} |
|
|
|
struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree) { |
|
// double *y0; |
|
// double x0,xn; |
|
// double ef; |
|
|
|
int i,j; // int i,j,rank; |
extern double M_x[]; |
extern double M_x[]; |
extern double *Beta; |
extern double *Beta; |
extern int M_2n; |
extern int M_2n; |
|
extern int M_mh_t_success; |
char swork[1024]; |
char swork[1024]; |
struct MH_RESULT *ans=NULL; |
struct MH_RESULT *ans=NULL; |
struct SFILE *ofp = NULL; |
struct SFILE *ofp = NULL; |
int idata=0; |
int idata=0; |
JK_byFile = 1; |
JK_byFile = 1; |
jk_initializeWorkArea(); |
jk_initializeWorkArea(); |
Debug = 0; |
|
UseTable = 1; |
UseTable = 1; |
Mapprox=6; |
Mapprox=6; |
for (i=1; i<argc; i++) { |
for (i=1; i<argc; i++) { |
if (strcmp(argv[i],"--idata")==0) { |
if (strcmp(argv[i],"--idata")==0) { |
inci(i); |
inci(i); |
setParam(argv[i]); idata=1; |
setParam(argv[i]); idata=1; |
}else if (strcmp(argv[i],"--degree")==0) { |
}else if (strcmp(argv[i],"--degree")==0) { |
inci(i); |
inci(i); |
sscanf(argv[i],"%d",&Mapprox); |
sscanf(argv[i],"%d",&Mapprox); |
}else if (strcmp(argv[i],"--x0")==0) { |
}else if (strcmp(argv[i],"--x0")==0) { |
inci(i); |
inci(i); |
sscanf(argv[i],"%lg",&X0g); |
sscanf(argv[i],"%lg",&X0g); |
}else if (strcmp(argv[i],"--notable")==0) { |
}else if (strcmp(argv[i],"--notable")==0) { |
UseTable = 0; |
UseTable = 0; |
}else if (strcmp(argv[i],"--debug")==0) { |
}else if (strcmp(argv[i],"--debug")==0) { |
Debug = 1; |
Debug = 1; |
}else if (strcmp(argv[i],"--help")==0) { |
}else if (strcmp(argv[i],"--help")==0) { |
usage(); return(0); |
usage(); return(0); |
}else if (strcmp(argv[i],"--bystring")==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; |
JK_byFile = 0; |
}else { |
}else if (strcmp(argv[i],"--automatic")==0) { |
fprintf(stderr,"Unknown option %s\n",argv[i]); |
inci(i); /* ignore, in this function */ |
usage(); |
}else if (strcmp(argv[i],"--assigned_series_error")==0) { |
return(NULL); |
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 { |
|
oxprintfe("Unknown option %s\n",argv[i]); |
|
usage(); |
|
return(NULL); |
|
} |
} |
} |
if (!idata) setParam(NULL); |
if (!idata) setParam(NULL); |
|
if (automode) { |
|
Mapprox = newDegree; |
|
X0g = newX0g; |
|
} |
|
|
/* Initialize global variables */ |
/* Initialize global variables */ |
M_n = Mg; |
M_n = Mg; |
HS_n=M_n; |
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; |
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 */ |
/* Output by a file=stdout */ |
ofp = mh_fopen("stdout","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,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp); |
sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp); |
sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp); |
if (M_n != Mg) { |
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); |
|
setM_x(); |
|
|
|
M_beta_i_x_o2_max=myabs(M_x[0]/2); |
|
if (M_n <= 1) M_beta_i_beta_j_min = myabs(Beta[0]); |
|
else M_beta_i_beta_j_min = myabs(Beta[1]-Beta[0]); |
for (i=0; i<M_n; i++) { |
for (i=0; i<M_n; i++) { |
M_x[i] = Beta[i]*X0g; |
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++) { |
a[0] = ((double)Mg+1.0)/2.0; |
if (myabs(Beta[i]-Beta[j]) < M_beta_i_beta_j_min) |
b[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */ |
M_beta_i_beta_j_min = myabs(Beta[i]-Beta[j]); |
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); |
} |
|
|
|
mh_t(A_pFq,B_pFq,M_n,Mapprox); |
|
if ((!M_mh_t_success) && M_automatic) { |
|
jk_freeWorkArea(); |
|
return NULL; |
|
} |
if (imypower(2,M_n) != M_2n) { |
if (imypower(2,M_n) != M_2n) { |
sprintf(swork,"M_n=%d,M_2n=%d\n",M_n,M_2n); mh_fputs(swork,ofp); |
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); |
myerror("2^M_n != M_2n\n"); mh_exit(-1); |
} |
} |
sprintf(swork,"%%%%M_rel_error=%lg\n",M_rel_error); mh_fputs(swork,ofp); |
sprintf(swork,"%%%%M_rel_error=%lg\n",M_rel_error); mh_fputs(swork,ofp); |
for (j=0; j<M_2n; j++) { |
for (j=0; j<M_2n; j++) { |
Iv[j] = mh_t2(j); |
Iv[j] = mh_t2(j); |
} |
} |
Ef = iv_factor(); |
Ef = iv_factor(); |
showParam(ofp); |
showParam(ofp,0); |
|
|
/* return the result */ |
/* return the result */ |
if (!JK_byFile) { |
if (!JK_byFile) { |
ans->x = X0g; |
ans->x = X0g; |
ans->rank = imypower(2,Mg); |
ans->rank = imypower(2,Mg); |
ans->y = (double *)mymalloc(sizeof(double)*(ans->rank)); |
ans->y = (double *)mymalloc(sizeof(double)*(ans->rank)); |
for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i]; |
for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i]; |
ans->size=1; |
ans->size=1; |
ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size)); |
ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size)); |
(ans->sfpp)[0] = ofp; |
(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) oxprintf("jk_freeWorkArea() starts\n"); |
jk_freeWorkArea(); |
jk_freeWorkArea(); |
|
if (Debug) oxprintf("jk_freeWorkArea() has finished.\n"); |
return(ans); |
return(ans); |
} |
} |
|
|
static usage() { |
static int usage() { |
fprintf(stderr,"Usages:\n"); |
oxprintfe("Usages:\n"); |
fprintf(stderr,"mh-m [--idata input_data_file --x0 x0 --degree approxm]\n"); |
#include "usage-jack-n.h" |
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"); |
return(0); |
fprintf(stderr,"The mh-m 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," 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," The line started with %% is a comment line.\n"); |
|
fprintf(stderr," With the --notable option, it does not use the Lemma 3.2 of Koev-Edelman (there is a typo: kappa'_r = mu'_r for 1<=r<=mu_k).\n"); |
|
fprintf(stderr," An example format of the input_data_file can be obtained by executing mh-m with no option.\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," gnuplot -persist <test-g-gp.txt\n"); |
|
} |
} |
|
|
static setParamDefault() { |
static int setParamDefault() { |
int rank; |
int rank; |
int i; |
int i; |
Mg = M_n_default ; |
Mg = M_n_default ; |
Line 1450 static setParamDefault() { |
|
Line 1738 static setParamDefault() { |
|
Hg = 0.001; |
Hg = 0.001; |
Dp = 1; |
Dp = 1; |
Xng = 10.0; |
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 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] = '%'; |
s[0] = '%'; |
while (s[0] == '%') { |
while ((s[0] == '%') || (s[0] == '#')) { |
if (!mh_fgets(s,SMAX,sfp)) { |
if (!mh_fgets(s,SMAX,sfp)) { |
fprintf(stderr,"Data format error at %s\n",msg); |
oxprintfe("Data format error at %s\n",msg); |
mh_exit(-1); |
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 (s[0] != '%') return(0); |
if (check && (strncmp(msg,ng,4)==0)) { |
|
if (strncmp(s,ng,5) != 0) { |
|
oxprintfe("Warning, there is no %%Ng= at the border of Beta's and Ng, s=%s\n",s); |
|
} |
|
/* check=0; */ |
|
} |
|
if (s[0] != '%') return(0); |
} |
} |
|
return(0); |
} |
} |
static setParam(char *fname) { |
static int setParam(char *fname) { |
int rank; |
int rank; |
char s[SMAX]; |
char s[SMAX]; |
struct SFILE *fp; |
struct SFILE *fp; |
int i; |
int i; |
|
struct mh_token tk; |
|
int version; |
if (fname == NULL) return(setParamDefault()); |
if (fname == NULL) return(setParamDefault()); |
|
|
Sample = 0; |
Sample = 0; |
if ((fp=mh_fopen(fname,"r",JK_byFile)) == NULL) { |
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); |
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); |
sscanf(s,"%d",&Mg); |
rank = imypower(2,Mg); |
rank = imypower(2,Mg); |
|
|
Beta = (double *)mymalloc(sizeof(double)*Mg); |
Beta = (double *)mymalloc(sizeof(double)*Mg); |
for (i=0; i<Mg; i++) { |
for (i=0; i<Mg; i++) { |
next(fp,s,"Beta"); |
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); |
sscanf(s,"%lf",Ng); |
|
|
next(fp,s,"X0g(initial point)"); |
next(fp,s,"X0g(initial point)"); |
sscanf(s,"%lf",&X0g); |
sscanf(s,"%lf",&X0g); |
|
|
Iv = (double *)mymalloc(sizeof(double)*rank); |
Iv = (double *)mymalloc(sizeof(double)*rank); |
for (i=0; i<rank; i++) { |
for (i=0; i<rank; i++) { |
next(fp,s,"Iv(initial values)"); |
next(fp,s,"Iv(initial values)"); |
sscanf(s,"%lg",&(Iv[i])); |
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)"); |
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)"); |
next(fp,s,"Hg (step size of rk)"); |
sscanf(s,"%lf",&Hg); |
sscanf(s,"%lf",&Hg); |
Line 1508 static setParam(char *fname) { |
|
Line 1834 static setParam(char *fname) { |
|
|
|
next(fp,s,"Xng (the last point, cf. --xmax)"); |
next(fp,s,"Xng (the last point, cf. --xmax)"); |
sscanf(s,"%lf",&Xng); |
sscanf(s,"%lf",&Xng); |
|
|
|
/* Reading the optional parameters */ |
|
myparse: |
|
while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) { |
|
/* expect ID */ |
|
if (tk.type != MH_TOKEN_ID) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if (strcmp(s,"automatic")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
M_automatic = tk.ival; |
|
continue; |
|
} |
|
if (strcmp(s,"assigned_series_error")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
M_assigned_series_error = tk.dval; |
|
continue; |
|
} |
|
if (strcmp(s,"x0value_min")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
M_x0value_min = tk.dval; |
|
continue; |
|
} |
|
if ((strcmp(s,"Mapprox")==0) || (strcmp(s,"degree")==0)) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
Mapprox = tk.ival; |
|
continue; |
|
} |
|
if (strcmp(s,"X0g_bound")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
M_X0g_bound = tk.dval; |
|
continue; |
|
} |
|
if (strcmp(s,"show_autosteps")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
M_show_autosteps = tk.ival; |
|
continue; |
|
} |
|
// Format: #p_pFq=2 1.5 3.2 |
|
if (strcmp(s,"p_pFq")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
setA(NULL,tk.ival); |
|
for (i=0; i<A_LEN; i++) { |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
A_pFq[i] = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
A_pFq[i] = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"q_pFq")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
setB(NULL,tk.ival); |
|
for (i=0; i<B_LEN; i++) { |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
B_pFq[i] = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
B_pFq[i] = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"ef_type")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
Ef_type = tk.ival; |
|
continue; |
|
} |
|
|
|
if (strcmp(s,"Mg")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
Mg = tk.ival; |
|
rank = imypower(2,Mg); |
|
continue; |
|
} |
|
if (strcmp(s,"Beta")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if (Mg <= 0) { |
|
oxprintfe("Mg should be set before reading Beta.\n"); mh_exit(-1); |
|
} |
|
Beta = (double *)mymalloc(sizeof(double)*Mg); |
|
for (i=0; i<Mg; i++) { |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
Beta[i] = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
Beta[i] = tk.ival; |
|
}else { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
} |
|
Iv = (double *)mymalloc(sizeof(double)*rank); |
|
for (i=0; i<rank; i++) { |
|
Iv[i] = 0.0; |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"Ng")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
|
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type== MH_TOKEN_DOUBLE) { |
|
*Ng = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
*Ng = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"X0g")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
X0g = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
X0g = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"Iv")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
for (i=0; i<rank; i++) { |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
Iv[i] = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
Iv[i] = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"Ef")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
Ef = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
Ef = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"Hg")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
Hg = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
Hg = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
continue; |
|
} |
|
if (strcmp(s,"Dp")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
Dp = tk.dval; |
|
continue; |
|
} |
|
if (strcmp(s,"Xng")==0) { |
|
if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) { |
|
Xng = tk.dval; |
|
}else if (tk.type == MH_TOKEN_INT) { |
|
Xng = tk.ival; |
|
}else{ |
|
oxprintfe("Syntax error at %s\n",s); mh_exit(-1); |
|
} |
|
continue; |
|
} |
|
|
|
oxprintfe("Unknown ID at %s\n",s); mh_exit(-1); |
|
} |
|
|
|
/* 1F1, original wishart case. */ |
|
if ((A_LEN <= 1) && (B_LEN <= 1) && (Ef_type==1) && (*Ng >= 0)) { |
|
if (A_LEN <1) setA(NULL,1); |
|
if (B_LEN <1) setB(NULL,1); |
|
A_pFq[0] = ((double)Mg+1.0)/2.0; |
|
B_pFq[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */ |
|
if (Debug) oxprintf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",A_pFq[0],B_pFq[0],M_n,Mapprox); |
|
} |
mh_fclose(fp); |
mh_fclose(fp); |
|
return(0); |
} |
} |
|
#ifdef STANDALONE |
static showParam(struct SFILE *fp) { |
/* may remove */ |
|
static int showParam_v1(struct SFILE *fp,int fd) { |
int rank,i; |
int rank,i; |
char swork[1024]; |
char swork[1024]; |
|
if (fd) { |
|
fp = mh_fopen("stdout","w",1); |
|
} |
rank = imypower(2,Mg); |
rank = imypower(2,Mg); |
sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp); |
sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp); |
for (i=0; i<Mg; i++) { |
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); |
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,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp); |
for (i=0; i<rank; i++) { |
for (i=0; i<rank; i++) { |
sprintf(swork,"%%Iv[%d]=\n%lg\n",i,Iv[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)) { |
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]-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,"%%Ef=\n%lg\n",Ef); mh_fputs(swork,fp); |
sprintf(swork,"%%Hg=\n%lf\n",Hg); 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,"%%Dp=\n%d\n",Dp); mh_fputs(swork,fp); |
sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp); |
sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp); |
|
|
|
sprintf(swork,"%%%% Optional paramters\n"); mh_fputs(swork,fp); |
|
sprintf(swork,"#success=%d\n",M_mh_t_success); mh_fputs(swork,fp); |
|
sprintf(swork,"#automatic=%d\n",M_automatic); mh_fputs(swork,fp); |
|
sprintf(swork,"#series_error=%lg\n",M_series_error); mh_fputs(swork,fp); |
|
sprintf(swork,"#recommended_abserr\n"); mh_fputs(swork,fp); |
|
sprintf(swork,"%%abserror=%lg\n",M_recommended_abserr); mh_fputs(swork,fp); |
|
if (M_recommended_relerr < MH_RELERR_DEFAULT) { |
|
sprintf(swork,"%%relerror=%lg\n",M_recommended_relerr); mh_fputs(swork,fp); |
|
} |
|
sprintf(swork,"#mh_t_value=%lg # Value of matrix hg at X0g.\n",M_mh_t_value); mh_fputs(swork,fp); |
|
sprintf(swork,"# M_m=%d # Approximation degree of matrix hg.\n",M_m); mh_fputs(swork,fp); |
|
sprintf(swork,"#beta_i_x_o2_max=%lg #max(|beta[i]*x|/2)\n",M_beta_i_x_o2_max); mh_fputs(swork,fp); |
|
sprintf(swork,"#beta_i_beta_j_min=%lg #min(|beta[i]-beta[j]|)\n",M_beta_i_beta_j_min); mh_fputs(swork,fp); |
|
sprintf(swork,"# change # to %% to read as an optional parameter.\n"); mh_fputs(swork,fp); |
|
sprintf(swork,"%%p_pFq=%d, ",A_LEN); mh_fputs(swork,fp); |
|
for (i=0; i<A_LEN; i++) { |
|
if (i != A_LEN-1) sprintf(swork," %lg,",A_pFq[i]); |
|
else sprintf(swork," %lg\n",A_pFq[i]); |
|
mh_fputs(swork,fp); |
|
} |
|
sprintf(swork,"%%q_pFq=%d, ",B_LEN); mh_fputs(swork,fp); |
|
for (i=0; i<B_LEN; i++) { |
|
if (i != B_LEN-1) sprintf(swork," %lg,",B_pFq[i]); |
|
else sprintf(swork," %lg\n",B_pFq[i]); |
|
mh_fputs(swork,fp); |
|
} |
|
sprintf(swork,"%%ef_type=%d\n",Ef_type); mh_fputs(swork,fp); |
|
return(0); |
} |
} |
|
#endif |
|
/* version2.0 format */ |
|
static int showParam(struct SFILE *fp,int fd) { |
|
int rank,i; |
|
char swork[1024]; |
|
if (fd) { |
|
fp = mh_fopen("stdout","w",1); |
|
} |
|
rank = imypower(2,Mg); |
|
sprintf(swork,"%s\n",VSTRING); mh_fputs(swork,fp); |
|
sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp); |
|
sprintf(swork,"%%p_pFq=%d, ",A_LEN); mh_fputs(swork,fp); |
|
for (i=0; i<A_LEN; i++) { |
|
if (i != A_LEN-1) sprintf(swork," %lg,",A_pFq[i]); |
|
else sprintf(swork," %lg\n",A_pFq[i]); |
|
mh_fputs(swork,fp); |
|
} |
|
sprintf(swork,"%%q_pFq=%d, ",B_LEN); mh_fputs(swork,fp); |
|
for (i=0; i<B_LEN; i++) { |
|
if (i != B_LEN-1) sprintf(swork," %lg,",B_pFq[i]); |
|
else sprintf(swork," %lg\n",B_pFq[i]); |
|
mh_fputs(swork,fp); |
|
} |
|
sprintf(swork,"%%ef_type=%d\n",Ef_type); mh_fputs(swork,fp); |
|
sprintf(swork,"%%Beta=\n"); mh_fputs(swork,fp); |
|
for (i=0; i<Mg; i++) { |
|
sprintf(swork,"#Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp); |
|
} |
|
if (*Ng >= 0) { |
|
sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp); |
|
} |
|
sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp); |
|
sprintf(swork,"%%Iv=\n"); mh_fputs(swork,fp); |
|
for (i=0; i<rank; i++) { |
|
sprintf(swork,"#Iv[%d]=\n%lg\n",i,Iv[i]); mh_fputs(swork,fp); |
|
if (Sample && (M_n == 2) && (X0g == 0.3)) { |
|
sprintf(swork,"%%Iv[%d]-Iv2[%d]=%lg\n",i,i,Iv[i]-Iv2[i]); mh_fputs(swork,fp); |
|
} |
|
} |
|
sprintf(swork,"%%Ef=\n%lg\n",Ef); mh_fputs(swork,fp); |
|
sprintf(swork,"%%Hg=\n%lf\n",Hg); mh_fputs(swork,fp); |
|
sprintf(swork,"%%Dp=\n%d\n",Dp); mh_fputs(swork,fp); |
|
sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp); |
|
|
|
sprintf(swork,"%%%% Optional paramters\n"); mh_fputs(swork,fp); |
|
sprintf(swork,"#success=%d\n",M_mh_t_success); mh_fputs(swork,fp); |
|
sprintf(swork,"#automatic=%d\n",M_automatic); mh_fputs(swork,fp); |
|
sprintf(swork,"#series_error=%lg\n",M_series_error); mh_fputs(swork,fp); |
|
sprintf(swork,"#recommended_abserr\n"); mh_fputs(swork,fp); |
|
sprintf(swork,"#abserror=%lg\n",M_recommended_abserr); mh_fputs(swork,fp); |
|
if (M_recommended_relerr < MH_RELERR_DEFAULT) { |
|
sprintf(swork,"%%relerror=%lg\n",M_recommended_relerr); mh_fputs(swork,fp); |
|
} |
|
sprintf(swork,"#mh_t_value=%lg # Value of matrix hg at X0g.\n",M_mh_t_value); mh_fputs(swork,fp); |
|
sprintf(swork,"# M_m=%d # Approximation degree of matrix hg.\n",M_m); mh_fputs(swork,fp); |
|
sprintf(swork,"#beta_i_x_o2_max=%lg #max(|beta[i]*x|/2)\n",M_beta_i_x_o2_max); mh_fputs(swork,fp); |
|
sprintf(swork,"#beta_i_beta_j_min=%lg #min(|beta[i]-beta[j]|)\n",M_beta_i_beta_j_min); mh_fputs(swork,fp); |
|
sprintf(swork,"# change # to %% to read as an optional parameter.\n"); mh_fputs(swork,fp); |
|
return(0); |
|
} |
|
|
static double gammam(double a,int n) { |
static double gammam(double a,int n) { |
double v,v2; |
double v,v2; |
int i; |
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; |
v2=0; |
for (i=1; i<=n; i++) { |
for (i=1; i<=n; i++) { |
v2 += lgamma(a-((double)(i-1))/2.0); /* not for big n */ |
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)); |
return(v*exp(v2)); |
} |
} |
|
|
static double iv_factor(void) { |
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 v1; |
double t; |
double t; |
double b; |
double b; |
Line 1562 static double iv_factor(void) { |
|
Line 2269 static double iv_factor(void) { |
|
detSigma = 1.0/(b*mypower(2.0,M_n)); |
detSigma = 1.0/(b*mypower(2.0,M_n)); |
|
|
c = gammam(((double)(M_n+1))/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); |
return( c*v1); |
} |
} |
|
|
|
static double iv_factor_ef_type2(void) { |
|
double ef; |
|
int i,m; |
|
double a,b,c; |
|
double t; |
|
a = A_pFq[0]; b = A_pFq[1]; c= B_pFq[0]; |
|
m = M_n; |
|
ef = 1.0; |
|
|
|
/*fprintf(stderr,"iv_factor_ef_type2: a=%lf, b=%lf, c=%lf, m=%d\n",a,b,c,m);*/ |
|
/* Ref: note 2016.02.04 */ |
|
if (X0g<0){ myerror("Negative X0g\n"); mh_exit(-1);} |
|
t = 0; |
|
for (i=0; i<m; i++) if (Beta[i]<0){ myerror("Negative beta\n"); mh_exit(-1);} |
|
for (i=0; i<m; i++) t += log(Beta[i]); |
|
ef *= exp((a+b-c)*t); |
|
|
|
t = 0; |
|
for (i=0; i<m; i++) t += log(Beta[i]+X0g); |
|
ef *= exp(-b*t); |
|
|
|
ef *= exp((c-a)*(2*a-1)*log(X0g)); |
|
|
|
ef *= gammam(b,m)/gammam(a+b-c,m); |
|
ef *= gammam(a,m)/gammam(c,m); |
|
return(ef); |
|
} |
|
static void setM_x(void) { |
|
if (Ef_type <= 1) return(setM_x_ef_type1()); |
|
else if (Ef_type==2) return(setM_x_ef_type2()); |
|
setM_x_ef_type1(); |
|
} |
|
static void setM_x_ef_type1(void) { |
|
int i; |
|
for (i=0; i<M_n; i++) { |
|
M_x[i] = Beta[i]*X0g; |
|
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); |
|
} |