Annotation of OpenXM/src/hgm/mh/src/jack-n.c, Revision 1.40
1.1 takayama 1: #include <stdio.h>
2: #include <stdlib.h>
3: #define _ISOC99_SOURCE
4: #include <math.h>
5: #include <string.h>
6: #include "sfile.h"
7: /*
1.40 ! takayama 8: $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.39 2016/02/04 02:52:19 takayama Exp $
1.12 takayama 9: Ref: copied from this11/misc-2011/A1/wishart/Prog
10: jack-n.c, translated from mh.rr or tk_jack.rr in the asir-contrib. License: LGPL
11: Koev-Edelman for higher order derivatives.
12: cf. 20-my-note-mh-jack-n.pdf, /Movies/oxvh/priv-j1/2011-12-14-ke-mh-jack.mov
13: Todo:
14: 1. Cash the transposes of partitions.
15: 2. Use the recurrence to obtain beta().
16: 3. Easier input data file format for mh-n.c
17: Changelog:
1.40 ! takayama 18: 2016.02.04 Ef_type (exponential or scalar factor type)
1.38 takayama 19: 2016.02.01 ifdef C_2F1 ...
1.36 takayama 20: 2016.01.12 2F1
1.22 takayama 21: 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.
1.18 takayama 22: 2014.03.14, --automatic option. Output estimation data.
1.12 takayama 23: 2012.02.21, porting to OpenXM/src/hgm
24: 2011.12.22, --table option, which is experimental.
25: 2011.12.19, bug fix. jjk() should return double. It can become more than max int.
26: 2011.12.15, mh.r --> jack-n.c
1.1 takayama 27: */
28:
1.3 takayama 29: /****** from mh-n.c *****/
30: static int JK_byFile=1;
1.6 takayama 31: static int JK_deallocate=0;
1.3 takayama 32: #define M_n_default 3
1.4 takayama 33: #define Sample_default 1
1.3 takayama 34: static int M_n=0;
1.1 takayama 35: /* global variables. They are set in setParam() */
1.2 takayama 36: static int Mg; /* n */
37: static int Mapprox; /* m, approximation degree */
38: static double *Beta; /* beta[0], ..., beta[m-1] */
39: static double *Ng; /* freedom n. c=(m+1)/2+n/2; Note that it is a pointer */
40: static double X0g; /* initial point */
41: static double *Iv; /* Initial values of mhg sorted by mhbase() in rd.rr at beta*x0 */
42: static double Ef; /* exponential factor at beta*x0 */
43: static double Hg; /* step size of rk defined in rk.c */
44: static int Dp; /* Data sampling period */
45: static double Xng=0.0; /* the last point */
1.4 takayama 46: static int Sample = Sample_default;
1.1 takayama 47:
48: /* for sample inputs */
1.2 takayama 49: static double *Iv2;
50: static double Ef2;
1.1 takayama 51:
52: #ifdef NAN
53: #else
1.2 takayama 54: #define NAN 3.40282e+38 /* for old 32 bit machines. Todo, configure */
1.1 takayama 55: #endif
56:
57: /* #define M_n 3 defined in the Makefile */ /* number of variables */
58: #define M_n0 3 /* used for tests. Must be equal to M_n */
59: #define M_m_MAX 200
1.3 takayama 60: #define M_nmx M_m_MAX /* maximal of M_n */
1.37 takayama 61: #ifdef C_2F1
62: #define A_LEN 2 /* (a_1) , (a_1, ..., a_p)*/
63: #define B_LEN 1 /* (b_1) */
64: static int P_pFq=2;
65: static int Q_pFq=1;
66: #else
1.1 takayama 67: #define A_LEN 1 /* (a_1) , (a_1, ..., a_p)*/
68: #define B_LEN 1 /* (b_1) */
1.36 takayama 69: static int P_pFq=1;
70: static int Q_pFq=1;
1.37 takayama 71: #endif
1.36 takayama 72: static double A_pFq[A_LEN];
73: static double B_pFq[B_LEN];
1.38 takayama 74: #ifndef C_2F1
1.36 takayama 75: static int Orig_1F1=1;
1.40 ! takayama 76: static int Ef_type=1; /* 1F1 for Wishart */
1.38 takayama 77: #else
78: static int Orig_1F1=0;
1.40 ! takayama 79: static int Ef_type=2; /* 2F1, Hashiguchi note (1) */
1.38 takayama 80: #endif
1.9 takayama 81: static int Debug = 0;
1.5 takayama 82: static int Alpha = 2; /* 2 implies the zonal polynomial */
1.2 takayama 83: static int *Darray = NULL;
84: static int **Parray = NULL; /* array of partitions of size M_n */
85: static int *ParraySize = NULL; /* length of each partitions */
1.3 takayama 86: static int M_kap[M_nmx];
1.2 takayama 87: static int M_m=M_m_MAX-2; /* | | <= M_m, bug do check of M_m <=M_m_MAX-2 */
1.1 takayama 88: void (*M_pExec)(void);
1.3 takayama 89: static int HS_mu[M_nmx];
90: static int HS_n=M_nmx; /* It is initialized to M_n in jk_main */
1.1 takayama 91: void (*HS_hsExec)(void);
1.3 takayama 92: static double M_x[M_nmx];
1.1 takayama 93:
94: /* They are used in pmn */
1.2 takayama 95: static int *P_pki=NULL;
96: static int P_pmn=0;
1.1 takayama 97:
98: /* It is used genDarray2(), list partitions... */
1.2 takayama 99: static int DR_parray=0;
1.1 takayama 100:
101: /* Used in genBeta() and enumeration of horizontal strip. */
1.2 takayama 102: static double *M_beta_0=NULL; /* M_beta[0][*] value of beta_{kappa,mu}, M_beta[1][*] N_mu */
103: static int *M_beta_1=NULL;
104: static int M_beta_pt=0;
1.3 takayama 105: static int M_beta_kap[M_nmx];
1.2 takayama 106: static int UseTable = 0;
107:
108: static double **M_jack;
109: static int M_df=1; /* Compute differentials? */
110: static int M_2n=0; /* 2^N */
1.1 takayama 111:
1.3 takayama 112: static double Xarray[M_nmx][M_m_MAX];
1.1 takayama 113: /* (x_1, ..., x_n) */
114: /* Xarray[i][0] x_{i+1}^0, Xarray[i][1], x_{i+1}^1, ... */
115:
1.2 takayama 116: static double *M_qk=NULL; /* saves pochhammerb */
117: static double M_rel_error=0.0; /* relative errors */
1.1 takayama 118:
1.16 takayama 119: /* For automatic degree and X0g setting. */
1.19 takayama 120: /*
121: If automatic == 1, then the series is reevaluated as long as t_success!=1
122: by increasing X0g (evaluation point) and M_m (approx degree);
123: */
1.26 takayama 124: static int M_automatic=1;
1.19 takayama 125: /* Estimated degree bound for series expansion. See mh_t */
1.16 takayama 126: static int M_m_estimated_approx_deg=0;
1.19 takayama 127: /* Let F(i) be the approximation up to degree i.
128: The i-th series_error is defined
129: by |(F(i)-F(i-1))/F(i-1)|.
130: */
131: static double M_series_error;
132: /*
133: M_series_error < M_assigend_series_error (A) is required for the
134: estimated_approx_deg.
135: */
1.26 takayama 136: static double M_assigned_series_error=M_ASSIGNED_SERIES_ERROR_DEFAULT;
1.19 takayama 137: /*
138: Let Ef be the exponential factor ( Ef=(4)/1F1 of [HNTT] )
139: If F(M_m)*Ef < x0value_min (B), the success=0 and X0g is increased.
140: Note that minimal double is about 2e-308
141: */
1.26 takayama 142: static double M_x0value_min=1e-60;
1.19 takayama 143: /*
144: estimated_X0g is the suggested value of X0g.
145: */
1.16 takayama 146: static double M_estimated_X0g=0.0;
1.19 takayama 147: /*
148: X0g should be less than M_X0g_bound.
149: */
150: static double M_X0g_bound = 1e+100;
151: /*
152: success is set to 1 when (A) and (B) are satisfied.
153: */
1.16 takayama 154: static int M_mh_t_success=1;
1.19 takayama 155: /*
156: recommended_abserr is the recommended value of the absolute error
157: for the Runge-Kutta method. It is defined as
158: assigend_series_error(standing for significant digits)*Ig[0](initial value)
159: */
160: static double M_recommended_abserr;
161: /*
1.27 takayama 162: recommended_relerr is the recommended value of the relative error
163: for the Runge-Kutta method..
164: */
165: static double M_recommended_relerr;
166: /*
1.19 takayama 167: max of beta(i)*x/2
168: */
169: static double M_beta_i_x_o2_max;
170: /*
171: minimum of |beta_i-beta_j|
172: */
173: static double M_beta_i_beta_j_min;
174: /*
175: Value of matrix hg
176: */
177: static double M_mh_t_value;
1.25 takayama 178: /*
179: Show the process of updating degree.
180: */
181: int M_show_autosteps=1;
1.16 takayama 182:
1.2 takayama 183: /* prototypes */
184: static void *mymalloc(int size);
1.10 takayama 185: static int myfree(void *p);
186: static int myerror(char *s);
1.2 takayama 187: static double jack1(int K);
188: static double jack1diff(int k);
189: static double xval(int i,int p); /* x_i^p */
190: static int mysum(int L[]);
191: static int plength(int P[]);
192: static int plength_t(int P[]);
1.3 takayama 193: static void ptrans(int P[M_nmx],int Pt[]);
1.2 takayama 194: static int huk(int K[],int I,int J);
195: static int hdk(int K[],int I,int J);
196: static double jjk(int K[]);
197: static double ppoch(double A,int K[]);
198: static double ppoch2(double A,double B,int K[]);
199: static double mypower(double x,int n);
200: static double qk(int K[],double A[A_LEN],double B[B_LEN]);
201: static int bb(int N[],int K[],int M[],int I,int J);
202: static double beta(int K[],int M[]);
1.10 takayama 203: static int printp(int kappa[]);
1.2 takayama 204: static double q3_10(int K[],int M[],int SK);
205: static int nk(int KK[]);
206: static int myeq(int P1[],int P2[]);
207: static int pListPartition(int M,int N);
208: static int pListPartition2(int Less,int From,int To, int M);
209: static void pExec_0();
210: static int pListHS(int Kap[],int N);
211: static int pListHS2(int From,int To,int Kap[]);
212: static void hsExec_0();
213: static int pmn(int M,int N);
214: static int *cloneP(int a[]);
1.10 takayama 215: static int copyP(int p[],int a[]);
1.2 takayama 216: static void pExec_darray(void);
1.10 takayama 217: static int genDarray2(int M,int N);
218: static int isHStrip(int Kap[],int Nu[]);
1.2 takayama 219: static void hsExec_beta(void);
1.10 takayama 220: static int genBeta(int Kap[]);
1.2 takayama 221: static int psublen(int Kap[],int Mu[]);
1.10 takayama 222: static int genJack(int M,int N);
1.2 takayama 223:
224: static int imypower(int x,int n);
1.10 takayama 225: static int usage();
226: static int setParamDefault();
227: static int next(struct SFILE *fp,char *s,char *msg);
1.2 takayama 228: static int setParam(char *fname);
1.8 takayama 229: static int showParam(struct SFILE *fp,int fd);
1.2 takayama 230: static double iv_factor(void);
231: static double gammam(double a,int n);
232: static double mypower(double a,int n);
233:
1.40 ! takayama 234: static double iv_factor_ef_type1(void);
! 235: static double iv_factor_ef_type2(void);
! 236: static void setM_x(void);
! 237: static void setM_x_ef_type1(void);
! 238: static void setM_x_ef_type2(void);
! 239:
1.32 takayama 240: #ifdef STANDALONE
241: static int test_ptrans();
242: static int printp2(int kappa[]);
243: static int test_beta();
244: static void mtest4();
245: static void mtest4b();
246: static int plength2(int P1[],int P2[]);
247: static int checkBeta1();
248: static int checkJack1(int M,int N);
249: static int checkJack2(int M,int N);
250: static int mtest1b();
251: static double q3_5(double A[],double B[],int K[],int I);
252: #endif
253:
1.2 takayama 254: double mh_t(double A[A_LEN],double B[B_LEN],int N,int M);
255: double mh_t2(int J);
1.3 takayama 256: struct MH_RESULT *jk_main(int argc,char *argv[]);
1.16 takayama 257: struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree);
1.3 takayama 258: int jk_freeWorkArea();
259: int jk_initializeWorkArea();
260:
261: int jk_freeWorkArea() {
262: /* bug, p in the cloneP will not be deallocated.
1.12 takayama 263: Nk in genDarray2 will not be deallocated.
264: */
1.3 takayama 265: int i;
1.6 takayama 266: JK_deallocate=1;
1.3 takayama 267: if (Darray) {myfree(Darray); Darray=NULL;}
268: if (Parray) {myfree(Parray); Parray=NULL;}
269: if (ParraySize) {myfree(ParraySize); ParraySize=NULL;}
270: if (M_beta_0) {myfree(M_beta_0); M_beta_0=NULL;}
271: if (M_beta_1) {myfree(M_beta_1); M_beta_1=NULL;}
272: if (M_jack) {
1.12 takayama 273: for (i=0; M_jack[i] != NULL; i++) {
1.29 takayama 274: if (Debug) oxprintf("Free M_jack[%d]\n",i);
1.12 takayama 275: myfree(M_jack[i]); M_jack[i] = NULL;
276: }
277: myfree(M_jack); M_jack=NULL;
1.3 takayama 278: }
279: if (M_qk) {myfree(M_qk); M_qk=NULL;}
280: if (P_pki) {myfree(P_pki); P_pki=NULL;}
1.6 takayama 281: JK_deallocate=0;
1.13 takayama 282: return(0);
1.3 takayama 283: }
284: int jk_initializeWorkArea() {
285: int i,j;
1.6 takayama 286: JK_deallocate=1;
287: xval(0,0);
288: JK_deallocate=0;
1.3 takayama 289: Darray=NULL;
290: Parray=NULL;
291: ParraySize=NULL;
292: M_beta_0=NULL;
293: M_beta_1=NULL;
294: M_jack=NULL;
295: M_qk=NULL;
296: for (i=0; i<M_nmx; i++) M_kap[i]=HS_mu[i]=0;
297: for (i=0; i<M_nmx; i++) M_x[i]=0;
298: for (i=0; i<M_nmx; i++) for (j=0; j<M_m_MAX; j++) Xarray[i][j]=0;
1.5 takayama 299: for (i=0; i<M_nmx; i++) M_beta_kap[i]=0;
1.3 takayama 300: M_m=M_m_MAX-2;
301: Alpha = 2;
302: HS_n=M_nmx;
303: P_pki=NULL;
304: P_pmn=0;
305: DR_parray=0;
306: M_beta_pt=0;
307: M_df=1;
308: M_2n=0;
309: M_rel_error=0.0;
1.4 takayama 310: Sample = Sample_default;
1.5 takayama 311: Xng=0.0;
312: M_n=0;
1.13 takayama 313: return(0);
1.3 takayama 314: }
1.2 takayama 315:
316: static void *mymalloc(int size) {
1.1 takayama 317: void *p;
1.29 takayama 318: if (Debug) oxprintf("mymalloc(%d)\n",size);
1.8 takayama 319: p = (void *)mh_malloc(size);
1.1 takayama 320: if (p == NULL) {
1.29 takayama 321: oxprintfe("No more memory.\n");
1.3 takayama 322: mh_exit(-1);
1.1 takayama 323: }
324: return(p);
325: }
1.10 takayama 326: static int myfree(void *p) {
1.29 takayama 327: if (Debug) oxprintf("myFree at %p\n",p);
1.13 takayama 328: return(mh_free(p));
1.7 takayama 329: }
1.29 takayama 330: static int myerror(char *s) { oxprintfe("%s: type in control-C\n",s); getchar(); getchar(); return(0);}
1.1 takayama 331:
1.2 takayama 332: static double jack1(int K) {
1.1 takayama 333: double F;
334: extern int Alpha;
1.32 takayama 335: int J,II,JJ,N; /* int I,J,L,II,JJ,N; */
1.1 takayama 336: N = 1;
337: if (K == 0) return((double)1);
338: F = xval(1,K);
339: for (J=0; J<K; J++) {
340: II = 1; JJ = J+1;
341: F *= (N-(II-1)+Alpha*(JJ-1));
342: }
343: return(F);
344: }
1.2 takayama 345: static double jack1diff(int K) {
1.1 takayama 346: double F;
347: extern int Alpha;
1.32 takayama 348: int J,II,JJ,N; /* int I,J,S,L,II,JJ,N; */
1.1 takayama 349: N = 1;
350: if (K == 0) return((double)1);
351: F = K*xval(1,K-1);
352: for (J=0; J<K; J++) {
353: II = 1; JJ = J+1;
354: F *= (N-(II-1)+Alpha*(JJ-1));
355: }
356: return(F);
357: }
358:
1.2 takayama 359: static double xval(int ii,int p) { /* x_i^p */
1.1 takayama 360: extern double M_x[];
361: int i,j;
1.10 takayama 362: static int init=0;
1.6 takayama 363: if (JK_deallocate) { init=0; return(0.0);}
1.1 takayama 364: if (!init) {
365: for (i=1; i<=M_n; i++) {
366: for (j=0; j<M_m_MAX; j++) {
1.12 takayama 367: if (j != 0) {
368: Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1];
369: }else{
370: Xarray[i-1][0] = 1;
371: }
1.1 takayama 372: }
373: }
374: init = 1;
375: }
376: if (ii < 1) myerror("xval, index out of bound.");
377: if (p > M_m_MAX-2) myerror("xval, p is too large.");
378: if (p < 0) {
379: myerror("xval, p is negative.");
1.29 takayama 380: oxprintf("ii=%d, p=%d\n",ii,p);
1.3 takayama 381: mh_exit(-1);
1.1 takayama 382: }
383: return(Xarray[ii-1][p]);
384: }
385:
1.2 takayama 386: static int mysum(int L[]) {
1.1 takayama 387: int S,I,N;
388: N=M_n;
389: S=0;
390: for (I=0; I<N; I++) S += L[I];
391: return(S);
392: }
393:
394: /*
1.12 takayama 395: (3,2,2,0,0) --> 3
1.1 takayama 396: */
1.2 takayama 397: static int plength(int P[]) {
1.1 takayama 398: int I;
399: for (I=0; I<M_n; I++) {
400: if (P[I] == 0) return(I);
401: }
402: return(M_n);
403: }
404: /* plength for transpose */
1.2 takayama 405: static int plength_t(int P[]) {
1.1 takayama 406: int I;
407: for (I=0; I<M_m; I++) {
408: if (P[I] == 0) return(I);
409: }
410: return(M_m);
411: }
412:
413: /*
1.12 takayama 414: ptrans(P) returns Pt
1.1 takayama 415: */
1.3 takayama 416: static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] */
1.1 takayama 417: extern int M_m;
418: int i,j,len;
1.3 takayama 419: int p[M_nmx];
1.1 takayama 420: for (i=0; i<M_n; i++) p[i] = P[i];
421: for (i=0; i<M_m+1; i++) Pt[i] = 0;
422: for (i=0; i<M_m; i++) {
423: len=plength(p); Pt[i] = len;
424: if (len == 0) return;
425: for (j=0; j<len; j++) p[j] -= 1;
426: }
427: }
428:
1.32 takayama 429: #ifdef STANDALONE
1.10 takayama 430: static int test_ptrans() {
1.1 takayama 431: extern int M_m;
432: int p[M_n0]={5,3,2};
433: int pt[10];
434: int i;
435: M_m = 10;
436: ptrans(p,pt);
1.29 takayama 437: if (Debug) {for (i=0; i<10; i++) oxprintf("%d,",pt[i]); oxprintf("\n");}
1.13 takayama 438: return(0);
1.1 takayama 439: }
1.32 takayama 440: #endif
1.1 takayama 441:
442: /*
443: upper hook length
444: h_kappa^*(K)
445: */
1.2 takayama 446: static int huk(int K[],int I,int J) {
1.1 takayama 447: extern int Alpha;
448: int Kp[M_m_MAX];
449: int A,H;
450: A=Alpha;
1.29 takayama 451: /*oxprintf("h^k(%a,%a,%a,%a)\n",K,I,J,A);*/
1.1 takayama 452: ptrans(K,Kp);
453: H=Kp[J-1]-I+A*(K[I-1]-J+1);
454: return(H);
455: }
456:
457: /*
458: lower hook length
459: h^kappa_*(K)
460: */
1.2 takayama 461: static int hdk(int K[],int I,int J) {
1.1 takayama 462: extern int Alpha;
463: int Kp[M_m_MAX];
464: int A,H;
465: A = Alpha;
1.29 takayama 466: /*oxprintf("h_k(%a,%a,%a,%a)\n",K,I,J,A);*/
1.1 takayama 467: ptrans(K,Kp);
468: H=Kp[J-1]-I+1+A*(K[I-1]-J);
469: return(H);
470: }
471: /*
472: j_kappa. cf. Stanley.
473: */
1.2 takayama 474: static double jjk(int K[]) {
1.1 takayama 475: extern int Alpha;
1.33 takayama 476: int L,I,J;
1.1 takayama 477: double V;
1.33 takayama 478:
1.1 takayama 479: V=1;
480: L=plength(K);
481: for (I=0; I<L; I++) {
482: for (J=0; J<K[I]; J++) {
483: V *= huk(K,I+1,J+1)*hdk(K,I+1,J+1);
484: }
485: }
1.29 takayama 486: if (Debug) {printp(K); oxprintf("<--K, jjk=%lg\n",V);}
1.1 takayama 487: return(V);
488: }
489: /*
490: (a)_kappa^\alpha, Pochhammer symbol
491: Note that ppoch(a,[n]) = (a)_n, Alpha=2
492: */
1.2 takayama 493: static double ppoch(double A,int K[]) {
1.1 takayama 494: extern int Alpha;
495: double V;
496: int L,I,J,II,JJ;
497: V = 1;
498: L=plength(K);
499: for (I=0; I<L; I++) {
500: for (J=0; J<K[I]; J++) {
501: II = I+1; JJ = J+1;
502: V *= (A-((double)(II-1))/((double)Alpha)+JJ-1);
503: }
504: }
505: return(V);
506: }
1.2 takayama 507: static double ppoch2(double A,double B,int K[]) {
1.1 takayama 508: extern int Alpha;
509: double V;
510: int L,I,J,II,JJ;
511: V = 1;
512: L=plength(K);
513: for (I=0; I<L; I++) {
514: for (J=0; J<K[I]; J++) {
515: II = I+1; JJ = J+1;
516: V *= (A-((double)(II-1))/((double)Alpha)+JJ-1);
517: V /= (B-((double)(II-1))/((double)Alpha)+JJ-1);
518: }
519: }
520: return(V);
521: }
1.2 takayama 522: static double mypower(double x,int n) {
1.1 takayama 523: int i;
524: double v;
525: if (n < 0) return(1/mypower(x,-n));
526: v = 1;
527: for (i=0; i<n; i++) v *= x;
528: return(v);
529: }
530: /* Q_kappa
1.12 takayama 531: */
1.2 takayama 532: static double qk(int K[],double A[A_LEN],double B[B_LEN]) {
1.1 takayama 533: extern int Alpha;
534: int P,Q,I;
535: double V;
536: P = A_LEN;
537: Q = B_LEN;
538: V = mypower((double) Alpha,mysum(K))/jjk(K);
1.2 takayama 539: /* to reduce numerical errors, temporary. */
1.1 takayama 540: if (P == Q) {
541: for (I=0; I<P; I++) V = V*ppoch2(A[I],B[I],K);
1.38 takayama 542: return(V);
1.1 takayama 543: }
544:
1.38 takayama 545: if (P > Q) {
546: for (I=0; I<Q; I++) V = V*ppoch2(A[I],B[I],K);
547: for (I=Q; I<P; I++) V = V*ppoch(A[I],K);
548: }else {
549: for (I=0; I<P; I++) V = V*ppoch2(A[I],B[I],K);
550: for (I=P; I<Q; I++) V = V/ppoch(B[I],K);
551: }
552: /* for debug
553: printf("K=");
554: for (I=0; I<3; I++) printf("%d, ",K[I]);
555: printf("qk=%lg\n",V);
556: */
1.1 takayama 557: return(V);
558: }
559:
560: /*
1.12 takayama 561: B^nu_{kappa,mu}(i,j)
562: bb(N,K,M,I,J)
1.1 takayama 563: */
1.2 takayama 564: static int bb(int N[],int K[],int M[],int I,int J) {
1.1 takayama 565: int Kp[M_m_MAX]; int Mp[M_m_MAX];
566: ptrans(K,Kp);
567: ptrans(M,Mp);
568:
569: /*
1.29 takayama 570: printp(K); oxprintf("K<--, "); printp2(Kp); oxprintf("<--Kp\n");
571: printp(M); oxprintf("M<--, "); printp2(Mp); oxprintf("<--Mp\n");
1.1 takayama 572: */
573:
574: if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J));
575: if (Kp[J-1] == Mp[J-1]) return(huk(N,I,J));
576: else return(hdk(N,I,J));
577: }
578: /*
579: beta_{kappa,mu}
580: beta(K,M)
581: */
1.2 takayama 582: static double beta(int K[],int M[]) {
1.1 takayama 583: double V;
584: int L,I,J,II,JJ;
585: V = 1;
586:
587: L=plength(K);
588: for (I=0; I<L; I++) {
589: for (J=0; J<K[I]; J++) {
590: II = I+1; JJ = J+1;
591: V *= (double)bb(K,K,M,II,JJ);
1.29 takayama 592: /* oxprintf("[%d,%d,%lf]\n",I,J,V); */
1.1 takayama 593: }
594: }
595:
596: L=plength(M);
597: for (I=0; I<L; I++) {
598: for (J=0; J<M[I]; J++) {
599: II = I+1; JJ = J+1;
600: V /= (double)bb(M,K,M,II,JJ);
1.29 takayama 601: /* oxprintf("[%d,%d,%lf]\n",I,J,V);*/
1.1 takayama 602: }
603: }
604:
605: return(V);
606: }
1.10 takayama 607: static int printp(int kappa[]) {
1.1 takayama 608: int i;
1.29 takayama 609: oxprintf("(");
1.1 takayama 610: for (i=0; i<M_n; i++) {
1.29 takayama 611: if (i <M_n-1) oxprintf("%d,",kappa[i]);
612: else oxprintf("%d)",kappa[i]);
1.1 takayama 613: }
1.13 takayama 614: return(0);
1.1 takayama 615: }
1.32 takayama 616: #ifdef STANDALONE
1.10 takayama 617: static int printp2(int kappa[]) {
1.1 takayama 618: int i,ell;
1.29 takayama 619: oxprintf("(");
1.1 takayama 620: ell = plength_t(kappa);
621: for (i=0; i<ell+1; i++) {
1.29 takayama 622: if (i <ell+1-1) oxprintf("%d,",kappa[i]);
623: else oxprintf("%d)",kappa[i]);
1.1 takayama 624: }
1.13 takayama 625: return(0);
1.1 takayama 626: }
627:
1.10 takayama 628: static int test_beta() {
1.1 takayama 629: int kappa[M_n0]={2,1,0};
630: int mu1[M_n0]={1,0,0};
631: int mu2[M_n0]={1,1,0};
632: int mu3[M_n0]={2,0,0};
1.29 takayama 633: printp(kappa); oxprintf(","); printp(mu3); oxprintf(": beta = %lf\n",beta(kappa,mu3));
634: printp(kappa); oxprintf(","); printp(mu1); oxprintf(": beta = %lf\n",beta(kappa,mu1));
1.30 takayama 635: printp(kappa); oxprintf(","); printp(mu2); oxprintf(": beta = %lf\n",beta(kappa,mu2)); return(0);
1.1 takayama 636: }
1.32 takayama 637: #endif
1.1 takayama 638: /* main() { test_beta(); } */
639:
640:
641: /*
1.12 takayama 642: cf. w1m.rr
1.1 takayama 643: matrix hypergeometric by jack
644: N variables, up to degree M.
645: */
646: /* todo
1.12 takayama 647: def mhgj(A,B,N,M) {
648: F = 0;
649: P = partition_a(N,M);
650: for (I=0; I<length(P); I++) {
651: K = P[I];
652: F += qk(K,A,B)*jack(K,N);
653: }
654: return(F);
655: }
1.1 takayama 656: */
657:
658:
659: /* The quotient of (3.10) of Koev-Edelman K=kappa, M=mu, SK=k */
1.2 takayama 660: static double q3_10(int K[],int M[],int SK) {
1.1 takayama 661: extern int Alpha;
662: int Mp[M_m_MAX];
1.33 takayama 663: // int ML[M_nmx];
1.3 takayama 664: int N[M_nmx];
1.1 takayama 665: int i,R;
666: double T,Q,V,Ur,Vr,Wr;
667: ptrans(M,Mp);
1.33 takayama 668: for (i=0; i<M_n; i++) {N[i] = M[i];}
1.1 takayama 669: N[SK-1] = N[SK-1]-1;
670:
671: T = SK-Alpha*M[SK-1];
672: Q = T+1;
673: V = Alpha;
674: for (R=1; R<=SK; R++) {
675: Ur = Q-R+Alpha*K[R-1];
676: V *= Ur/(Ur+Alpha-1);
677: }
678: for (R=1; R<=SK-1; R++) {
679: Vr = T-R+Alpha*M[R-1];
680: V *= (Vr+Alpha)/Vr;
681: }
682: for (R=1; R<=M[SK-1]-1; R++) {
683: Wr = Mp[R-1]-T-Alpha*R;
684: V *= (Wr+Alpha)/Wr;
685: }
686: return(V);
687: }
688:
1.32 takayama 689: #ifdef STANDALONE
1.2 takayama 690: static double q3_5(double A[],double B[],int K[],int I) {
1.1 takayama 691: extern int Alpha;
692: int Kp[M_m_MAX];
693: double C,D,V,Ej,Fj,Gj,Hj,Lj;
694: int J,P,Q;
695: ptrans(K,Kp);
696: P=A_LEN;; Q = B_LEN;
697: C = -((double)(I-1))/Alpha+K[I-1]-1;
698: D = K[I-1]*Alpha-I;
699:
700: V=1;
701:
702: for (J=1; J<=P; J++) {
1.12 takayama 703: V *= (A[J-1]+C);
1.1 takayama 704: }
705: for (J=1; J<=Q; J++) {
1.12 takayama 706: V /= (B[J-1]+C);
1.1 takayama 707: }
708:
709: for (J=1; J<=K[I-1]-1; J++) {
1.12 takayama 710: Ej = D-J*Alpha+Kp[J-1];
711: Gj = Ej+1;
712: V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha));
1.1 takayama 713: }
714: for (J=1; J<=I-1; J++) {
1.12 takayama 715: Fj=K[J-1]*Alpha-J-D;
716: Hj=Fj+Alpha;
717: Lj=Hj*Fj;
718: V *= (Lj-Fj)/(Lj+Hj);
1.1 takayama 719: }
720: return(V);
721: }
1.32 takayama 722: #endif
723: #ifdef STANDALONE
1.2 takayama 724: static void mtest4() {
1.1 takayama 725: double A[A_LEN] = {1.5};
726: double B[B_LEN]={6.5};
727: int K[M_n0] = {3,2,0};
728: int I=2;
729: int Ki[M_n0]={3,1,0};
730: double V1,V2;
731: V1=q3_5(A,B,K,I);
732: V2=qk(K,A,B)/qk(Ki,A,B);
1.29 takayama 733: oxprintf("%lf== %lf?\n",V1,V2);
1.1 takayama 734: }
1.2 takayama 735: static void mtest4b() {
1.1 takayama 736: int K[M_n0]={3,2,0};
737: int M[M_n0]={2,1,0};
738: int N[M_n0]={2,0};
739: int SK=2;
740: double V1,V2;
741: V1=q3_10(K,M,SK);
742: V2=beta(K,N)/beta(K,M);
1.29 takayama 743: oxprintf("%lf== %lf?\n",V1,V2);
1.1 takayama 744: }
1.32 takayama 745: #endif
1.1 takayama 746:
747: /* main() { mtest4(); mtest4b(); } */
748:
749: /* nk in (4.1),
1.12 takayama 750: */
1.2 takayama 751: static int nk(int KK[]) {
1.1 takayama 752: extern int *Darray;
753: int N,I,Ki;
1.3 takayama 754: int Kpp[M_nmx];
1.1 takayama 755: int i;
756: N = plength(KK);
757: if (N == 0) return(0);
758: if (N == 1) return(KK[0]);
759: for (i=0; i<M_n; i++) Kpp[i] = 0;
760: for (I=0; I<N-1; I++) Kpp[I] = KK[I];
761: Ki = KK[N-1];
762: /* K = (Kpp,Ki) */
763: return(Darray[nk(Kpp)]+Ki-1);
764: }
1.32 takayama 765: #ifdef STANDALONE
1.2 takayama 766: static int plength2(int P1[],int P2[]) {
1.1 takayama 767: int S1,S2;
768: S1 = plength(P1); S2 = plength(P2);
769: if (S1 > S2) return(1);
770: else if (S1 == S2) {
771: S1=mysum(P1); S2=mysum(P2);
772: if(S1 > S2) return(1);
773: else if (S1 == S2) return(0);
774: else return(-1);
775: }
776: else return(-1);
777: }
1.32 takayama 778: #endif
1.2 takayama 779: static int myeq(int P1[],int P2[]) {
1.1 takayama 780: int I,L1;
781: if ((L1=plength(P1)) != plength(P2)) return(0);
782: for (I=0; I<L1; I++) {
783: if (P1[I] != P2[I]) return(0);
784: }
785: return(1);
786: }
787: /*
788: M is a degree, N is a number of variables
789: genDarray(3,3);
790: N(0)=0;
791: N(1)=1;
792: N(2)=2;
793: N(3)=3;
794: N(1,1)=4; D[1] = 4
795: N(2,1)=5; D[2] = 5;
796: N(1,1,1)=6; D[4] = 6;
797: still buggy.
798: */
799:
1.2 takayama 800: static int pListPartition(int M,int N) {
1.1 takayama 801: extern int M_m;
802: extern int M_kap[];
803: int I;
804: /* initialize */
805: if (M_n != N) {
1.29 takayama 806: oxprintfe("M_n != N\n"); mh_exit(-1);
1.1 takayama 807: }
808: M_m = M;
809: /* M_plist = []; */
810: /* end of initialize */
811: (*M_pExec)(); /* exec for 0 */
812: for (I=1; I<=M_n; I++) {
813: pListPartition2(M_m,1,I,M_m);
814: }
815: /* M_plist = reverse(M_plist); */
816: return(1);
817: }
818:
819: /*
820: Enumerate all such that
821: Less >= M_kap[From], ..., M_kap[To], |(M_kap[From],...,M_kap[To])|<=M,
822: */
1.2 takayama 823: static int pListPartition2(int Less,int From,int To, int M) {
1.33 takayama 824: int I;
1.11 takayama 825: mh_check_intr(100);
1.1 takayama 826: if (To < From) {
1.12 takayama 827: (*M_pExec)(); return(0);
1.1 takayama 828: }
829: for (I=1; (I<=Less) && (I<=M) ; I++) {
830: M_kap[From-1] = I;
1.33 takayama 831: pListPartition2(I,From+1,To,M-I);
1.1 takayama 832: }
833: return(1);
834: }
835:
836: /*
1.33 takayama 837: Commands to do for each partition are given here.
1.1 takayama 838: */
1.2 takayama 839: static void pExec_0() {
1.1 takayama 840: if (Debug) {
1.29 takayama 841: oxprintf("M_kap=");
1.12 takayama 842: printp(M_kap);
1.29 takayama 843: oxprintf("\n");
1.1 takayama 844: }
845: }
846:
847: /* Test.
1.12 takayama 848: Compare pListPartition(4,3); genDarray(4,3);
849: Compare pListPartition(5,3); genDarray(5,3);
1.1 takayama 850:
851: */
852:
853: /*
1.12 takayama 854: main() {
1.1 takayama 855: M_pExec = pExec_0;
856: pListPartition(5,3);
1.12 takayama 857: }
1.1 takayama 858: */
859:
860:
861: /*
862: List all horizontal strips.
863: Kap[0] is not a dummy in C code. !(Start from Kap[1].)
864: */
1.2 takayama 865: static int pListHS(int Kap[],int N) {
1.1 takayama 866: extern int HS_n;
867: extern int HS_mu[];
868: int i;
869: HS_n = N;
870: /* Clear HS_mu. Do not forget when N < M_n */
871: for (i=0; i<M_n; i++) HS_mu[i] = 0;
1.13 takayama 872: return(pListHS2(1,N,Kap));
1.1 takayama 873: }
874:
1.2 takayama 875: static int pListHS2(int From,int To,int Kap[]) {
1.1 takayama 876: int More,I;
877: if (To <From) {(*HS_hsExec)(); return(0);}
878: if (From == HS_n) More=0; else More=Kap[From];
879: for (I=Kap[From-1]; I>= More; I--) {
880: HS_mu[From-1] = I;
881: pListHS2(From+1,To,Kap);
882: }
883: return(1);
884: }
885:
1.2 takayama 886: static void hsExec_0() {
1.32 takayama 887: /* int i; */
1.29 takayama 888: if(Debug) {oxprintf("hsExec: "); printp(HS_mu); oxprintf("\n");}
1.1 takayama 889: }
890:
891: /*
892: pListHS([0,4,2,1],3);
893: */
894: /*
1.12 takayama 895: main() {
1.1 takayama 896: int Kap[3]={4,2,1};
897: HS_hsExec = hsExec_0;
898: pListHS(Kap,3);
1.12 takayama 899: }
1.1 takayama 900: */
901:
902: /* The number of partitions <= M, with N parts.
1.12 takayama 903: (0,0,...,0) is excluded.
1.1 takayama 904: */
905: #define aP_pki(i,j) P_pki[(i)*(M+1)+(j)]
1.2 takayama 906: static int pmn(int M,int N) {
1.1 takayama 907: int Min_m_n,I,K,S,T,i,j;
908: extern int P_pmn;
909: extern int *P_pki;
910: Min_m_n = (M>N?N:M);
911: /* P_pki=newmat(Min_m_n+1,M+1); */
912: P_pki = (int *) mymalloc(sizeof(int)*(Min_m_n+1)*(M+1));
913: for (i=0; i<Min_m_n+1; i++) for (j=0; j<M+1; j++) aP_pki(i,j) = 0;
914: for (I=1; I<=M; I++) aP_pki(1,I) = 1;
915: for (K=1; K<=Min_m_n; K++) aP_pki(K,0) = 0;
916: S = M;
917: for (K=2; K<=Min_m_n; K++) {
918: for (I=1; I<=M; I++) {
919: if (I-K < 0) T=0; else T=aP_pki(K,I-K);
920: aP_pki(K,I) = aP_pki(K-1,I-1)+T;
921: S += aP_pki(K,I);
922: }
923: }
924: P_pmn=S;
925: if (Debug) {
1.29 takayama 926: oxprintf("P_pmn=%d\n",P_pmn);
1.1 takayama 927: for (i=0; i<=Min_m_n; i++) {
1.29 takayama 928: for (j=0; j<=M; j++) oxprintf("%d,",aP_pki(i,j));
929: oxprintf("\n");
1.1 takayama 930: }
931: }
932: myfree(P_pki); P_pki=NULL;
933: return(S);
934: }
935:
936: /*
1.29 takayama 937: main() {pmn(4,3); oxprintf("P_pmn=%d\n",P_pmn);}
1.1 takayama 938: */
939:
1.2 takayama 940: static int *cloneP(int a[]) {
1.1 takayama 941: int *p;
942: int i;
943: p = (int *) mymalloc(sizeof(int)*M_n);
944: for (i=0; i<M_n; i++) p[i] = a[i];
945: return(p);
946: }
1.10 takayama 947: static int copyP(int p[],int a[]) {
1.1 takayama 948: int i;
949: for (i=0; i<M_n; i++) p[i] = a[i];
1.13 takayama 950: return(0);
1.1 takayama 951: }
952:
1.2 takayama 953: static void pExec_darray(void) {
1.1 takayama 954: extern int DR_parray;
955: extern int M_kap[];
956: extern int **Parray;
957: extern int *ParraySize;
958: int *K;
959: pExec_0();
960: K = cloneP(M_kap);
961: Parray[DR_parray] = K;
962: ParraySize[DR_parray] = mysum(K);
963: DR_parray++;
964: }
1.10 takayama 965: static int genDarray2(int M,int N) {
1.1 takayama 966: extern int *Darray;
967: extern int **Parray;
968: extern int DR_parray;
969: extern int M_m;
970: int Pmn,I,J,Ksize,i;
971: int **L;
972: int *Nk;
973: int *K;
1.3 takayama 974: int Kone[M_nmx];
1.1 takayama 975:
976: M_m = M;
977: Pmn = pmn(M,N)+1;
1.29 takayama 978: if (Debug) oxprintf("Degree M = %d, N of vars N = %d, Pmn+1=%d\n",M,N,Pmn);
1.1 takayama 979: Darray=(int *) mymalloc(sizeof(int)*Pmn);
980: for (i=0; i<Pmn; i++) Darray[i] = 0;
981: Parray=(int **) mymalloc(sizeof(int *)*Pmn);
982: for (i=0; i<Pmn; i++) Parray[i] = NULL;
983: ParraySize=(int *) mymalloc(sizeof(int *)*Pmn);
984: for (i=0; i<Pmn; i++) ParraySize[i] = 0;
985: DR_parray=0;
986: M_pExec = pExec_darray;
987: pListPartition(M,N); /* pExec_darray() is executed for all partitions */
988: L = Parray;
989:
990: Nk = (int *) mymalloc(sizeof(int)*(Pmn+1));
991: for (I=0; I<Pmn; I++) Nk[I] = I;
992: for (I=0; I<Pmn; I++) {
1.11 takayama 993: mh_check_intr(100);
1.12 takayama 994: K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */
995: Ksize = plength(K);
996: if (Ksize >= M_n) {
1.29 takayama 997: if (Debug) {oxprintfe("Ksize >= M_n\n");}
1.12 takayama 998: continue;
999: }
1000: for (i=0; i<M_n; i++) Kone[i] = 0;
1001: for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1;
1002: for (J=0; J<Pmn; J++) {
1003: if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */
1004: }
1.1 takayama 1005: }
1006: if (Debug) {
1.29 takayama 1007: oxprintf("Darray=\n");
1008: for (i=0; i<Pmn; i++) oxprintf("%d\n",Darray[i]);
1009: oxprintf("-----------\n");
1.1 takayama 1010: }
1.13 takayama 1011: return(0);
1.1 takayama 1012: }
1013:
1014:
1015: /* main() { genDarray2(4,3);} */
1016:
1017: /* M_beta_0[*] value of beta_{kappa,mu}, M_beta_1[*] N_mu */
1.10 takayama 1018: static int isHStrip(int Kap[],int Nu[]) {
1.1 takayama 1019: int N1,N2,I,P;
1020: N1 = plength(Kap); N2 = plength(Nu);
1021: if (N2 > N1) return(0);
1022: for (I=0; I<N2; I++) {
1023: if (I >= N1-1) P = 0; else P=Kap[I+1];
1024: if (Kap[I] < Nu[I]) return(0);
1025: if (Nu[I] < P) return(0);
1026: }
1027: return(1);
1028: }
1029:
1.2 takayama 1030: static void hsExec_beta(void) {
1.1 takayama 1031: int *Mu;
1032: int N,Nmu,Nnu,Done,J,K,OK,I,RR;
1033: int Kapt[M_m_MAX];
1034: int Nut[M_m_MAX];
1.3 takayama 1035: int Nu[M_nmx];
1.1 takayama 1036: int rrMax;
1037: hsExec_0();
1.29 takayama 1038: /* oxprintf("M_beta_pt=%a\n",M_beta_pt); */
1.1 takayama 1039: /* Mu = cdr(vtol(HS_mu)); */
1040: Mu = HS_mu; /* buggy? need cloneP */
1041: if (M_beta_pt == 0) {
1042: M_beta_0[0] = 1; M_beta_1[0] = nk(Mu);
1043: M_beta_pt++; return;
1044: }
1045:
1046: N = HS_n;
1047: Nmu = nk(Mu);
1048: M_beta_1[M_beta_pt] = Nmu;
1049: ptrans(M_beta_kap,Kapt);
1050: /* Mu, Nu is exchanged in this code. cf. the K-E paper */
1051: copyP(Nu,Mu); /* buggy need clone? */
1052: for (I=0; I<N; I++) {
1053: Nu[I]++;
1054: if (!isHStrip(M_beta_kap,Nu)) {Nu[I]--; continue;}
1055: Nnu = nk(Nu);
1056: ptrans(Nu,Nut);
1057: Done=0;
1058: for (J=M_beta_pt-1; J>=0; J--) {
1059: if (M_beta_1[J] == Nnu) {
1.12 takayama 1060: K=I+1;
1061: if (Debug) {
1.29 takayama 1062: oxprintf("Found at J=%d, K=%d, q3_10(Kap,Nu,K)=%lf,Nu,Mu= \n",
1.12 takayama 1063: J,K,q3_10(M_beta_kap,Nu,K));
1.29 takayama 1064: printp(Nu); oxprintf("\n");
1065: printp(Mu); oxprintf("\n");
1.12 takayama 1066: }
1067: /* Check other conditions. See Numata's mail on Dec 24, 2011. */
1068: rrMax = Nu[I]-1;
1069: if ((plength_t(Kapt) < rrMax) || (plength_t(Nut) < rrMax)) {
1.29 takayama 1070: if (Debug) oxprintf(" is not taken (length). \n");
1.12 takayama 1071: break;
1072: }
1073: OK=1;
1074: for (RR=0; RR<rrMax; RR++) {
1075: if (Kapt[RR] != Nut[RR]) { OK=0; break;}
1076: }
1.29 takayama 1077: if (!OK) { if (Debug) oxprintf(" is not taken.\n"); break; }
1.12 takayama 1078: /* check done. */
1079: M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K);
1080: Done = 1; break;
1.1 takayama 1081: }
1082: }
1083: if (Done) break; else Nu[I]--;
1084: }
1085: if (!Done) {
1.29 takayama 1086: if (Debug) oxprintf("BUG: not found M_beta_pt=%d.\n",M_beta_pt);
1.13 takayama 1087: /* M_beta_0[M_beta_pt] = NAN; error("Not found."); */
1.1 takayama 1088: M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu);
1089: }
1090: /* Fix the bug of mh.rr */
1091: M_beta_pt++;
1092: }
1.10 takayama 1093: static int genBeta(int Kap[]) {
1.1 takayama 1094: extern double *M_beta_0;
1095: extern int *M_beta_1;
1096: extern int M_beta_pt;
1097: extern int M_beta_kap[];
1098: extern int P_pmn;
1.32 takayama 1099: int I,N;
1.29 takayama 1100: if (Debug) {printp(Kap); oxprintf("<-Kappa, P_pmn=%d\n",P_pmn);}
1.1 takayama 1101: /* M_beta = newmat(2,P_pmn+1); */
1102: M_beta_0 = (double *)mymalloc(sizeof(double)*(P_pmn+1));
1.8 takayama 1103: M_beta_1 = (int *)mymalloc(sizeof(int)*(P_pmn+1));
1.1 takayama 1104: M_beta_pt = 0;
1105: for (I=0; I<=P_pmn; I++) {M_beta_0[I] = NAN; M_beta_1[I] = -1;}
1106: N = plength(Kap);
1107: HS_hsExec = hsExec_beta;
1108: copyP(M_beta_kap,Kap);
1.30 takayama 1109: pListHS(Kap,N); return(0);
1.1 takayama 1110: }
1111: /*
1112: genDarray2(4,3);
1113: genBeta([2,2,0]);
1114: genBeta([2,1,1]);
1115: */
1.32 takayama 1116: #ifdef STANDALONE
1.10 takayama 1117: static int checkBeta1() {
1.1 takayama 1118: int Kap[3] = {2,2,0};
1119: int Kap2[3] = {2,1,0};
1120: int I;
1121: int *Mu;
1122: double Beta_km;
1123: genDarray2(4,3);
1124: genBeta(Kap);
1125: for (I=0; I<M_beta_pt; I++) {
1126: Mu = Parray[M_beta_1[I]];
1127: Beta_km = M_beta_0[I];
1.12 takayama 1128: if (Debug) {
1.29 takayama 1129: printp(Kap); oxprintf("<--Kap, ");
1130: printp(Mu); oxprintf("<--Mu,");
1131: oxprintf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap,Mu));
1.12 takayama 1132: }
1.1 takayama 1133: }
1.29 takayama 1134: if (Debug) oxprintf("-------------------------------------\n");
1.1 takayama 1135: genBeta(Kap2);
1136: for (I=0; I<M_beta_pt; I++) {
1137: Mu = Parray[M_beta_1[I]];
1138: Beta_km = M_beta_0[I];
1.12 takayama 1139: if (Debug) {
1.29 takayama 1140: printp(Kap2); oxprintf("<--Kap, ");
1141: printp(Mu); oxprintf("<--Mu,");
1142: oxprintf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap2,Mu));
1.12 takayama 1143: }
1.1 takayama 1144: }
1.13 takayama 1145: return(0);
1.1 takayama 1146: }
1.32 takayama 1147: #endif
1.1 takayama 1148: /*
1.12 takayama 1149: def checkBeta2() {
1.1 takayama 1150: genDarray2(3,3);
1151: Kap = [2,1,0];
1.29 takayama 1152: oxprintf("Kap=%a\n",Kap);
1.1 takayama 1153: genBeta(Kap);
1154: for (I=0; I<M_beta_pt; I++) {
1.12 takayama 1155: Mu = Parray[M_beta[1][I]];
1156: Beta_km = M_beta[0][I];
1.29 takayama 1157: oxprintf("Mu=%a,",Mu);
1158: oxprintf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));
1.12 takayama 1159: }
1.1 takayama 1160: }
1161: */
1162:
1163: /* main() { checkBeta1(); } */
1164:
1.2 takayama 1165: static int psublen(int Kap[],int Mu[]) {
1.1 takayama 1166: int L1,L2,A,I;
1167: L1 = plength(Kap);
1168: L2 = plength(Mu);
1169: if (L2 > L1) myerror("psub, length mismatches.");
1170: A = 0;
1171: for (I=0; I<L2; I++) {
1172: if (Kap[I] < Mu[I]) myerror("psub, not Kap >= Mu");
1173: A += Kap[I]-Mu[I];
1174: }
1175: for (I=L2; I<L1; I++) A += Kap[I];
1176: return(A);
1177: }
1178:
1179:
1180: /* Table of Jack polynomials
1.12 takayama 1181: Jack[1]* one variable.
1182: Jack[2]* two variables.
1.1 takayama 1183: ...
1.12 takayama 1184: Jack[M_n]* n variables.
1185: Jack[P][J]*
1186: D^J(P variables jack of p variables). Example. J=001 d_1, 010 d_2, 100 d_3
1187: 0<=J<=2^{M_n}-1
1188: Jack[P][J][nk(Kappa)] Jack_Kappa, Kappa is a partition.
1189: 0<=nk(Kappa)<=pmn(M_m,M_n)
1.1 takayama 1190: */
1191:
1192: #define aM_jack(i,j,k) ((M_jack[i])[(j)*(Pmn+1)+(k)])
1.10 takayama 1193: static int genJack(int M,int N) {
1.1 takayama 1194: extern double **M_jack;
1195: extern int M_2n;
1196: extern int P_pmn;
1197: extern int *M_beta_1;
1198: int Pmn,I,J,K,L,Nv,H,P;
1199: int *Kap,*Mu;
1200: double Jack,Beta_km;
1.32 takayama 1201: int Nk,JJ, two_to_I;
1.29 takayama 1202: if (Debug) oxprintf("genJack(%d,%d)\n",M,N);
1.3 takayama 1203: M_jack = (double **) mymalloc(sizeof(double *)*(N+2));
1.1 takayama 1204: M_2n = imypower(2,N);
1205: Pmn = pmn(M,N); /*P_pmn is initializeded.
1206: Warning. It is reset when pmn is called.*/
1207: for (I=0; I<=N; I++) M_jack[I] = (double *)mymalloc(sizeof(double)*(M_2n*(Pmn+1))); /* newmat(M_2n,Pmn+1); */
1.3 takayama 1208: M_jack[N+1] = NULL;
1.1 takayama 1209: genDarray2(M,N); /* Darray, Parray is initialized */
1210: for (I=1; I<=N; I++) aM_jack(I,0,0) = 1;
1211: if (M_df) {
1212: for (I=1; I<=N; I++) {
1213: for (J=1; J<M_2n; J++) aM_jack(I,J,0) = 0;
1214: }
1215: }
1216:
1.3 takayama 1217: /* N must satisfies N > 0 */
1.1 takayama 1218: for (K=1; K<=M; K++) {
1219: aM_jack(1,0,K) = jack1(K);
1220: if (M_df) {
1221: aM_jack(1,1,K) = jack1diff(K); /* diff(jack([K],1),x_1); */
1222: for (J=2; J<M_2n; J++) aM_jack(1,J,K) = 0;
1223: }
1224: }
1.32 takayama 1225: for (I=1; I<=N; I++) { two_to_I = imypower(2,I);
1.1 takayama 1226: for (K=M+1; K<Pmn+1; K++) {
1227: aM_jack(I,0,K) = NAN;
1228: if (M_df) {
1229: for (J=1; J<M_2n; J++) {
1.32 takayama 1230: if (J >= two_to_I) aM_jack(I,J,K) = 0; /* J >= 2^I */
1.12 takayama 1231: else aM_jack(I,J,K) = NAN;
1.1 takayama 1232: }
1233: }
1234: }
1235: }
1236:
1237: /* Start to evaluate the entries of the table */
1238: for (K=1; K<=Pmn; K++) {
1239: Kap = Parray[K]; /* bug. need copy? */
1240: L = plength(Kap);
1241: for (I=1; I<=L-1; I++) {
1242: aM_jack(I,0,K) = 0;
1243: if (M_df) {
1244: for (J=1; J<M_2n; J++) aM_jack(I,J,K) = 0;
1245: }
1246: }
1.29 takayama 1247: if (Debug) {oxprintf("Kappa="); printp(Kap);}
1.1 takayama 1248: /* Enumerate horizontal strip of Kappa */
1249: genBeta(Kap); /* M_beta_pt stores the number of hs */
1250: /* Nv is the number of variables */
1251: for (Nv = (L==1?2:L); Nv <= N; Nv++) {
1252: Jack = 0;
1253: for (H=0; H<M_beta_pt; H++) {
1254: Nk = M_beta_1[H];
1255: Mu = Parray[Nk];
1.12 takayama 1256: if (UseTable) {
1257: Beta_km = M_beta_0[H];
1258: }else{
1259: Beta_km = beta(Kap,Mu);
1260: /* do not use the M_beta table. It's buggy. UseTable is experimental.*/
1261: }
1.29 takayama 1262: if (Debug) {oxprintf("Nv(number of variables)=%d, Beta_km=%lf, Mu=",Nv,Beta_km);
1263: printp(Mu); oxprintf("\n");}
1.1 takayama 1264: P = psublen(Kap,Mu);
1265: Jack += aM_jack(Nv-1,0,Nk)*Beta_km*xval(Nv,P); /* util_v(x,[Nv])^P;*/
1.29 takayama 1266: if (Debug) oxprintf("xval(%d,%d)=%lf\n",Nv,P,xval(Nv,P));
1.1 takayama 1267: }
1268: aM_jack(Nv,0,K) = Jack;
1269: if (M_df) {
1270: /* The case of M_df > 0. */
1271: for (J=1; J<M_2n; J++) {
1.12 takayama 1272: mh_check_intr(100);
1273: Jack = 0;
1274: for (H=0; H<M_beta_pt; H++) {
1275: Nk = M_beta_1[H];
1276: Mu = Parray[Nk];
1277: if (UseTable) {
1278: Beta_km = M_beta_0[H];
1279: }else{
1280: Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */
1281: }
1.29 takayama 1282: if (Debug) {oxprintf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);
1283: printp(Mu); oxprintf("\n"); }
1.12 takayama 1284: P = psublen(Kap,Mu);
1285: if (J & (1 << (Nv-1))) {
1286: JJ = J & ((1 << (Nv-1)) ^ 0xffff); /* NOTE!! Up to 16 bits. mh-15 */
1287: if (P != 0) {
1288: Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1);
1289: }
1290: }else{
1291: Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P);
1292: }
1293: }
1294: aM_jack(Nv,J,K) = Jack;
1.29 takayama 1295: if (Debug) oxprintf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);
1.1 takayama 1296: } /* end of J loop */
1297: }
1298: }
1.30 takayama 1299: } return(0);
1.1 takayama 1300: }
1301:
1.32 takayama 1302: #ifdef STANDALONE
1.1 takayama 1303: /* checkJack1(3,3)
1.12 takayama 1304: */
1.10 takayama 1305: static int checkJack1(int M,int N) {
1.1 takayama 1306: int I,K;
1307: extern int P_pmn;
1308: extern double M_x[];
1309: int Pmn; /* used in aM_jack */
1310: /* initialize x vars. */
1311: for (I=1; I<=N; I++) {
1312: M_x[I-1] = ((double)I)/10.0;
1313: }
1314: genJack(M,N);
1315: Pmn = P_pmn;
1316: for (I=1; I<=N; I++) {
1317: for (K=0; K<=P_pmn; K++) {
1318: printp(Parray[K]);
1.29 takayama 1319: oxprintf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
1.1 takayama 1320: }
1321: }
1.29 takayama 1322: for (I=1; I<=N; I++) oxprintf("%lf, ",M_x[I-1]);
1323: oxprintf("<--x\n");
1.13 takayama 1324: return(0);
1.1 takayama 1325: }
1326: /*main() { checkJack1(3,3); }*/
1327:
1328:
1.10 takayama 1329: static int checkJack2(int M,int N) {
1.1 takayama 1330: int I,K,J;
1331: extern int P_pmn;
1332: extern double M_x[];
1.10 takayama 1333: extern int M_df;
1.1 takayama 1334: int Pmn; /* used in aM_jack */
1335: M_df=1;
1336: /* initialize x vars. */
1337: for (I=1; I<=N; I++) {
1338: M_x[I-1] = ((double)I)/10.0;
1339: }
1340: genJack(M,N);
1341: Pmn = P_pmn;
1342: for (I=1; I<=N; I++) {
1343: for (K=0; K<=P_pmn; K++) {
1344: printp(Parray[K]);
1.29 takayama 1345: oxprintf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
1.1 takayama 1346: }
1347: }
1.29 takayama 1348: for (I=1; I<=N; I++) oxprintf("%lf, ",M_x[I-1]);
1349: oxprintf("<--x\n");
1.1 takayama 1350:
1351: for (I=1; I<=N; I++) {
1352: for (K=0; K<=P_pmn; K++) {
1353: for (J=0; J<M_2n; J++) {
1.12 takayama 1354: printp(Parray[K]);
1.29 takayama 1355: oxprintf("<--Kap, Nv=%d,J(diff)=%d, D^J Jack=%lf\n",
1.12 takayama 1356: I,J,aM_jack(I,J,K));
1357: }
1358: }
1.1 takayama 1359: }
1.13 takayama 1360: return(0);
1.1 takayama 1361: }
1362:
1363: /* main() { checkJack2(3,3); } */
1.32 takayama 1364: #endif
1.1 takayama 1365:
1366: double mh_t(double A[A_LEN],double B[B_LEN],int N,int M) {
1367: double F,F2;
1368: extern int M_df;
1369: extern int P_pmn;
1370: extern double *M_qk;
1371: extern double M_rel_error;
1.15 takayama 1372: extern int M_m;
1.16 takayama 1373: extern int M_m_estimated_approx_deg;
1374: extern double M_assigned_series_error;
1.1 takayama 1375: int Pmn;
1376: int K;
1377: int *Kap;
1378: int size;
1.15 takayama 1379: int i;
1380: double partial_sum[M_m_MAX+1];
1.16 takayama 1381: double iv;
1382: double serror;
1.1 takayama 1383: F = 0; F2=0;
1384: M_df=1;
1385: genJack(M,N);
1.8 takayama 1386: M_qk = (double *)mymalloc(sizeof(double)*(P_pmn+1)); /* found a bug. */
1.1 takayama 1387: Pmn = P_pmn;
1388: size = ParraySize[P_pmn];
1389: for (K=0; K<=P_pmn; K++) {
1.11 takayama 1390: mh_check_intr(100);
1.12 takayama 1391: Kap = Parray[K];
1392: M_qk[K] = qk(Kap,A,B);
1393: F += M_qk[K]*aM_jack(N,0,K);
1394: if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K);
1.29 takayama 1395: if (Debug) oxprintf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size);
1396: if (Debug && (ParraySize[K] == size)) oxprintf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K));
1.1 takayama 1397: }
1398: M_rel_error = F-F2;
1.15 takayama 1399:
1.16 takayama 1400: M_m_estimated_approx_deg = -1; serror=1;
1.15 takayama 1401: for (i=0; i<=M_m; i++) {
1402: partial_sum[i] = 0.0; partial_sum[i+1] = 0.0;
1403: for (K=0; K<=P_pmn; K++) {
1404: if (ParraySize[K] == i) partial_sum[i] += M_qk[K]*aM_jack(N,0,K);
1405: }
1406: if (i>0) partial_sum[i] += partial_sum[i-1];
1.31 takayama 1407: if (i>0) serror = myabs((partial_sum[i]-partial_sum[i-1])/partial_sum[i-1]);
1.16 takayama 1408: if ((i>0)&&(M_m_estimated_approx_deg < 0)&&(serror<M_assigned_series_error)) {
1409: M_m_estimated_approx_deg = i; break;
1410: }
1411: }
1412: if (M_m_estimated_approx_deg < 0) {
1413: M_m_estimated_approx_deg = M_m+mymin(5,mymax(1,(int)log(serror/M_assigned_series_error))); /* Heuristic */
1.15 takayama 1414: }
1415: /*
1416: for (K=0; K<=P_pmn; K++) {
1.29 takayama 1417: oxprintf("Kappa="); for (i=0; i<N; i++) oxprintf("%d ",Parray[K][i]); oxprintf("\n");
1418: oxprintf("ParraySize(%d)=%d (|kappa|), M_m=%d\n",K,ParraySize[K],M_m);
1.15 takayama 1419: }
1420: for (i=0; i<=M_m; i++) {
1.29 takayama 1421: oxprintf("partial_sum[%d]=%lg\n",i,partial_sum[i]);
1.15 takayama 1422: }
1423: */
1.16 takayama 1424: M_estimated_X0g = X0g;
1425: iv=myabs(F*iv_factor());
1426: if (iv < M_x0value_min) M_estimated_X0g = X0g*mymax(2,log(log(1/iv))); /* This is heuristic */
1.19 takayama 1427: M_estimated_X0g = mymin(M_estimated_X0g,M_X0g_bound);
1.16 takayama 1428: M_mh_t_success = 1;
1.19 takayama 1429: if (M_estimated_X0g != X0g) M_mh_t_success=0;
1430: if (M_m_estimated_approx_deg > M_m) M_mh_t_success=0;
1.16 takayama 1431:
1.19 takayama 1432: M_series_error = serror;
1433: M_recommended_abserr = iv*M_assigned_series_error;
1.27 takayama 1434: M_recommended_relerr = M_series_error;
1.25 takayama 1435:
1436: if (M_show_autosteps) {
1.29 takayama 1437: 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);
1438: oxprintf("%%%%x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
1439: oxprintf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
1440: 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);
1441: oxprintfe("%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
1.40 ! takayama 1442: 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);
1.25 takayama 1443: }
1.16 takayama 1444:
1.19 takayama 1445: M_mh_t_value=F;
1.1 takayama 1446: return(F);
1447: }
1448: double mh_t2(int J) {
1449: extern double *M_qk;
1450: double F;
1451: int K;
1452: int Pmn;
1453: extern int P_pmn;
1.3 takayama 1454: if (M_qk == NULL) {myerror("Call mh_t first."); mh_exit(-1); }
1.1 takayama 1455: F = 0;
1456: Pmn = P_pmn;
1457: for (K=0; K<P_pmn; K++) {
1.12 takayama 1458: F += M_qk[K]*aM_jack(M_n,J,K);
1.1 takayama 1459: }
1460: return(F);
1461: }
1462:
1.32 takayama 1463: #ifdef STANDALONE
1.10 takayama 1464: static int mtest1b() {
1.1 takayama 1465: double A[1] = {1.5};
1466: double B[1] = {1.5+5};
1467: int I,N,M,J;
1468: double F;
1469: N=3; M=6;
1470: for (I=1; I<=N; I++) {
1471: M_x[I-1] = ((double)I)/10.0;
1472: }
1473: mh_t(A,B,N,M);
1474: for (J=0; J<M_2n; J++) {
1.12 takayama 1475: F=mh_t2(J);
1.29 takayama 1476: oxprintf("J=%d, D^J mh_t=%lf\n",J,F);
1.1 takayama 1477: }
1.30 takayama 1478: return(0);
1.1 takayama 1479: }
1480:
1481: /* main() { mtest1b(); }*/
1.32 takayama 1482: #endif
1.1 takayama 1483:
1484:
1485:
1486: #define TEST 1
1487: #ifndef TEST
1488:
1489: #endif
1490:
1491: /****** from mh-n.c *****/
1492:
1493: #define SMAX 4096
1.29 takayama 1494: #define inci(i) { i++; if (i >= argc) { oxprintfe("Option argument is not given.\n"); return(NULL); }}
1.1 takayama 1495:
1.2 takayama 1496: static int imypower(int x,int n) {
1.1 takayama 1497: int i;
1498: int v;
1.3 takayama 1499: if (n < 0) {myerror("imypower"); mh_exit(-1);}
1.1 takayama 1500: v = 1;
1501: for (i=0; i<n; i++) v *= x;
1502: return(v);
1503: }
1504:
1.11 takayama 1505: #ifdef STANDALONE2
1.32 takayama 1506: int main(int argc,char *argv[]) {
1.3 takayama 1507: mh_exit(MH_RESET_EXIT);
1.12 takayama 1508: /* jk_main(argc,argv);
1.29 takayama 1509: oxprintf("second run.\n"); */
1.2 takayama 1510: jk_main(argc,argv);
1.17 takayama 1511: return(0);
1.2 takayama 1512: }
1.3 takayama 1513: #endif
1514:
1515: struct MH_RESULT *jk_main(int argc,char *argv[]) {
1.16 takayama 1516: int i;
1517: struct MH_RESULT *ans;
1518: extern int M_automatic;
1519: extern int M_mh_t_success;
1520: extern double M_estimated_X0g;
1521: extern int M_m_estimated_approx_deg;
1522: for (i=1; i<argc; i++) {
1523: if (strcmp(argv[i],"--automatic")==0) {
1524: inci(i);
1525: sscanf(argv[i],"%d",&M_automatic);
1526: break;
1527: }
1528: }
1529: ans=jk_main2(argc,argv,0,0.0,0);
1530: if (!M_automatic) return(ans);
1531: if (M_mh_t_success) return(ans);
1532: while (!M_mh_t_success) {
1533: ans=jk_main2(argc,argv,1,M_estimated_X0g,M_m_estimated_approx_deg);
1534: }
1535: return(ans);
1536: }
1537:
1538: struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree) {
1.32 takayama 1539: // double *y0;
1540: // double x0,xn;
1541: // double ef;
1542:
1543: int i,j; // int i,j,rank;
1.35 takayama 1544: double a[A_LEN]; double b[B_LEN];
1.1 takayama 1545: extern double M_x[];
1546: extern double *Beta;
1547: extern int M_2n;
1.16 takayama 1548: extern int M_mh_t_success;
1.3 takayama 1549: char swork[1024];
1550: struct MH_RESULT *ans=NULL;
1551: struct SFILE *ofp = NULL;
1552: int idata=0;
1553: JK_byFile = 1;
1554: jk_initializeWorkArea();
1.1 takayama 1555: UseTable = 1;
1.3 takayama 1556: Mapprox=6;
1.1 takayama 1557: for (i=1; i<argc; i++) {
1.12 takayama 1558: if (strcmp(argv[i],"--idata")==0) {
1559: inci(i);
1560: setParam(argv[i]); idata=1;
1561: }else if (strcmp(argv[i],"--degree")==0) {
1562: inci(i);
1563: sscanf(argv[i],"%d",&Mapprox);
1564: }else if (strcmp(argv[i],"--x0")==0) {
1565: inci(i);
1566: sscanf(argv[i],"%lg",&X0g);
1567: }else if (strcmp(argv[i],"--notable")==0) {
1568: UseTable = 0;
1569: }else if (strcmp(argv[i],"--debug")==0) {
1570: Debug = 1;
1571: }else if (strcmp(argv[i],"--help")==0) {
1572: usage(); return(0);
1573: }else if (strcmp(argv[i],"--bystring")==0) {
1.29 takayama 1574: if (idata) {oxprintfe("--bystring must come before --idata option.\n"); mh_exit(-1);}
1.12 takayama 1575: JK_byFile = 0;
1.16 takayama 1576: }else if (strcmp(argv[i],"--automatic")==0) {
1577: inci(i); /* ignore, in this function */
1.23 takayama 1578: }else if (strcmp(argv[i],"--assigned_series_error")==0) {
1.16 takayama 1579: inci(i);
1580: sscanf(argv[i],"%lg",&M_assigned_series_error);
1581: }else if (strcmp(argv[i],"--x0value_min")==0) {
1582: inci(i);
1583: sscanf(argv[i],"%lg",&M_x0value_min);
1.12 takayama 1584: }else {
1.29 takayama 1585: oxprintfe("Unknown option %s\n",argv[i]);
1.12 takayama 1586: usage();
1587: return(NULL);
1588: }
1.1 takayama 1589: }
1.16 takayama 1590: if (!idata) setParam(NULL);
1591: if (automode) {
1592: Mapprox = newDegree;
1593: X0g = newX0g;
1594: }
1.3 takayama 1595:
1596: /* Initialize global variables */
1597: M_n = Mg;
1598: HS_n=M_n;
1.16 takayama 1599: if (!JK_byFile) {
1600: ans = (struct MH_RESULT *)mymalloc(sizeof(struct MH_RESULT));
1601: ans->message = NULL;
1.19 takayama 1602: ans->t_success = 0;
1603: ans->series_error = 1.0e+10;
1604: ans->recommended_abserr = 1.0e-10;
1.16 takayama 1605: }
1.3 takayama 1606: else ans = NULL;
1.26 takayama 1607: if (M_automatic) {
1608: /* Differentiation can be M_m in the bit pattern in the M_n variable case.*/
1609: if (M_n > Mapprox) Mapprox=M_n;
1610: }
1.3 takayama 1611: /* Output by a file=stdout */
1.34 takayama 1612: ofp = mh_fopen("stdout","w",JK_byFile);
1.1 takayama 1613:
1.3 takayama 1614: sprintf(swork,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp);
1615: sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp);
1.1 takayama 1616: if (M_n != Mg) {
1.12 takayama 1617: myerror("Mg must be equal to M_n\n"); mh_exit(-1);
1.1 takayama 1618: }
1.8 takayama 1619: if (Debug) showParam(NULL,1);
1.40 ! takayama 1620: setM_x();
1.19 takayama 1621:
1622: M_beta_i_x_o2_max=myabs(M_x[0]/2);
1623: if (M_n <= 1) M_beta_i_beta_j_min = myabs(Beta[0]);
1624: else M_beta_i_beta_j_min = myabs(Beta[1]-Beta[0]);
1625: for (i=0; i<M_n; i++) {
1626: if (myabs(M_x[i]/2) > M_beta_i_x_o2_max) M_beta_i_x_o2_max = myabs(M_x[i]/2);
1627: for (j=i+1; j<M_n; j++) {
1628: if (myabs(Beta[i]-Beta[j]) < M_beta_i_beta_j_min)
1629: M_beta_i_beta_j_min = myabs(Beta[i]-Beta[j]);
1630: }
1631: }
1632:
1.36 takayama 1633: if ((P_pFq != A_LEN) || (Q_pFq != B_LEN)) {
1634: oxprintfe("It must be P_pFq == A_LEN and Q_pFq == B_LEN in this version. %s\n","");
1635: mh_exit(-1);
1636: }
1637: oxprintfe("%%%%(stderr) Orig_1F1=%d\n",Orig_1F1);
1638: if ((P_pFq == 1) && (Q_pFq == 1) && (Orig_1F1)) {
1639: A_pFq[0] = a[0] = ((double)Mg+1.0)/2.0;
1640: B_pFq[0] = b[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */
1641: if (Debug) oxprintf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",a[0],b[0],M_n,Mapprox);
1642: }else{
1643: for (i=0; i<P_pFq; i++) a[i] = A_pFq[i];
1644: for (i=0; i<Q_pFq; i++) b[i] = B_pFq[i];
1645: }
1.1 takayama 1646: mh_t(a,b,M_n,Mapprox);
1.19 takayama 1647: if ((!M_mh_t_success) && M_automatic) {
1.16 takayama 1648: jk_freeWorkArea();
1649: return NULL;
1650: }
1.1 takayama 1651: if (imypower(2,M_n) != M_2n) {
1.12 takayama 1652: sprintf(swork,"M_n=%d,M_2n=%d\n",M_n,M_2n); mh_fputs(swork,ofp);
1653: myerror("2^M_n != M_2n\n"); mh_exit(-1);
1.1 takayama 1654: }
1.3 takayama 1655: sprintf(swork,"%%%%M_rel_error=%lg\n",M_rel_error); mh_fputs(swork,ofp);
1.1 takayama 1656: for (j=0; j<M_2n; j++) {
1.12 takayama 1657: Iv[j] = mh_t2(j);
1.1 takayama 1658: }
1659: Ef = iv_factor();
1.8 takayama 1660: showParam(ofp,0);
1.3 takayama 1661:
1662: /* return the result */
1.16 takayama 1663: if (!JK_byFile) {
1.12 takayama 1664: ans->x = X0g;
1665: ans->rank = imypower(2,Mg);
1666: ans->y = (double *)mymalloc(sizeof(double)*(ans->rank));
1667: for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i];
1668: ans->size=1;
1669: ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size));
1670: (ans->sfpp)[0] = ofp;
1.19 takayama 1671:
1672: ans->t_success = M_mh_t_success;
1673: ans->series_error = M_series_error;
1674: ans->recommended_abserr = M_recommended_abserr;
1.3 takayama 1675: }
1.29 takayama 1676: if (Debug) oxprintf("jk_freeWorkArea() starts\n");
1.3 takayama 1677: jk_freeWorkArea();
1.29 takayama 1678: if (Debug) oxprintf("jk_freeWorkArea() has finished.\n");
1.3 takayama 1679: return(ans);
1.1 takayama 1680: }
1681:
1.10 takayama 1682: static int usage() {
1.29 takayama 1683: oxprintfe("Usages:\n");
1.37 takayama 1684: #ifdef C_2F1
1685: oxprintfe("hgm_jack-n-2f1");
1686: #else
1687: oxprintfe("hgm_jack-n ");
1688: #endif
1689: oxprintfe(" [--idata input_data_file --x0 x0 --degree approxm]\n");
1.29 takayama 1690: oxprintfe(" [--automatic n --assigned_series_error e --x0value_min e2]\n");
1691: oxprintfe("\nThe command hgm_jack-n [options] generates an input for hgm_w-n, Pr({y | y<xmax}), which is the cumulative distribution function of the largest root of the m by m Wishart matrices with n degrees of freedom and the covariantce matrix sigma.\n");
1692: oxprintfe("The hgm_jack-n uses the Koev-Edelman algorithm to evalute the matrix hypergeometric function.\n");
1693: oxprintfe("The degree of the approximation (Mapprox) is given by the --degree option.\n");
1694: oxprintfe("Parameters are specified by the input_data_file. Otherwise, default values are used.\n\n");
1695: oxprintfe("The format of the input_data_file: (The orders of the input data must be kept.)\n");
1696: oxprintfe(" Mg: m(the number of variables), Beta: beta=sigma^(-1)/2 (diagonized), Ng: n,\n");
1697: oxprintfe(" (Add a comment line %%Ng= before the data Ng to check the number of beta.)\n");
1698: oxprintfe(" X0g: starting value of x(when --x0 option is used, this value is used)\n");
1699: oxprintfe(" 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");
1700: oxprintfe(" Ef: a scalar factor to the initial value. It is calculated by this program. Give the zero.\n");
1701: oxprintfe(" Hg: h (step size) which is for hgm_w-n, \n");
1702: oxprintfe(" Dp: output data is stored in every Dp steps when output_data_file is specified. This is for hgm_w-n.\n");
1703: oxprintfe(" Xng: terminating value of x. This is for hgm_w-n.\n");
1704: oxprintfe("Optional parameters automatic, ... are interpreted by a parser. See setParam() in jack-n.c and Testdata/tmp-idata2.txt as an example. Optional paramters are given as %%parameter_name=value Lines starting with %%%% or # are comment lines.\n");
1705: oxprintfe("Parameters are redefined when they appear more than once in the idata file and the command line options.\n\n");
1706: oxprintfe("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");
1707: oxprintfe("An example format of the input_data_file can be obtained by executing hgm_jack-n with no option. When there is no --idata file, all options are ignored.\n");
1708: oxprintfe("By --automatic option, X0g and degree are automatically determined from assigend_series_error. The current strategy is described in mh_t in jack-n.c\n");
1709: oxprintfe("Default values for the papameters of the automatic mode: automatic=%d, assigned_series_error=%lg, x0value_min=%lg\n",M_automatic,M_assigned_series_error,M_x0value_min);
1.37 takayama 1710: #ifdef C_2F1
1711: oxprintfe("The parameters a,b,c of 2F1 are given by %%p_pFq=2, a,b and %%q_pFq=1, c\nNg is ignored.\n");
1712: #endif
1.29 takayama 1713: oxprintfe("Todo: automatic mode throws away the table of Jack polynomials of the previous degrees and reevaluate them. They should be kept.\n");
1714: oxprintfe("\nExamples:\n");
1715: oxprintfe("[1] ./hgm_jack-n \n");
1716: oxprintfe("[2] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --degree 15 --automatic 0\n");
1717: oxprintfe("[3] ./hgm_jack-n --idata Testdata/tmp-idata2.txt --degree 15 >test2.txt\n");
1718: oxprintfe(" ./hgm_w-n --idata test2.txt --gnuplotf test-g\n");
1719: oxprintfe(" gnuplot -persist <test-g-gp.txt\n");
1720: oxprintfe("[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --automatic 1 --assigned_series_error=1e-12\n");
1721: oxprintfe("[5] ./hgm_jack-n --idata Testdata/tmp-idata4.txt\n");
1.37 takayama 1722: #ifdef C_2F1
1723: oxprintfe("Todo for 2F1: example for hgm_jack-n-2f1 has not been written.\niv_factor? Ef?");
1724: #endif
1.13 takayama 1725: return(0);
1.1 takayama 1726: }
1727:
1.10 takayama 1728: static int setParamDefault() {
1.1 takayama 1729: int rank;
1730: int i;
1.3 takayama 1731: Mg = M_n_default ;
1.1 takayama 1732: rank = imypower(2,Mg);
1733: Beta = (double *)mymalloc(sizeof(double)*Mg);
1734: for (i=0; i<Mg; i++) Beta[i] = 1.0+i;
1735: Ng = (double *)mymalloc(sizeof(double)); *Ng = 3.0;
1736: Iv = (double *)mymalloc(sizeof(double)*rank);
1737: Iv2 = (double *)mymalloc(sizeof(double)*rank);
1738: for (i=0; i<rank; i++) Iv[i] = 0;
1739: Ef = 0;
1740: Ef2 = 0.01034957388338225707;
1741: if (M_n == 2) {
1742: Iv2[0] = 1.58693;
1743: Iv2[1] = 0.811369;
1744: Iv2[2] = 0.846874;
1745: Iv2[3] = 0.413438;
1746: }
1747: X0g = (Beta[0]/Beta[Mg-1])*0.5;
1748: Hg = 0.001;
1749: Dp = 1;
1.38 takayama 1750: if ((P_pFq == 1) && (Q_pFq == 1)) {
1751: Xng = 10.0;
1752: }else {
1753: Xng=0.25;
1754: for (i=0; i<A_LEN; i++) A_pFq[i] = (i+1)/5.0;
1755: for (i=0; i<B_LEN; i++) B_pFq[i] = (A_LEN+i+1)/7.0;
1756: }
1.13 takayama 1757: return(0);
1.1 takayama 1758: }
1759:
1.10 takayama 1760: static int next(struct SFILE *sfp,char *s,char *msg) {
1.14 takayama 1761: static int check=1;
1.31 takayama 1762: char *ng="%%Ng=";
1.32 takayama 1763: // int i;
1.1 takayama 1764: s[0] = '%';
1765: while (s[0] == '%') {
1.12 takayama 1766: if (!mh_fgets(s,SMAX,sfp)) {
1.29 takayama 1767: oxprintfe("Data format error at %s\n",msg);
1.12 takayama 1768: mh_exit(-1);
1769: }
1.14 takayama 1770: if (check && (strncmp(msg,ng,4)==0)) {
1.31 takayama 1771: if (strncmp(s,ng,5) != 0) {
1.29 takayama 1772: oxprintfe("Warning, there is no %%Ng= at the border of Beta's and Ng, s=%s\n",s);
1.14 takayama 1773: }
1.31 takayama 1774: /* check=0; */
1.14 takayama 1775: }
1.12 takayama 1776: if (s[0] != '%') return(0);
1.1 takayama 1777: }
1.13 takayama 1778: return(0);
1.1 takayama 1779: }
1.10 takayama 1780: static int setParam(char *fname) {
1.1 takayama 1781: int rank;
1782: char s[SMAX];
1.3 takayama 1783: struct SFILE *fp;
1.1 takayama 1784: int i;
1.19 takayama 1785: struct mh_token tk;
1.1 takayama 1786: if (fname == NULL) return(setParamDefault());
1787:
1788: Sample = 0;
1.3 takayama 1789: if ((fp=mh_fopen(fname,"r",JK_byFile)) == NULL) {
1.29 takayama 1790: if (JK_byFile) oxprintfe("File %s is not found.\n",fp->s);
1.12 takayama 1791: mh_exit(-1);
1.1 takayama 1792: }
1793: next(fp,s,"Mg(m)");
1794: sscanf(s,"%d",&Mg);
1795: rank = imypower(2,Mg);
1796:
1797: Beta = (double *)mymalloc(sizeof(double)*Mg);
1798: for (i=0; i<Mg; i++) {
1799: next(fp,s,"Beta");
1.12 takayama 1800: sscanf(s,"%lf",&(Beta[i]));
1.1 takayama 1801: }
1802:
1803: Ng = (double *)mymalloc(sizeof(double));
1.14 takayama 1804: next(fp,s,"%Ng= (freedom parameter n)");
1.1 takayama 1805: sscanf(s,"%lf",Ng);
1806:
1807: next(fp,s,"X0g(initial point)");
1808: sscanf(s,"%lf",&X0g);
1.14 takayama 1809:
1.1 takayama 1810: Iv = (double *)mymalloc(sizeof(double)*rank);
1811: for (i=0; i<rank; i++) {
1.12 takayama 1812: next(fp,s,"Iv(initial values)");
1.14 takayama 1813: if (strncmp(s,"*",1)==0) {
1814: for (i=0; i<rank; i++) Iv[i] = 0.0;
1815: break;
1816: }
1.12 takayama 1817: sscanf(s,"%lg",&(Iv[i]));
1.1 takayama 1818: }
1819:
1820: next(fp,s,"Ef(exponential factor)");
1.14 takayama 1821: if (strncmp(s,"*",1)==0) Ef=0.0;
1822: else sscanf(s,"%lg",&Ef);
1.1 takayama 1823:
1824: next(fp,s,"Hg (step size of rk)");
1825: sscanf(s,"%lf",&Hg);
1826:
1827: next(fp,s,"Dp (data sampling period)");
1828: sscanf(s,"%d",&Dp);
1829:
1830: next(fp,s,"Xng (the last point, cf. --xmax)");
1831: sscanf(s,"%lf",&Xng);
1.19 takayama 1832:
1833: /* Reading the optional parameters */
1834: while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) {
1835: /* expect ID */
1836: if (tk.type != MH_TOKEN_ID) {
1.29 takayama 1837: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1838: }
1839: if (strcmp(s,"automatic")==0) {
1840: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1.29 takayama 1841: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1842: }
1843: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1.29 takayama 1844: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1845: }
1846: M_automatic = tk.ival;
1847: continue;
1848: }
1849: if (strcmp(s,"assigned_series_error")==0) {
1850: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1.29 takayama 1851: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1852: }
1853: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
1.29 takayama 1854: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1855: }
1856: M_assigned_series_error = tk.dval;
1857: continue;
1858: }
1859: if (strcmp(s,"x0value_min")==0) {
1860: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1.29 takayama 1861: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1862: }
1863: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
1.29 takayama 1864: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1865: }
1866: M_x0value_min = tk.dval;
1867: continue;
1868: }
1869: if ((strcmp(s,"Mapprox")==0) || (strcmp(s,"degree")==0)) {
1870: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1.29 takayama 1871: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1872: }
1873: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1.29 takayama 1874: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1875: }
1876: Mapprox = tk.ival;
1877: continue;
1878: }
1879: if (strcmp(s,"X0g_bound")==0) {
1880: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1.29 takayama 1881: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1882: }
1883: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
1.29 takayama 1884: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.19 takayama 1885: }
1886: M_X0g_bound = tk.dval;
1887: continue;
1888: }
1.25 takayama 1889: if (strcmp(s,"show_autosteps")==0) {
1890: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1.29 takayama 1891: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.25 takayama 1892: }
1893: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1.29 takayama 1894: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1.25 takayama 1895: }
1896: M_show_autosteps = tk.ival;
1897: continue;
1898: }
1.36 takayama 1899: // Format: #p_pFq=2 1.5 3.2
1900: if (strcmp(s,"p_pFq")==0) {
1901: Orig_1F1=0;
1902: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1903: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1904: }
1905: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1906: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1907: }
1908: P_pFq = tk.ival;
1909: for (i=0; i<P_pFq; i++) {
1910: if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
1911: A_pFq[i] = tk.dval;
1912: }else if (tk.type == MH_TOKEN_INT) {
1913: A_pFq[i] = tk.ival;
1914: }else{
1915: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1916: }
1917: }
1918: continue;
1919: }
1920: if (strcmp(s,"q_pFq")==0) {
1921: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1922: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1923: }
1924: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1925: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1926: }
1927: Q_pFq = tk.ival;
1928: for (i=0; i<Q_pFq; i++) {
1929: if ((tk=mh_getoken(s,SMAX-1,fp)).type == MH_TOKEN_DOUBLE) {
1930: B_pFq[i] = tk.dval;
1931: }else if (tk.type == MH_TOKEN_INT) {
1932: B_pFq[i] = tk.ival;
1933: }else{
1934: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
1935: }
1936: }
1937: continue;
1938: }
1.40 ! takayama 1939: if (strcmp(s,"ef_type")==0) {
! 1940: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
! 1941: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
! 1942: }
! 1943: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
! 1944: oxprintfe("Syntax error at %s\n",s); mh_exit(-1);
! 1945: }
! 1946: Ef_type = tk.ival;
! 1947: continue;
! 1948: }
1.29 takayama 1949: oxprintfe("Unknown ID at %s\n",s); mh_exit(-1);
1.19 takayama 1950: }
1.3 takayama 1951: mh_fclose(fp);
1.13 takayama 1952: return(0);
1.1 takayama 1953: }
1954:
1.10 takayama 1955: static int showParam(struct SFILE *fp,int fd) {
1.1 takayama 1956: int rank,i;
1.3 takayama 1957: char swork[1024];
1.8 takayama 1958: if (fd) {
1.34 takayama 1959: fp = mh_fopen("stdout","w",1);
1.8 takayama 1960: }
1.1 takayama 1961: rank = imypower(2,Mg);
1.3 takayama 1962: sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp);
1.1 takayama 1963: for (i=0; i<Mg; i++) {
1.12 takayama 1964: sprintf(swork,"%%Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp);
1.1 takayama 1965: }
1.3 takayama 1966: sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp);
1967: sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp);
1.1 takayama 1968: for (i=0; i<rank; i++) {
1.12 takayama 1969: sprintf(swork,"%%Iv[%d]=\n%lg\n",i,Iv[i]); mh_fputs(swork,fp);
1970: if (Sample && (M_n == 2) && (X0g == 0.3)) {
1971: sprintf(swork,"%%Iv[%d]-Iv2[%d]=%lg\n",i,i,Iv[i]-Iv2[i]); mh_fputs(swork,fp);
1972: }
1.3 takayama 1973: }
1974: sprintf(swork,"%%Ef=\n%lg\n",Ef); mh_fputs(swork,fp);
1975: sprintf(swork,"%%Hg=\n%lf\n",Hg); mh_fputs(swork,fp);
1976: sprintf(swork,"%%Dp=\n%d\n",Dp); mh_fputs(swork,fp);
1977: sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp);
1.19 takayama 1978:
1979: sprintf(swork,"%%%% Optional paramters\n"); mh_fputs(swork,fp);
1980: sprintf(swork,"#success=%d\n",M_mh_t_success); mh_fputs(swork,fp);
1981: sprintf(swork,"#automatic=%d\n",M_automatic); mh_fputs(swork,fp);
1982: sprintf(swork,"#series_error=%lg\n",M_series_error); mh_fputs(swork,fp);
1983: sprintf(swork,"#recommended_abserr\n"); mh_fputs(swork,fp);
1.36 takayama 1984: sprintf(swork,"#abserror=%lg\n",M_recommended_abserr); mh_fputs(swork,fp);
1.27 takayama 1985: if (M_recommended_relerr < MH_RELERR_DEFAULT) {
1986: sprintf(swork,"%%relerror=%lg\n",M_recommended_relerr); mh_fputs(swork,fp);
1987: }
1.19 takayama 1988: sprintf(swork,"#mh_t_value=%lg # Value of matrix hg at X0g.\n",M_mh_t_value); mh_fputs(swork,fp);
1989: sprintf(swork,"# M_m=%d # Approximation degree of matrix hg.\n",M_m); mh_fputs(swork,fp);
1990: sprintf(swork,"#beta_i_x_o2_max=%lg #max(|beta[i]*x|/2)\n",M_beta_i_x_o2_max); mh_fputs(swork,fp);
1991: sprintf(swork,"#beta_i_beta_j_min=%lg #min(|beta[i]-beta[j]|)\n",M_beta_i_beta_j_min); mh_fputs(swork,fp);
1.36 takayama 1992: sprintf(swork,"# change # to %% to read as an optional parameter.\n"); mh_fputs(swork,fp);
1.39 takayama 1993: sprintf(swork,"%%p_pFq=%d, ",P_pFq); mh_fputs(swork,fp);
1.36 takayama 1994: for (i=0; i<P_pFq; i++) {
1995: if (i != P_pFq-1) sprintf(swork," %lg,",A_pFq[i]);
1996: else sprintf(swork," %lg\n",A_pFq[i]);
1997: mh_fputs(swork,fp);
1998: }
1.39 takayama 1999: sprintf(swork,"%%q_pFq=%d, ",Q_pFq); mh_fputs(swork,fp);
1.36 takayama 2000: for (i=0; i<Q_pFq; i++) {
2001: if (i != Q_pFq-1) sprintf(swork," %lg,",B_pFq[i]);
2002: else sprintf(swork," %lg\n",B_pFq[i]);
2003: mh_fputs(swork,fp);
2004: }
1.40 ! takayama 2005: sprintf(swork,"%%ef_type=%d\n",Ef_type); mh_fputs(swork,fp);
1.13 takayama 2006: return(0);
1.1 takayama 2007: }
2008:
1.2 takayama 2009: static double gammam(double a,int n) {
1.1 takayama 2010: double v,v2;
2011: int i;
1.39 takayama 2012: v=mypower(sqrt(M_PI),(n*(n-1))/2); /* pi^(n*(n-1)/2) */
1.1 takayama 2013: v2=0;
2014: for (i=1; i<=n; i++) {
2015: v2 += lgamma(a-((double)(i-1))/2.0); /* not for big n */
2016: }
1.29 takayama 2017: if (Debug) oxprintf("gammam(%lg,%d)=%lg\n",a,n,v*exp(v2));
1.1 takayama 2018: return(v*exp(v2));
2019: }
2020:
1.2 takayama 2021: static double iv_factor(void) {
1.40 ! takayama 2022: if (Ef_type < 1) return(1.0);
! 2023: if (Ef_type == 1) return(iv_factor_ef_type1());
! 2024: else if (Ef_type==2) return(iv_factor_ef_type2());
! 2025: return(1.0);
! 2026: }
! 2027:
! 2028: static double iv_factor_ef_type1(void) {
1.1 takayama 2029: double v1;
2030: double t;
2031: double b;
2032: double detSigma;
2033: double c;
2034: int i,n;
2035: n = (int) (*Ng);
2036: v1= mypower(sqrt(X0g),n*M_n);
2037: t = 0.0;
2038: for (i=0; i<M_n; i++) t += -X0g*Beta[i];
2039: v1 = v1*exp(t);
2040:
2041: b = 1; for (i=0; i<M_n; i++) b *= Beta[i];
2042: detSigma = 1.0/(b*mypower(2.0,M_n));
2043:
2044: c = gammam(((double)(M_n+1))/2.0, M_n)/
1.12 takayama 2045: ( mypower(sqrt(2), M_n*n)*mypower(sqrt(detSigma),n)*gammam(((double)(M_n+n+1))/2.0,M_n) );
1.1 takayama 2046: return( c*v1);
2047: }
2048:
1.40 ! takayama 2049: static double iv_factor_ef_type2(void) {
! 2050: double ef;
! 2051: int i,m;
! 2052: double a,b,c;
! 2053: double t;
! 2054: a = A_pFq[0]; b = A_pFq[1]; c= B_pFq[0];
! 2055: m = M_n;
! 2056: ef = 1.0;
! 2057:
! 2058: /*fprintf(stderr,"iv_factor_ef_type2: a=%lf, b=%lf, c=%lf, m=%d\n",a,b,c,m);*/
! 2059: /* Ref: note 2016.02.04 */
! 2060: if (X0g<0){ myerror("Negative X0g\n"); mh_exit(-1);}
! 2061: t = 0;
! 2062: for (i=0; i<m; i++) if (Beta[i]<0){ myerror("Negative beta\n"); mh_exit(-1);}
! 2063: for (i=0; i<m; i++) t += log(Beta[i]);
! 2064: ef *= exp((a+b-c)*t);
! 2065:
! 2066: t = 0;
! 2067: for (i=0; i<m; i++) t += log(Beta[i]+X0g);
! 2068: ef *= exp(-b*t);
! 2069:
! 2070: ef *= exp((c-a)*(2*a-1)*log(X0g));
! 2071:
! 2072: ef *= gammam(b,m)/gammam(a+b-c,m);
! 2073: ef *= gammam(a,m)/gammam(c,m);
! 2074: return(ef);
! 2075: }
! 2076: static void setM_x(void) {
! 2077: if (Ef_type <= 1) return(setM_x_ef_type1());
! 2078: else if (Ef_type==2) return(setM_x_ef_type2());
! 2079: setM_x_ef_type1();
! 2080: }
! 2081: static void setM_x_ef_type1(void) {
! 2082: int i;
! 2083: for (i=0; i<M_n; i++) {
! 2084: M_x[i] = Beta[i]*X0g;
! 2085: }
! 2086: }
! 2087: static void setM_x_ef_type2(void) {
! 2088: int i;
! 2089: for (i=0; i<M_n; i++) {
! 2090: M_x[i] = X0g/(Beta[i]+X0g);
! 2091: }
! 2092: }
1.1 takayama 2093:
2094:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to jack-n.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / hgm / mh / src