Annotation of OpenXM/src/hgm/mh/src/jack-n.c, Revision 1.13
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.13 ! takayama 8: $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.12 2013/03/07 03:00:43 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:
18: 2012.02.21, porting to OpenXM/src/hgm
19: 2011.12.22, --table option, which is experimental.
20: 2011.12.19, bug fix. jjk() should return double. It can become more than max int.
21: 2011.12.15, mh.r --> jack-n.c
1.1 takayama 22: */
23:
1.3 takayama 24: /****** from mh-n.c *****/
25: static int JK_byFile=1;
1.6 takayama 26: static int JK_deallocate=0;
1.3 takayama 27: #define M_n_default 3
1.4 takayama 28: #define Sample_default 1
1.3 takayama 29: static int M_n=0;
1.1 takayama 30: /* global variables. They are set in setParam() */
1.2 takayama 31: static int Mg; /* n */
32: static int Mapprox; /* m, approximation degree */
33: static double *Beta; /* beta[0], ..., beta[m-1] */
34: static double *Ng; /* freedom n. c=(m+1)/2+n/2; Note that it is a pointer */
35: static double X0g; /* initial point */
36: static double *Iv; /* Initial values of mhg sorted by mhbase() in rd.rr at beta*x0 */
37: static double Ef; /* exponential factor at beta*x0 */
38: static double Hg; /* step size of rk defined in rk.c */
39: static int Dp; /* Data sampling period */
40: static double Xng=0.0; /* the last point */
1.4 takayama 41: static int Sample = Sample_default;
1.1 takayama 42:
43: /* for sample inputs */
1.2 takayama 44: static double *Iv2;
45: static double Ef2;
1.1 takayama 46:
47: #ifdef NAN
48: #else
1.2 takayama 49: #define NAN 3.40282e+38 /* for old 32 bit machines. Todo, configure */
1.1 takayama 50: #endif
51:
52: /* #define M_n 3 defined in the Makefile */ /* number of variables */
53: #define M_n0 3 /* used for tests. Must be equal to M_n */
54: #define M_m_MAX 200
1.3 takayama 55: #define M_nmx M_m_MAX /* maximal of M_n */
1.1 takayama 56: #define A_LEN 1 /* (a_1) , (a_1, ..., a_p)*/
57: #define B_LEN 1 /* (b_1) */
1.9 takayama 58: static int Debug = 0;
1.5 takayama 59: static int Alpha = 2; /* 2 implies the zonal polynomial */
1.2 takayama 60: static int *Darray = NULL;
61: static int **Parray = NULL; /* array of partitions of size M_n */
62: static int *ParraySize = NULL; /* length of each partitions */
1.3 takayama 63: static int M_kap[M_nmx];
1.2 takayama 64: static int M_m=M_m_MAX-2; /* | | <= M_m, bug do check of M_m <=M_m_MAX-2 */
1.1 takayama 65: void (*M_pExec)(void);
1.3 takayama 66: static int HS_mu[M_nmx];
67: static int HS_n=M_nmx; /* It is initialized to M_n in jk_main */
1.1 takayama 68: void (*HS_hsExec)(void);
1.3 takayama 69: static double M_x[M_nmx];
1.1 takayama 70:
71: /* They are used in pmn */
1.2 takayama 72: static int *P_pki=NULL;
73: static int P_pmn=0;
1.1 takayama 74:
75: /* It is used genDarray2(), list partitions... */
1.2 takayama 76: static int DR_parray=0;
1.1 takayama 77:
78: /* Used in genBeta() and enumeration of horizontal strip. */
1.2 takayama 79: static double *M_beta_0=NULL; /* M_beta[0][*] value of beta_{kappa,mu}, M_beta[1][*] N_mu */
80: static int *M_beta_1=NULL;
81: static int M_beta_pt=0;
1.3 takayama 82: static int M_beta_kap[M_nmx];
1.2 takayama 83: static int UseTable = 0;
84:
85: static double **M_jack;
86: static int M_df=1; /* Compute differentials? */
87: static int M_2n=0; /* 2^N */
1.1 takayama 88:
1.3 takayama 89: static double Xarray[M_nmx][M_m_MAX];
1.1 takayama 90: /* (x_1, ..., x_n) */
91: /* Xarray[i][0] x_{i+1}^0, Xarray[i][1], x_{i+1}^1, ... */
92:
1.2 takayama 93: static double *M_qk=NULL; /* saves pochhammerb */
94: static double M_rel_error=0.0; /* relative errors */
1.1 takayama 95:
1.2 takayama 96: /* prototypes */
97: static void *mymalloc(int size);
1.10 takayama 98: static int myfree(void *p);
99: static int myerror(char *s);
1.2 takayama 100: static double jack1(int K);
101: static double jack1diff(int k);
102: static double xval(int i,int p); /* x_i^p */
103: static int mysum(int L[]);
104: static int plength(int P[]);
105: static int plength_t(int P[]);
1.3 takayama 106: static void ptrans(int P[M_nmx],int Pt[]);
1.10 takayama 107: static int test_ptrans();
1.2 takayama 108: static int huk(int K[],int I,int J);
109: static int hdk(int K[],int I,int J);
110: static double jjk(int K[]);
111: static double ppoch(double A,int K[]);
112: static double ppoch2(double A,double B,int K[]);
113: static double mypower(double x,int n);
114: static double qk(int K[],double A[A_LEN],double B[B_LEN]);
115: static int bb(int N[],int K[],int M[],int I,int J);
116: static double beta(int K[],int M[]);
1.10 takayama 117: static int printp(int kappa[]);
118: static int printp2(int kappa[]);
119: static int test_beta();
1.2 takayama 120: static double q3_10(int K[],int M[],int SK);
121: static double q3_5(double A[],double B[],int K[],int I);
122: static void mtest4();
123: static void mtest4b();
124: static int nk(int KK[]);
125: static int plength2(int P1[],int P2[]);
126: static int myeq(int P1[],int P2[]);
127: static int pListPartition(int M,int N);
128: static int pListPartition2(int Less,int From,int To, int M);
129: static void pExec_0();
130: static int pListHS(int Kap[],int N);
131: static int pListHS2(int From,int To,int Kap[]);
132: static void hsExec_0();
133: static int pmn(int M,int N);
134: static int *cloneP(int a[]);
1.10 takayama 135: static int copyP(int p[],int a[]);
1.2 takayama 136: static void pExec_darray(void);
1.10 takayama 137: static int genDarray2(int M,int N);
138: static int isHStrip(int Kap[],int Nu[]);
1.2 takayama 139: static void hsExec_beta(void);
1.10 takayama 140: static int genBeta(int Kap[]);
141: static int checkBeta1();
1.2 takayama 142: static int psublen(int Kap[],int Mu[]);
1.10 takayama 143: static int genJack(int M,int N);
144: static int checkJack1(int M,int N);
145: static int checkJack2(int M,int N);
146: static int mtest1b();
1.2 takayama 147:
148: static int imypower(int x,int n);
1.10 takayama 149: static int usage();
150: static int setParamDefault();
151: static int next(struct SFILE *fp,char *s,char *msg);
1.2 takayama 152: static int gopen_file(void);
153: static int setParam(char *fname);
1.8 takayama 154: static int showParam(struct SFILE *fp,int fd);
1.2 takayama 155: static double iv_factor(void);
156: static double gammam(double a,int n);
157: static double mypower(double a,int n);
158:
159: double mh_t(double A[A_LEN],double B[B_LEN],int N,int M);
160: double mh_t2(int J);
1.3 takayama 161: struct MH_RESULT *jk_main(int argc,char *argv[]);
162: int jk_freeWorkArea();
163: int jk_initializeWorkArea();
164:
165: int jk_freeWorkArea() {
166: /* bug, p in the cloneP will not be deallocated.
1.12 takayama 167: Nk in genDarray2 will not be deallocated.
168: */
1.3 takayama 169: int i;
1.6 takayama 170: JK_deallocate=1;
1.3 takayama 171: if (Darray) {myfree(Darray); Darray=NULL;}
172: if (Parray) {myfree(Parray); Parray=NULL;}
173: if (ParraySize) {myfree(ParraySize); ParraySize=NULL;}
174: if (M_beta_0) {myfree(M_beta_0); M_beta_0=NULL;}
175: if (M_beta_1) {myfree(M_beta_1); M_beta_1=NULL;}
176: if (M_jack) {
1.12 takayama 177: for (i=0; M_jack[i] != NULL; i++) {
178: if (Debug) printf("Free M_jack[%d]\n",i);
179: myfree(M_jack[i]); M_jack[i] = NULL;
180: }
181: myfree(M_jack); M_jack=NULL;
1.3 takayama 182: }
183: if (M_qk) {myfree(M_qk); M_qk=NULL;}
184: if (P_pki) {myfree(P_pki); P_pki=NULL;}
1.6 takayama 185: JK_deallocate=0;
1.13 ! takayama 186: return(0);
1.3 takayama 187: }
188: int jk_initializeWorkArea() {
189: int i,j;
1.6 takayama 190: JK_deallocate=1;
191: xval(0,0);
192: JK_deallocate=0;
1.3 takayama 193: Darray=NULL;
194: Parray=NULL;
195: ParraySize=NULL;
196: M_beta_0=NULL;
197: M_beta_1=NULL;
198: M_jack=NULL;
199: M_qk=NULL;
200: for (i=0; i<M_nmx; i++) M_kap[i]=HS_mu[i]=0;
201: for (i=0; i<M_nmx; i++) M_x[i]=0;
202: for (i=0; i<M_nmx; i++) for (j=0; j<M_m_MAX; j++) Xarray[i][j]=0;
1.5 takayama 203: for (i=0; i<M_nmx; i++) M_beta_kap[i]=0;
1.3 takayama 204: M_m=M_m_MAX-2;
205: Alpha = 2;
206: HS_n=M_nmx;
207: P_pki=NULL;
208: P_pmn=0;
209: DR_parray=0;
210: M_beta_pt=0;
211: M_df=1;
212: M_2n=0;
213: M_rel_error=0.0;
1.4 takayama 214: Sample = Sample_default;
1.5 takayama 215: Xng=0.0;
216: M_n=0;
1.13 ! takayama 217: return(0);
1.3 takayama 218: }
1.2 takayama 219:
220: static void *mymalloc(int size) {
1.1 takayama 221: void *p;
222: if (Debug) printf("mymalloc(%d)\n",size);
1.8 takayama 223: p = (void *)mh_malloc(size);
1.1 takayama 224: if (p == NULL) {
225: fprintf(stderr,"No more memory.\n");
1.3 takayama 226: mh_exit(-1);
1.1 takayama 227: }
228: return(p);
229: }
1.10 takayama 230: static int myfree(void *p) {
1.8 takayama 231: if (Debug) printf("myFree at %p\n",p);
1.13 ! takayama 232: return(mh_free(p));
1.7 takayama 233: }
1.13 ! takayama 234: static int myerror(char *s) { fprintf(stderr,"%s: type in control-C\n",s); getchar(); getchar(); return(0);}
1.1 takayama 235:
1.2 takayama 236: static double jack1(int K) {
1.1 takayama 237: double F;
238: extern int Alpha;
239: int I,J,L,II,JJ,N;
240: N = 1;
241: if (K == 0) return((double)1);
242: F = xval(1,K);
243: for (J=0; J<K; J++) {
244: II = 1; JJ = J+1;
245: F *= (N-(II-1)+Alpha*(JJ-1));
246: }
247: return(F);
248: }
1.2 takayama 249: static double jack1diff(int K) {
1.1 takayama 250: double F;
251: extern int Alpha;
252: int I,J,S,L,II,JJ,N;
253: N = 1;
254: if (K == 0) return((double)1);
255: F = K*xval(1,K-1);
256: for (J=0; J<K; J++) {
257: II = 1; JJ = J+1;
258: F *= (N-(II-1)+Alpha*(JJ-1));
259: }
260: return(F);
261: }
262:
1.2 takayama 263: static double xval(int ii,int p) { /* x_i^p */
1.1 takayama 264: extern double M_x[];
265: double F;
266: int i,j;
1.10 takayama 267: static int init=0;
1.6 takayama 268: if (JK_deallocate) { init=0; return(0.0);}
1.1 takayama 269: if (!init) {
270: for (i=1; i<=M_n; i++) {
271: for (j=0; j<M_m_MAX; j++) {
1.12 takayama 272: if (j != 0) {
273: Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1];
274: }else{
275: Xarray[i-1][0] = 1;
276: }
1.1 takayama 277: }
278: }
279: init = 1;
280: }
281: if (ii < 1) myerror("xval, index out of bound.");
282: if (p > M_m_MAX-2) myerror("xval, p is too large.");
283: if (p < 0) {
284: myerror("xval, p is negative.");
285: printf("ii=%d, p=%d\n",ii,p);
1.3 takayama 286: mh_exit(-1);
1.1 takayama 287: }
288: return(Xarray[ii-1][p]);
289: }
290:
1.2 takayama 291: static int mysum(int L[]) {
1.1 takayama 292: int S,I,N;
293: N=M_n;
294: S=0;
295: for (I=0; I<N; I++) S += L[I];
296: return(S);
297: }
298:
299: /*
1.12 takayama 300: (3,2,2,0,0) --> 3
1.1 takayama 301: */
1.2 takayama 302: static int plength(int P[]) {
1.1 takayama 303: int I;
304: for (I=0; I<M_n; I++) {
305: if (P[I] == 0) return(I);
306: }
307: return(M_n);
308: }
309: /* plength for transpose */
1.2 takayama 310: static int plength_t(int P[]) {
1.1 takayama 311: int I;
312: for (I=0; I<M_m; I++) {
313: if (P[I] == 0) return(I);
314: }
315: return(M_m);
316: }
317:
318: /*
1.12 takayama 319: ptrans(P) returns Pt
1.1 takayama 320: */
1.3 takayama 321: static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] */
1.1 takayama 322: extern int M_m;
323: int i,j,len;
1.3 takayama 324: int p[M_nmx];
1.1 takayama 325: for (i=0; i<M_n; i++) p[i] = P[i];
326: for (i=0; i<M_m+1; i++) Pt[i] = 0;
327: for (i=0; i<M_m; i++) {
328: len=plength(p); Pt[i] = len;
329: if (len == 0) return;
330: for (j=0; j<len; j++) p[j] -= 1;
331: }
332: }
333:
1.10 takayama 334: static int test_ptrans() {
1.1 takayama 335: extern int M_m;
336: int p[M_n0]={5,3,2};
337: int pt[10];
338: int i;
339: M_m = 10;
340: ptrans(p,pt);
341: if (Debug) {for (i=0; i<10; i++) printf("%d,",pt[i]); printf("\n");}
1.13 ! takayama 342: return(0);
1.1 takayama 343: }
344:
345:
346: /*
347: upper hook length
348: h_kappa^*(K)
349: */
1.2 takayama 350: static int huk(int K[],int I,int J) {
1.1 takayama 351: extern int Alpha;
352: int Kp[M_m_MAX];
353: int A,H;
354: A=Alpha;
355: /*printf("h^k(%a,%a,%a,%a)\n",K,I,J,A);*/
356: ptrans(K,Kp);
357: H=Kp[J-1]-I+A*(K[I-1]-J+1);
358: return(H);
359: }
360:
361: /*
362: lower hook length
363: h^kappa_*(K)
364: */
1.2 takayama 365: static int hdk(int K[],int I,int J) {
1.1 takayama 366: extern int Alpha;
367: int Kp[M_m_MAX];
368: int A,H;
369: A = Alpha;
370: /*printf("h_k(%a,%a,%a,%a)\n",K,I,J,A);*/
371: ptrans(K,Kp);
372: H=Kp[J-1]-I+1+A*(K[I-1]-J);
373: return(H);
374: }
375: /*
376: j_kappa. cf. Stanley.
377: */
1.2 takayama 378: static double jjk(int K[]) {
1.1 takayama 379: extern int Alpha;
380: int A,L,I,J;
381: double V;
382: A=Alpha;
383: V=1;
384: L=plength(K);
385: for (I=0; I<L; I++) {
386: for (J=0; J<K[I]; J++) {
387: V *= huk(K,I+1,J+1)*hdk(K,I+1,J+1);
388: }
389: }
390: if (Debug) {printp(K); printf("<--K, jjk=%lg\n",V);}
391: return(V);
392: }
393: /*
394: (a)_kappa^\alpha, Pochhammer symbol
395: Note that ppoch(a,[n]) = (a)_n, Alpha=2
396: */
1.2 takayama 397: static double ppoch(double A,int K[]) {
1.1 takayama 398: extern int Alpha;
399: double V;
400: int L,I,J,II,JJ;
401: V = 1;
402: L=plength(K);
403: for (I=0; I<L; I++) {
404: for (J=0; J<K[I]; J++) {
405: II = I+1; JJ = J+1;
406: V *= (A-((double)(II-1))/((double)Alpha)+JJ-1);
407: }
408: }
409: return(V);
410: }
1.2 takayama 411: static double ppoch2(double A,double B,int K[]) {
1.1 takayama 412: extern int Alpha;
413: double V;
414: int L,I,J,II,JJ;
415: V = 1;
416: L=plength(K);
417: for (I=0; I<L; I++) {
418: for (J=0; J<K[I]; J++) {
419: II = I+1; JJ = J+1;
420: V *= (A-((double)(II-1))/((double)Alpha)+JJ-1);
421: V /= (B-((double)(II-1))/((double)Alpha)+JJ-1);
422: }
423: }
424: return(V);
425: }
1.2 takayama 426: static double mypower(double x,int n) {
1.1 takayama 427: int i;
428: double v;
429: if (n < 0) return(1/mypower(x,-n));
430: v = 1;
431: for (i=0; i<n; i++) v *= x;
432: return(v);
433: }
434: /* Q_kappa
1.12 takayama 435: */
1.2 takayama 436: static double qk(int K[],double A[A_LEN],double B[B_LEN]) {
1.1 takayama 437: extern int Alpha;
438: int P,Q,I;
439: double V;
440: P = A_LEN;
441: Q = B_LEN;
442: V = mypower((double) Alpha,mysum(K))/jjk(K);
1.2 takayama 443: /* to reduce numerical errors, temporary. */
1.1 takayama 444: if (P == Q) {
445: for (I=0; I<P; I++) V = V*ppoch2(A[I],B[I],K);
446: }
447: return(V);
448:
449: for (I=0; I<P; I++) V = V*ppoch(A[I],K);
450: for (I=0; I<Q; I++) V = V/ppoch(B[I],K);
451: return(V);
452: }
453:
454: /*
1.12 takayama 455: B^nu_{kappa,mu}(i,j)
456: bb(N,K,M,I,J)
1.1 takayama 457: */
1.2 takayama 458: static int bb(int N[],int K[],int M[],int I,int J) {
1.1 takayama 459: int Kp[M_m_MAX]; int Mp[M_m_MAX];
460: ptrans(K,Kp);
461: ptrans(M,Mp);
462:
463: /*
1.12 takayama 464: printp(K); printf("K<--, "); printp2(Kp); printf("<--Kp\n");
465: printp(M); printf("M<--, "); printp2(Mp); printf("<--Mp\n");
1.1 takayama 466: */
467:
468: if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J));
469: if (Kp[J-1] == Mp[J-1]) return(huk(N,I,J));
470: else return(hdk(N,I,J));
471: }
472: /*
473: beta_{kappa,mu}
474: beta(K,M)
475: */
1.2 takayama 476: static double beta(int K[],int M[]) {
1.1 takayama 477: double V;
478: int L,I,J,II,JJ;
479: V = 1;
480:
481: L=plength(K);
482: for (I=0; I<L; I++) {
483: for (J=0; J<K[I]; J++) {
484: II = I+1; JJ = J+1;
485: V *= (double)bb(K,K,M,II,JJ);
486: /* printf("[%d,%d,%lf]\n",I,J,V); */
487: }
488: }
489:
490: L=plength(M);
491: for (I=0; I<L; I++) {
492: for (J=0; J<M[I]; J++) {
493: II = I+1; JJ = J+1;
494: V /= (double)bb(M,K,M,II,JJ);
495: /* printf("[%d,%d,%lf]\n",I,J,V);*/
496: }
497: }
498:
499: return(V);
500: }
1.10 takayama 501: static int printp(int kappa[]) {
1.1 takayama 502: int i;
503: printf("(");
504: for (i=0; i<M_n; i++) {
505: if (i <M_n-1) printf("%d,",kappa[i]);
506: else printf("%d)",kappa[i]);
507: }
1.13 ! takayama 508: return(0);
1.1 takayama 509: }
1.10 takayama 510: static int printp2(int kappa[]) {
1.1 takayama 511: int i,ell;
512: printf("(");
513: ell = plength_t(kappa);
514: for (i=0; i<ell+1; i++) {
515: if (i <ell+1-1) printf("%d,",kappa[i]);
516: else printf("%d)",kappa[i]);
517: }
1.13 ! takayama 518: return(0);
1.1 takayama 519: }
520:
1.10 takayama 521: static int test_beta() {
1.1 takayama 522: int kappa[M_n0]={2,1,0};
523: int mu1[M_n0]={1,0,0};
524: int mu2[M_n0]={1,1,0};
525: int mu3[M_n0]={2,0,0};
526: printp(kappa); printf(","); printp(mu3); printf(": beta = %lf\n",beta(kappa,mu3));
527: printp(kappa); printf(","); printp(mu1); printf(": beta = %lf\n",beta(kappa,mu1));
528: printp(kappa); printf(","); printp(mu2); printf(": beta = %lf\n",beta(kappa,mu2));
529: }
530:
531: /* main() { test_beta(); } */
532:
533:
534: /*
1.12 takayama 535: cf. w1m.rr
1.1 takayama 536: matrix hypergeometric by jack
537: N variables, up to degree M.
538: */
539: /* todo
1.12 takayama 540: def mhgj(A,B,N,M) {
541: F = 0;
542: P = partition_a(N,M);
543: for (I=0; I<length(P); I++) {
544: K = P[I];
545: F += qk(K,A,B)*jack(K,N);
546: }
547: return(F);
548: }
1.1 takayama 549: */
550:
551:
552: /* The quotient of (3.10) of Koev-Edelman K=kappa, M=mu, SK=k */
1.2 takayama 553: static double q3_10(int K[],int M[],int SK) {
1.1 takayama 554: extern int Alpha;
555: int Mp[M_m_MAX];
1.3 takayama 556: int ML[M_nmx];
557: int N[M_nmx];
1.1 takayama 558: int i,R;
559: double T,Q,V,Ur,Vr,Wr;
560: ptrans(M,Mp);
561: for (i=0; i<M_n; i++) {ML[i] = M[i]; N[i] = M[i];}
562: N[SK-1] = N[SK-1]-1;
563:
564: T = SK-Alpha*M[SK-1];
565: Q = T+1;
566: V = Alpha;
567: for (R=1; R<=SK; R++) {
568: Ur = Q-R+Alpha*K[R-1];
569: V *= Ur/(Ur+Alpha-1);
570: }
571: for (R=1; R<=SK-1; R++) {
572: Vr = T-R+Alpha*M[R-1];
573: V *= (Vr+Alpha)/Vr;
574: }
575: for (R=1; R<=M[SK-1]-1; R++) {
576: Wr = Mp[R-1]-T-Alpha*R;
577: V *= (Wr+Alpha)/Wr;
578: }
579: return(V);
580: }
581:
1.2 takayama 582: static double q3_5(double A[],double B[],int K[],int I) {
1.1 takayama 583: extern int Alpha;
584: int Kp[M_m_MAX];
585: double C,D,V,Ej,Fj,Gj,Hj,Lj;
586: int J,P,Q;
587: ptrans(K,Kp);
588: P=A_LEN;; Q = B_LEN;
589: C = -((double)(I-1))/Alpha+K[I-1]-1;
590: D = K[I-1]*Alpha-I;
591:
592: V=1;
593:
594: for (J=1; J<=P; J++) {
1.12 takayama 595: V *= (A[J-1]+C);
1.1 takayama 596: }
597: for (J=1; J<=Q; J++) {
1.12 takayama 598: V /= (B[J-1]+C);
1.1 takayama 599: }
600:
601: for (J=1; J<=K[I-1]-1; J++) {
1.12 takayama 602: Ej = D-J*Alpha+Kp[J-1];
603: Gj = Ej+1;
604: V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha));
1.1 takayama 605: }
606: for (J=1; J<=I-1; J++) {
1.12 takayama 607: Fj=K[J-1]*Alpha-J-D;
608: Hj=Fj+Alpha;
609: Lj=Hj*Fj;
610: V *= (Lj-Fj)/(Lj+Hj);
1.1 takayama 611: }
612: return(V);
613: }
614:
1.2 takayama 615: static void mtest4() {
1.1 takayama 616: double A[A_LEN] = {1.5};
617: double B[B_LEN]={6.5};
618: int K[M_n0] = {3,2,0};
619: int I=2;
620: int Ki[M_n0]={3,1,0};
621: double V1,V2;
622: V1=q3_5(A,B,K,I);
623: V2=qk(K,A,B)/qk(Ki,A,B);
624: printf("%lf== %lf?\n",V1,V2);
625: }
1.2 takayama 626: static void mtest4b() {
1.1 takayama 627: int K[M_n0]={3,2,0};
628: int M[M_n0]={2,1,0};
629: int N[M_n0]={2,0};
630: int SK=2;
631: double V1,V2;
632: V1=q3_10(K,M,SK);
633: V2=beta(K,N)/beta(K,M);
634: printf("%lf== %lf?\n",V1,V2);
635: }
636:
637: /* main() { mtest4(); mtest4b(); } */
638:
639: /* nk in (4.1),
1.12 takayama 640: */
1.2 takayama 641: static int nk(int KK[]) {
1.1 takayama 642: extern int *Darray;
643: int N,I,Ki;
1.3 takayama 644: int Kpp[M_nmx];
1.1 takayama 645: int i;
646: N = plength(KK);
647: if (N == 0) return(0);
648: if (N == 1) return(KK[0]);
649: for (i=0; i<M_n; i++) Kpp[i] = 0;
650: for (I=0; I<N-1; I++) Kpp[I] = KK[I];
651: Ki = KK[N-1];
652: /* K = (Kpp,Ki) */
653: return(Darray[nk(Kpp)]+Ki-1);
654: }
1.2 takayama 655: static int plength2(int P1[],int P2[]) {
1.1 takayama 656: int S1,S2;
657: S1 = plength(P1); S2 = plength(P2);
658: if (S1 > S2) return(1);
659: else if (S1 == S2) {
660: S1=mysum(P1); S2=mysum(P2);
661: if(S1 > S2) return(1);
662: else if (S1 == S2) return(0);
663: else return(-1);
664: }
665: else return(-1);
666: }
1.2 takayama 667: static int myeq(int P1[],int P2[]) {
1.1 takayama 668: int I,L1;
669: if ((L1=plength(P1)) != plength(P2)) return(0);
670: for (I=0; I<L1; I++) {
671: if (P1[I] != P2[I]) return(0);
672: }
673: return(1);
674: }
675: /*
676: M is a degree, N is a number of variables
677: genDarray(3,3);
678: N(0)=0;
679: N(1)=1;
680: N(2)=2;
681: N(3)=3;
682: N(1,1)=4; D[1] = 4
683: N(2,1)=5; D[2] = 5;
684: N(1,1,1)=6; D[4] = 6;
685: still buggy.
686: */
687:
1.2 takayama 688: static int pListPartition(int M,int N) {
1.1 takayama 689: extern int M_m;
690: extern int M_kap[];
691: int I;
692: /* initialize */
693: if (M_n != N) {
1.3 takayama 694: fprintf(stderr,"M_n != N\n"); mh_exit(-1);
1.1 takayama 695: }
696: M_m = M;
697: /* M_plist = []; */
698: /* end of initialize */
699: (*M_pExec)(); /* exec for 0 */
700: for (I=1; I<=M_n; I++) {
701: pListPartition2(M_m,1,I,M_m);
702: }
703: /* M_plist = reverse(M_plist); */
704: return(1);
705: }
706:
707: /*
708: Enumerate all such that
709: Less >= M_kap[From], ..., M_kap[To], |(M_kap[From],...,M_kap[To])|<=M,
710: */
1.2 takayama 711: static int pListPartition2(int Less,int From,int To, int M) {
1.1 takayama 712: int I,R;
1.11 takayama 713: mh_check_intr(100);
1.1 takayama 714: if (To < From) {
1.12 takayama 715: (*M_pExec)(); return(0);
1.1 takayama 716: }
717: for (I=1; (I<=Less) && (I<=M) ; I++) {
718: M_kap[From-1] = I;
719: R = pListPartition2(I,From+1,To,M-I);
720: }
721: return(1);
722: }
723:
724: /*
725: partition に対してやる仕事をこちらへ書く.
726: */
1.2 takayama 727: static void pExec_0() {
1.1 takayama 728: if (Debug) {
1.12 takayama 729: printf("M_kap=");
730: printp(M_kap);
731: printf("\n");
1.1 takayama 732: }
733: }
734:
735: /* Test.
1.12 takayama 736: Compare pListPartition(4,3); genDarray(4,3);
737: Compare pListPartition(5,3); genDarray(5,3);
1.1 takayama 738:
739: */
740:
741: /*
1.12 takayama 742: main() {
1.1 takayama 743: M_pExec = pExec_0;
744: pListPartition(5,3);
1.12 takayama 745: }
1.1 takayama 746: */
747:
748:
749: /*
750: List all horizontal strips.
751: Kap[0] is not a dummy in C code. !(Start from Kap[1].)
752: */
1.2 takayama 753: static int pListHS(int Kap[],int N) {
1.1 takayama 754: extern int HS_n;
755: extern int HS_mu[];
756: int i;
757: HS_n = N;
758: /* Clear HS_mu. Do not forget when N < M_n */
759: for (i=0; i<M_n; i++) HS_mu[i] = 0;
1.13 ! takayama 760: return(pListHS2(1,N,Kap));
1.1 takayama 761: }
762:
1.2 takayama 763: static int pListHS2(int From,int To,int Kap[]) {
1.1 takayama 764: int More,I;
765: if (To <From) {(*HS_hsExec)(); return(0);}
766: if (From == HS_n) More=0; else More=Kap[From];
767: for (I=Kap[From-1]; I>= More; I--) {
768: HS_mu[From-1] = I;
769: pListHS2(From+1,To,Kap);
770: }
771: return(1);
772: }
773:
1.2 takayama 774: static void hsExec_0() {
1.1 takayama 775: int i;
776: if(Debug) {printf("hsExec: "); printp(HS_mu); printf("\n");}
777: }
778:
779: /*
780: pListHS([0,4,2,1],3);
781: */
782: /*
1.12 takayama 783: main() {
1.1 takayama 784: int Kap[3]={4,2,1};
785: HS_hsExec = hsExec_0;
786: pListHS(Kap,3);
1.12 takayama 787: }
1.1 takayama 788: */
789:
790: /* The number of partitions <= M, with N parts.
1.12 takayama 791: (0,0,...,0) is excluded.
1.1 takayama 792: */
793: #define aP_pki(i,j) P_pki[(i)*(M+1)+(j)]
1.2 takayama 794: static int pmn(int M,int N) {
1.1 takayama 795: int Min_m_n,I,K,S,T,i,j;
796: extern int P_pmn;
797: extern int *P_pki;
798: Min_m_n = (M>N?N:M);
799: /* P_pki=newmat(Min_m_n+1,M+1); */
800: P_pki = (int *) mymalloc(sizeof(int)*(Min_m_n+1)*(M+1));
801: for (i=0; i<Min_m_n+1; i++) for (j=0; j<M+1; j++) aP_pki(i,j) = 0;
802: for (I=1; I<=M; I++) aP_pki(1,I) = 1;
803: for (K=1; K<=Min_m_n; K++) aP_pki(K,0) = 0;
804: S = M;
805: for (K=2; K<=Min_m_n; K++) {
806: for (I=1; I<=M; I++) {
807: if (I-K < 0) T=0; else T=aP_pki(K,I-K);
808: aP_pki(K,I) = aP_pki(K-1,I-1)+T;
809: S += aP_pki(K,I);
810: }
811: }
812: P_pmn=S;
813: if (Debug) {
814: printf("P_pmn=%d\n",P_pmn);
815: for (i=0; i<=Min_m_n; i++) {
816: for (j=0; j<=M; j++) printf("%d,",aP_pki(i,j));
817: printf("\n");
818: }
819: }
820: myfree(P_pki); P_pki=NULL;
821: return(S);
822: }
823:
824: /*
1.12 takayama 825: main() {pmn(4,3); printf("P_pmn=%d\n",P_pmn);}
1.1 takayama 826: */
827:
1.2 takayama 828: static int *cloneP(int a[]) {
1.1 takayama 829: int *p;
830: int i;
831: p = (int *) mymalloc(sizeof(int)*M_n);
832: for (i=0; i<M_n; i++) p[i] = a[i];
833: return(p);
834: }
1.10 takayama 835: static int copyP(int p[],int a[]) {
1.1 takayama 836: int i;
837: for (i=0; i<M_n; i++) p[i] = a[i];
1.13 ! takayama 838: return(0);
1.1 takayama 839: }
840:
1.2 takayama 841: static void pExec_darray(void) {
1.1 takayama 842: extern int DR_parray;
843: extern int M_kap[];
844: extern int **Parray;
845: extern int *ParraySize;
846: int *K;
847: pExec_0();
848: K = cloneP(M_kap);
849: Parray[DR_parray] = K;
850: ParraySize[DR_parray] = mysum(K);
851: DR_parray++;
852: }
1.10 takayama 853: static int genDarray2(int M,int N) {
1.1 takayama 854: extern int *Darray;
855: extern int **Parray;
856: extern int DR_parray;
857: extern int M_m;
858: int Pmn,I,J,Ksize,i;
859: int **L;
860: int *Nk;
861: int *K;
1.3 takayama 862: int Kone[M_nmx];
1.1 takayama 863:
864: M_m = M;
865: Pmn = pmn(M,N)+1;
866: if (Debug) printf("Degree M = %d, N of vars N = %d, Pmn+1=%d\n",M,N,Pmn);
867: Darray=(int *) mymalloc(sizeof(int)*Pmn);
868: for (i=0; i<Pmn; i++) Darray[i] = 0;
869: Parray=(int **) mymalloc(sizeof(int *)*Pmn);
870: for (i=0; i<Pmn; i++) Parray[i] = NULL;
871: ParraySize=(int *) mymalloc(sizeof(int *)*Pmn);
872: for (i=0; i<Pmn; i++) ParraySize[i] = 0;
873: DR_parray=0;
874: M_pExec = pExec_darray;
875: pListPartition(M,N); /* pExec_darray() is executed for all partitions */
876: L = Parray;
877:
878: Nk = (int *) mymalloc(sizeof(int)*(Pmn+1));
879: for (I=0; I<Pmn; I++) Nk[I] = I;
880: for (I=0; I<Pmn; I++) {
1.11 takayama 881: mh_check_intr(100);
1.12 takayama 882: K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */
883: Ksize = plength(K);
884: if (Ksize >= M_n) {
885: if (Debug) {fprintf(stderr,"Ksize >= M_n\n");}
886: continue;
887: }
888: for (i=0; i<M_n; i++) Kone[i] = 0;
889: for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1;
890: for (J=0; J<Pmn; J++) {
891: if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */
892: }
1.1 takayama 893: }
894: if (Debug) {
1.12 takayama 895: printf("Darray=\n");
896: for (i=0; i<Pmn; i++) printf("%d\n",Darray[i]);
897: printf("-----------\n");
1.1 takayama 898: }
1.13 ! takayama 899: return(0);
1.1 takayama 900: }
901:
902:
903: /* main() { genDarray2(4,3);} */
904:
905: /* M_beta_0[*] value of beta_{kappa,mu}, M_beta_1[*] N_mu */
1.10 takayama 906: static int isHStrip(int Kap[],int Nu[]) {
1.1 takayama 907: int N1,N2,I,P;
908: N1 = plength(Kap); N2 = plength(Nu);
909: if (N2 > N1) return(0);
910: for (I=0; I<N2; I++) {
911: if (I >= N1-1) P = 0; else P=Kap[I+1];
912: if (Kap[I] < Nu[I]) return(0);
913: if (Nu[I] < P) return(0);
914: }
915: return(1);
916: }
917:
1.2 takayama 918: static void hsExec_beta(void) {
1.1 takayama 919: int *Mu;
920: int N,Nmu,Nnu,Done,J,K,OK,I,RR;
921: int Kapt[M_m_MAX];
922: int Nut[M_m_MAX];
1.3 takayama 923: int Nu[M_nmx];
1.1 takayama 924: int rrMax;
925: hsExec_0();
926: /* printf("M_beta_pt=%a\n",M_beta_pt); */
927: /* Mu = cdr(vtol(HS_mu)); */
928: Mu = HS_mu; /* buggy? need cloneP */
929: if (M_beta_pt == 0) {
930: M_beta_0[0] = 1; M_beta_1[0] = nk(Mu);
931: M_beta_pt++; return;
932: }
933:
934: N = HS_n;
935: Nmu = nk(Mu);
936: M_beta_1[M_beta_pt] = Nmu;
937: ptrans(M_beta_kap,Kapt);
938: /* Mu, Nu is exchanged in this code. cf. the K-E paper */
939: copyP(Nu,Mu); /* buggy need clone? */
940: for (I=0; I<N; I++) {
941: Nu[I]++;
942: if (!isHStrip(M_beta_kap,Nu)) {Nu[I]--; continue;}
943: Nnu = nk(Nu);
944: ptrans(Nu,Nut);
945: Done=0;
946: for (J=M_beta_pt-1; J>=0; J--) {
947: if (M_beta_1[J] == Nnu) {
1.12 takayama 948: K=I+1;
949: if (Debug) {
950: printf("Found at J=%d, K=%d, q3_10(Kap,Nu,K)=%lf,Nu,Mu= \n",
951: J,K,q3_10(M_beta_kap,Nu,K));
952: printp(Nu); printf("\n");
953: printp(Mu); printf("\n");
954: }
955: /* Check other conditions. See Numata's mail on Dec 24, 2011. */
956: rrMax = Nu[I]-1;
957: if ((plength_t(Kapt) < rrMax) || (plength_t(Nut) < rrMax)) {
958: if (Debug) printf(" is not taken (length). \n");
959: break;
960: }
961: OK=1;
962: for (RR=0; RR<rrMax; RR++) {
963: if (Kapt[RR] != Nut[RR]) { OK=0; break;}
964: }
965: if (!OK) { if (Debug) printf(" is not taken.\n"); break; }
966: /* check done. */
967: M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K);
968: Done = 1; break;
1.1 takayama 969: }
970: }
971: if (Done) break; else Nu[I]--;
972: }
973: if (!Done) {
974: if (Debug) printf("BUG: not found M_beta_pt=%d.\n",M_beta_pt);
1.13 ! takayama 975: /* M_beta_0[M_beta_pt] = NAN; error("Not found."); */
1.1 takayama 976: M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu);
977: }
978: /* Fix the bug of mh.rr */
979: M_beta_pt++;
980: }
1.10 takayama 981: static int genBeta(int Kap[]) {
1.1 takayama 982: extern double *M_beta_0;
983: extern int *M_beta_1;
984: extern int M_beta_pt;
985: extern int M_beta_kap[];
986: extern int P_pmn;
987: int I,J,N;
988: if (Debug) {printp(Kap); printf("<-Kappa, P_pmn=%d\n",P_pmn);}
989: /* M_beta = newmat(2,P_pmn+1); */
990: M_beta_0 = (double *)mymalloc(sizeof(double)*(P_pmn+1));
1.8 takayama 991: M_beta_1 = (int *)mymalloc(sizeof(int)*(P_pmn+1));
1.1 takayama 992: M_beta_pt = 0;
993: for (I=0; I<=P_pmn; I++) {M_beta_0[I] = NAN; M_beta_1[I] = -1;}
994: N = plength(Kap);
995: HS_hsExec = hsExec_beta;
996: copyP(M_beta_kap,Kap);
997: pListHS(Kap,N);
998: }
999: /*
1000: genDarray2(4,3);
1001: genBeta([2,2,0]);
1002: genBeta([2,1,1]);
1003: */
1004:
1.10 takayama 1005: static int checkBeta1() {
1.1 takayama 1006: int Kap[3] = {2,2,0};
1007: int Kap2[3] = {2,1,0};
1008: int I;
1009: int *Mu;
1010: double Beta_km;
1011: genDarray2(4,3);
1012: genBeta(Kap);
1013: for (I=0; I<M_beta_pt; I++) {
1014: Mu = Parray[M_beta_1[I]];
1015: Beta_km = M_beta_0[I];
1.12 takayama 1016: if (Debug) {
1017: printp(Kap); printf("<--Kap, ");
1018: printp(Mu); printf("<--Mu,");
1019: printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap,Mu));
1020: }
1.1 takayama 1021: }
1022: if (Debug) printf("-------------------------------------\n");
1023: genBeta(Kap2);
1024: for (I=0; I<M_beta_pt; I++) {
1025: Mu = Parray[M_beta_1[I]];
1026: Beta_km = M_beta_0[I];
1.12 takayama 1027: if (Debug) {
1028: printp(Kap2); printf("<--Kap, ");
1029: printp(Mu); printf("<--Mu,");
1030: printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap2,Mu));
1031: }
1.1 takayama 1032: }
1.13 ! takayama 1033: return(0);
1.1 takayama 1034: }
1035:
1036: /*
1.12 takayama 1037: def checkBeta2() {
1.1 takayama 1038: genDarray2(3,3);
1039: Kap = [2,1,0];
1040: printf("Kap=%a\n",Kap);
1041: genBeta(Kap);
1042: for (I=0; I<M_beta_pt; I++) {
1.12 takayama 1043: Mu = Parray[M_beta[1][I]];
1044: Beta_km = M_beta[0][I];
1045: printf("Mu=%a,",Mu);
1046: printf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));
1047: }
1.1 takayama 1048: }
1049: */
1050:
1051: /* main() { checkBeta1(); } */
1052:
1.2 takayama 1053: static int psublen(int Kap[],int Mu[]) {
1.1 takayama 1054: int L1,L2,A,I;
1055: L1 = plength(Kap);
1056: L2 = plength(Mu);
1057: if (L2 > L1) myerror("psub, length mismatches.");
1058: A = 0;
1059: for (I=0; I<L2; I++) {
1060: if (Kap[I] < Mu[I]) myerror("psub, not Kap >= Mu");
1061: A += Kap[I]-Mu[I];
1062: }
1063: for (I=L2; I<L1; I++) A += Kap[I];
1064: return(A);
1065: }
1066:
1067:
1068: /* Table of Jack polynomials
1.12 takayama 1069: Jack[1]* one variable.
1070: Jack[2]* two variables.
1.1 takayama 1071: ...
1.12 takayama 1072: Jack[M_n]* n variables.
1073: Jack[P][J]*
1074: D^J(P variables jack of p variables). Example. J=001 d_1, 010 d_2, 100 d_3
1075: 0<=J<=2^{M_n}-1
1076: Jack[P][J][nk(Kappa)] Jack_Kappa, Kappa is a partition.
1077: 0<=nk(Kappa)<=pmn(M_m,M_n)
1.1 takayama 1078: */
1079:
1080: #define aM_jack(i,j,k) ((M_jack[i])[(j)*(Pmn+1)+(k)])
1.10 takayama 1081: static int genJack(int M,int N) {
1.1 takayama 1082: extern double **M_jack;
1083: extern int M_2n;
1084: extern int P_pmn;
1085: extern int *M_beta_1;
1086: int Pmn,I,J,K,L,Nv,H,P;
1087: int *Kap,*Mu;
1088: double Jack,Beta_km;
1089: int Nk,JJ;
1.5 takayama 1090: if (Debug) printf("genJack(%d,%d)\n",M,N);
1.3 takayama 1091: M_jack = (double **) mymalloc(sizeof(double *)*(N+2));
1.1 takayama 1092: M_2n = imypower(2,N);
1093: Pmn = pmn(M,N); /*P_pmn is initializeded.
1094: Warning. It is reset when pmn is called.*/
1095: 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 1096: M_jack[N+1] = NULL;
1.1 takayama 1097: genDarray2(M,N); /* Darray, Parray is initialized */
1098: for (I=1; I<=N; I++) aM_jack(I,0,0) = 1;
1099: if (M_df) {
1100: for (I=1; I<=N; I++) {
1101: for (J=1; J<M_2n; J++) aM_jack(I,J,0) = 0;
1102: }
1103: }
1104:
1.3 takayama 1105: /* N must satisfies N > 0 */
1.1 takayama 1106: for (K=1; K<=M; K++) {
1107: aM_jack(1,0,K) = jack1(K);
1108: if (M_df) {
1109: aM_jack(1,1,K) = jack1diff(K); /* diff(jack([K],1),x_1); */
1110: for (J=2; J<M_2n; J++) aM_jack(1,J,K) = 0;
1111: }
1112: }
1113: for (I=1; I<=N; I++) {
1114: for (K=M+1; K<Pmn+1; K++) {
1115: aM_jack(I,0,K) = NAN;
1116: if (M_df) {
1117: for (J=1; J<M_2n; J++) {
1118: if (J >= 2^I) aM_jack(I,J,K) = 0;
1.12 takayama 1119: else aM_jack(I,J,K) = NAN;
1.1 takayama 1120: }
1121: }
1122: }
1123: }
1124:
1125: /* Start to evaluate the entries of the table */
1126: for (K=1; K<=Pmn; K++) {
1127: Kap = Parray[K]; /* bug. need copy? */
1128: L = plength(Kap);
1129: for (I=1; I<=L-1; I++) {
1130: aM_jack(I,0,K) = 0;
1131: if (M_df) {
1132: for (J=1; J<M_2n; J++) aM_jack(I,J,K) = 0;
1133: }
1134: }
1135: if (Debug) {printf("Kappa="); printp(Kap);}
1136: /* Enumerate horizontal strip of Kappa */
1137: genBeta(Kap); /* M_beta_pt stores the number of hs */
1138: /* Nv is the number of variables */
1139: for (Nv = (L==1?2:L); Nv <= N; Nv++) {
1140: Jack = 0;
1141: for (H=0; H<M_beta_pt; H++) {
1142: Nk = M_beta_1[H];
1143: Mu = Parray[Nk];
1.12 takayama 1144: if (UseTable) {
1145: Beta_km = M_beta_0[H];
1146: }else{
1147: Beta_km = beta(Kap,Mu);
1148: /* do not use the M_beta table. It's buggy. UseTable is experimental.*/
1149: }
1.1 takayama 1150: if (Debug) {printf("Nv(number of variables)=%d, Beta_km=%lf, Mu=",Nv,Beta_km);
1.12 takayama 1151: printp(Mu); printf("\n");}
1.1 takayama 1152: P = psublen(Kap,Mu);
1153: Jack += aM_jack(Nv-1,0,Nk)*Beta_km*xval(Nv,P); /* util_v(x,[Nv])^P;*/
1.12 takayama 1154: if (Debug) printf("xval(%d,%d)=%lf\n",Nv,P,xval(Nv,P));
1.1 takayama 1155: }
1156: aM_jack(Nv,0,K) = Jack;
1157: if (M_df) {
1158: /* The case of M_df > 0. */
1159: for (J=1; J<M_2n; J++) {
1.12 takayama 1160: mh_check_intr(100);
1161: Jack = 0;
1162: for (H=0; H<M_beta_pt; H++) {
1163: Nk = M_beta_1[H];
1164: Mu = Parray[Nk];
1165: if (UseTable) {
1166: Beta_km = M_beta_0[H];
1167: }else{
1168: Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */
1169: }
1170: if (Debug) {printf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);
1171: printp(Mu); printf("\n"); }
1172: P = psublen(Kap,Mu);
1173: if (J & (1 << (Nv-1))) {
1174: JJ = J & ((1 << (Nv-1)) ^ 0xffff); /* NOTE!! Up to 16 bits. mh-15 */
1175: if (P != 0) {
1176: Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1);
1177: }
1178: }else{
1179: Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P);
1180: }
1181: }
1182: aM_jack(Nv,J,K) = Jack;
1183: if (Debug) printf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);
1.1 takayama 1184: } /* end of J loop */
1185: }
1186: }
1187: }
1188: }
1189:
1190:
1191: /* checkJack1(3,3)
1.12 takayama 1192: */
1.10 takayama 1193: static int checkJack1(int M,int N) {
1.1 takayama 1194: int I,K;
1195: extern int P_pmn;
1196: extern double M_x[];
1197: int Pmn; /* used in aM_jack */
1198: /* initialize x vars. */
1199: for (I=1; I<=N; I++) {
1200: M_x[I-1] = ((double)I)/10.0;
1201: }
1202: genJack(M,N);
1203: Pmn = P_pmn;
1204: for (I=1; I<=N; I++) {
1205: for (K=0; K<=P_pmn; K++) {
1206: printp(Parray[K]);
1207: printf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
1208: }
1209: }
1210: for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);
1211: printf("<--x\n");
1.13 ! takayama 1212: return(0);
1.1 takayama 1213: }
1214: /*main() { checkJack1(3,3); }*/
1215:
1216:
1.10 takayama 1217: static int checkJack2(int M,int N) {
1.1 takayama 1218: int I,K,J;
1219: extern int P_pmn;
1220: extern double M_x[];
1.10 takayama 1221: extern int M_df;
1.1 takayama 1222: int Pmn; /* used in aM_jack */
1223: M_df=1;
1224: /* initialize x vars. */
1225: for (I=1; I<=N; I++) {
1226: M_x[I-1] = ((double)I)/10.0;
1227: }
1228: genJack(M,N);
1229: Pmn = P_pmn;
1230: for (I=1; I<=N; I++) {
1231: for (K=0; K<=P_pmn; K++) {
1232: printp(Parray[K]);
1233: printf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
1234: }
1235: }
1236: for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);
1237: printf("<--x\n");
1238:
1239: for (I=1; I<=N; I++) {
1240: for (K=0; K<=P_pmn; K++) {
1241: for (J=0; J<M_2n; J++) {
1.12 takayama 1242: printp(Parray[K]);
1243: printf("<--Kap, Nv=%d,J(diff)=%d, D^J Jack=%lf\n",
1244: I,J,aM_jack(I,J,K));
1245: }
1246: }
1.1 takayama 1247: }
1.13 ! takayama 1248: return(0);
1.1 takayama 1249: }
1250:
1251: /* main() { checkJack2(3,3); } */
1252:
1253: double mh_t(double A[A_LEN],double B[B_LEN],int N,int M) {
1254: double F,F2;
1255: extern int M_df;
1256: extern int P_pmn;
1257: extern double *M_qk;
1258: extern double M_rel_error;
1259: int Pmn;
1260: int K;
1261: int *Kap;
1262: int size;
1263: F = 0; F2=0;
1264: M_df=1;
1265: genJack(M,N);
1.8 takayama 1266: M_qk = (double *)mymalloc(sizeof(double)*(P_pmn+1)); /* found a bug. */
1.1 takayama 1267: Pmn = P_pmn;
1268: size = ParraySize[P_pmn];
1269: for (K=0; K<=P_pmn; K++) {
1.11 takayama 1270: mh_check_intr(100);
1.12 takayama 1271: Kap = Parray[K];
1272: M_qk[K] = qk(Kap,A,B);
1273: F += M_qk[K]*aM_jack(N,0,K);
1274: if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K);
1275: if (Debug) printf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size);
1276: if (Debug && (ParraySize[K] == size)) printf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K));
1.1 takayama 1277: }
1278: M_rel_error = F-F2;
1279: return(F);
1280: }
1281: double mh_t2(int J) {
1282: extern double *M_qk;
1283: double F;
1284: int K;
1285: int Pmn;
1286: extern int P_pmn;
1.3 takayama 1287: if (M_qk == NULL) {myerror("Call mh_t first."); mh_exit(-1); }
1.1 takayama 1288: F = 0;
1289: Pmn = P_pmn;
1290: for (K=0; K<P_pmn; K++) {
1.12 takayama 1291: F += M_qk[K]*aM_jack(M_n,J,K);
1.1 takayama 1292: }
1293: return(F);
1294: }
1295:
1.10 takayama 1296: static int mtest1b() {
1.1 takayama 1297: double A[1] = {1.5};
1298: double B[1] = {1.5+5};
1299: int I,N,M,J;
1300: double F;
1301: N=3; M=6;
1302: for (I=1; I<=N; I++) {
1303: M_x[I-1] = ((double)I)/10.0;
1304: }
1305: mh_t(A,B,N,M);
1306: for (J=0; J<M_2n; J++) {
1.12 takayama 1307: F=mh_t2(J);
1308: printf("J=%d, D^J mh_t=%lf\n",J,F);
1.1 takayama 1309: }
1310: }
1311:
1312: /* main() { mtest1b(); }*/
1313:
1314:
1315:
1316:
1317: #define TEST 1
1318: #ifndef TEST
1319:
1320: #endif
1321:
1322: /****** from mh-n.c *****/
1323:
1324: #define SMAX 4096
1.3 takayama 1325: #define inci(i) { i++; if (i >= argc) { fprintf(stderr,"Option argument is not given.\n"); return(NULL); }}
1.1 takayama 1326:
1.2 takayama 1327: static int imypower(int x,int n) {
1.1 takayama 1328: int i;
1329: int v;
1.3 takayama 1330: if (n < 0) {myerror("imypower"); mh_exit(-1);}
1.1 takayama 1331: v = 1;
1332: for (i=0; i<n; i++) v *= x;
1333: return(v);
1334: }
1335:
1.11 takayama 1336: #ifdef STANDALONE2
1.1 takayama 1337: main(int argc,char *argv[]) {
1.3 takayama 1338: mh_exit(MH_RESET_EXIT);
1.12 takayama 1339: /* jk_main(argc,argv);
1340: printf("second run.\n"); */
1.2 takayama 1341: jk_main(argc,argv);
1342: }
1.3 takayama 1343: #endif
1344:
1345: struct MH_RESULT *jk_main(int argc,char *argv[]) {
1.1 takayama 1346: double *y0;
1347: double x0,xn;
1348: double ef;
1349: int i,j,rank;
1350: double a[1]; double b[1];
1351: extern double M_x[];
1352: extern double *Beta;
1353: extern int M_2n;
1.3 takayama 1354: char swork[1024];
1355: struct MH_RESULT *ans=NULL;
1356: struct SFILE *ofp = NULL;
1357: int idata=0;
1358: JK_byFile = 1;
1359: jk_initializeWorkArea();
1.1 takayama 1360: UseTable = 1;
1.3 takayama 1361: Mapprox=6;
1.1 takayama 1362: for (i=1; i<argc; i++) {
1.12 takayama 1363: if (strcmp(argv[i],"--idata")==0) {
1364: inci(i);
1365: setParam(argv[i]); idata=1;
1366: }else if (strcmp(argv[i],"--degree")==0) {
1367: inci(i);
1368: sscanf(argv[i],"%d",&Mapprox);
1369: }else if (strcmp(argv[i],"--x0")==0) {
1370: inci(i);
1371: sscanf(argv[i],"%lg",&X0g);
1372: }else if (strcmp(argv[i],"--notable")==0) {
1373: UseTable = 0;
1374: }else if (strcmp(argv[i],"--debug")==0) {
1375: Debug = 1;
1376: }else if (strcmp(argv[i],"--help")==0) {
1377: usage(); return(0);
1378: }else if (strcmp(argv[i],"--bystring")==0) {
1379: if (idata) {fprintf(stderr,"--bystring must come before --idata option.\n"); mh_exit(-1);}
1380: JK_byFile = 0;
1381: }else {
1382: fprintf(stderr,"Unknown option %s\n",argv[i]);
1383: usage();
1384: return(NULL);
1385: }
1.1 takayama 1386: }
1.3 takayama 1387: if (!idata) setParam(NULL);
1388:
1389: /* Initialize global variables */
1390: M_n = Mg;
1391: HS_n=M_n;
1392: if (!JK_byFile) ans = (struct MH_RESULT *)mymalloc(sizeof(struct MH_RESULT));
1393: else ans = NULL;
1394: /* Output by a file=stdout */
1395: ofp = mh_fopen("stdout","w",JK_byFile);
1.1 takayama 1396:
1.3 takayama 1397: sprintf(swork,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp);
1398: sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp);
1.1 takayama 1399: if (M_n != Mg) {
1.12 takayama 1400: myerror("Mg must be equal to M_n\n"); mh_exit(-1);
1.1 takayama 1401: }
1.8 takayama 1402: if (Debug) showParam(NULL,1);
1.1 takayama 1403: for (i=0; i<M_n; i++) {
1.12 takayama 1404: M_x[i] = Beta[i]*X0g;
1.1 takayama 1405: }
1406: a[0] = ((double)Mg+1.0)/2.0;
1407: b[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */
1408: if (Debug) printf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",a[0],b[0],M_n,Mapprox);
1409: mh_t(a,b,M_n,Mapprox);
1410: if (imypower(2,M_n) != M_2n) {
1.12 takayama 1411: sprintf(swork,"M_n=%d,M_2n=%d\n",M_n,M_2n); mh_fputs(swork,ofp);
1412: myerror("2^M_n != M_2n\n"); mh_exit(-1);
1.1 takayama 1413: }
1.3 takayama 1414: sprintf(swork,"%%%%M_rel_error=%lg\n",M_rel_error); mh_fputs(swork,ofp);
1.1 takayama 1415: for (j=0; j<M_2n; j++) {
1.12 takayama 1416: Iv[j] = mh_t2(j);
1.1 takayama 1417: }
1418: Ef = iv_factor();
1.8 takayama 1419: showParam(ofp,0);
1.3 takayama 1420:
1421: /* return the result */
1422: if (!JK_byFile) {
1.12 takayama 1423: ans->x = X0g;
1424: ans->rank = imypower(2,Mg);
1425: ans->y = (double *)mymalloc(sizeof(double)*(ans->rank));
1426: for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i];
1427: ans->size=1;
1428: ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size));
1429: (ans->sfpp)[0] = ofp;
1.3 takayama 1430: }
1.7 takayama 1431: if (Debug) printf("jk_freeWorkArea() starts\n");
1.3 takayama 1432: jk_freeWorkArea();
1.7 takayama 1433: if (Debug) printf("jk_freeWorkArea() has finished.\n");
1.3 takayama 1434: return(ans);
1.1 takayama 1435: }
1436:
1.10 takayama 1437: static int usage() {
1.1 takayama 1438: fprintf(stderr,"Usages:\n");
1439: fprintf(stderr,"mh-m [--idata input_data_file --x0 x0 --degree approxm]\n");
1440: 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");
1441: fprintf(stderr,"The mh-m uses the Koev-Edelman algorithm to evalute the matrix hypergeometric function.\n");
1442: fprintf(stderr,"The degree of the approximation (Mapprox) is given by the --degree option.\n");
1443: fprintf(stderr,"Parameters are specified by the input_data_file. Otherwise, default values are used.\n");
1444: fprintf(stderr,"The format of the input_data_file. The orders of the input numbers must be kept.\n");
1445: fprintf(stderr," Mg: m(the number of variables), Beta: beta=sigma^(-1)/2 (diagonized), Ng: n,\n");
1446: fprintf(stderr," Iv: initial values at X0g*Beta (see our paper how to order them), are evaluated in this program. Give zeros. \n");
1447: fprintf(stderr," Ef: a scalar factor to the initial value. It is calculated by this program. Give the zero.\n");
1448: 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");
1449: fprintf(stderr," Dp: output data is stored in every Dp steps when output_data_file is specified, which is for w-m.\n");
1450: fprintf(stderr," The line started with %% is a comment line.\n");
1451: 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");
1452: fprintf(stderr," An example format of the input_data_file can be obtained by executing mh-m with no option.\n");
1453: fprintf(stderr,"\nExamples:\n");
1454: fprintf(stderr,"[1] ./mh-2 \n");
1455: fprintf(stderr,"[2] ./mh-2 --x0 0.3 \n");
1456: fprintf(stderr,"[3] ./mh-3 --x0 0.1 \n");
1457: fprintf(stderr,"[4] ./mh-3 --x0 0.1 --degree 15 \n");
1458: fprintf(stderr,"[5] ./mh-3 --idata test.txt --degree 15 \n");
1459: fprintf(stderr,"[6] ./mh-3 --degree 15 >test2.txt\n");
1460: fprintf(stderr," ./w-3 --idata test2.txt --gnuplotf test-g\n");
1461: fprintf(stderr," gnuplot -persist <test-g-gp.txt\n");
1.13 ! takayama 1462: return(0);
1.1 takayama 1463: }
1464:
1.10 takayama 1465: static int setParamDefault() {
1.1 takayama 1466: int rank;
1467: int i;
1.3 takayama 1468: Mg = M_n_default ;
1.1 takayama 1469: rank = imypower(2,Mg);
1470: Beta = (double *)mymalloc(sizeof(double)*Mg);
1471: for (i=0; i<Mg; i++) Beta[i] = 1.0+i;
1472: Ng = (double *)mymalloc(sizeof(double)); *Ng = 3.0;
1473: Iv = (double *)mymalloc(sizeof(double)*rank);
1474: Iv2 = (double *)mymalloc(sizeof(double)*rank);
1475: for (i=0; i<rank; i++) Iv[i] = 0;
1476: Ef = 0;
1477: Ef2 = 0.01034957388338225707;
1478: if (M_n == 2) {
1479: Iv2[0] = 1.58693;
1480: Iv2[1] = 0.811369;
1481: Iv2[2] = 0.846874;
1482: Iv2[3] = 0.413438;
1483: }
1484: X0g = (Beta[0]/Beta[Mg-1])*0.5;
1485: Hg = 0.001;
1486: Dp = 1;
1487: Xng = 10.0;
1.13 ! takayama 1488: return(0);
1.1 takayama 1489: }
1490:
1.10 takayama 1491: static int next(struct SFILE *sfp,char *s,char *msg) {
1.1 takayama 1492: s[0] = '%';
1493: while (s[0] == '%') {
1.12 takayama 1494: if (!mh_fgets(s,SMAX,sfp)) {
1495: fprintf(stderr,"Data format error at %s\n",msg);
1496: mh_exit(-1);
1497: }
1498: if (s[0] != '%') return(0);
1.1 takayama 1499: }
1.13 ! takayama 1500: return(0);
1.1 takayama 1501: }
1.10 takayama 1502: static int setParam(char *fname) {
1.1 takayama 1503: int rank;
1504: char s[SMAX];
1.3 takayama 1505: struct SFILE *fp;
1.1 takayama 1506: int i;
1507: if (fname == NULL) return(setParamDefault());
1508:
1509: Sample = 0;
1.3 takayama 1510: if ((fp=mh_fopen(fname,"r",JK_byFile)) == NULL) {
1.12 takayama 1511: if (JK_byFile) fprintf(stderr,"File %s is not found.\n",fp->s);
1512: mh_exit(-1);
1.1 takayama 1513: }
1514: next(fp,s,"Mg(m)");
1515: sscanf(s,"%d",&Mg);
1516: rank = imypower(2,Mg);
1517:
1518: Beta = (double *)mymalloc(sizeof(double)*Mg);
1519: for (i=0; i<Mg; i++) {
1520: next(fp,s,"Beta");
1.12 takayama 1521: sscanf(s,"%lf",&(Beta[i]));
1.1 takayama 1522: }
1523:
1524: Ng = (double *)mymalloc(sizeof(double));
1525: next(fp,s,"Ng(freedom parameter n)");
1526: sscanf(s,"%lf",Ng);
1527:
1528: next(fp,s,"X0g(initial point)");
1529: sscanf(s,"%lf",&X0g);
1530:
1531: Iv = (double *)mymalloc(sizeof(double)*rank);
1532: for (i=0; i<rank; i++) {
1.12 takayama 1533: next(fp,s,"Iv(initial values)");
1534: sscanf(s,"%lg",&(Iv[i]));
1.1 takayama 1535: }
1536:
1537: next(fp,s,"Ef(exponential factor)");
1538: sscanf(s,"%lg",&Ef);
1539:
1540: next(fp,s,"Hg (step size of rk)");
1541: sscanf(s,"%lf",&Hg);
1542:
1543: next(fp,s,"Dp (data sampling period)");
1544: sscanf(s,"%d",&Dp);
1545:
1546: next(fp,s,"Xng (the last point, cf. --xmax)");
1547: sscanf(s,"%lf",&Xng);
1.3 takayama 1548: mh_fclose(fp);
1.13 ! takayama 1549: return(0);
1.1 takayama 1550: }
1551:
1.10 takayama 1552: static int showParam(struct SFILE *fp,int fd) {
1.1 takayama 1553: int rank,i;
1.3 takayama 1554: char swork[1024];
1.8 takayama 1555: if (fd) {
1556: fp = mh_fopen("stdout","w",1);
1557: }
1.1 takayama 1558: rank = imypower(2,Mg);
1.3 takayama 1559: sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp);
1.1 takayama 1560: for (i=0; i<Mg; i++) {
1.12 takayama 1561: sprintf(swork,"%%Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp);
1.1 takayama 1562: }
1.3 takayama 1563: sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp);
1564: sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp);
1.1 takayama 1565: for (i=0; i<rank; i++) {
1.12 takayama 1566: sprintf(swork,"%%Iv[%d]=\n%lg\n",i,Iv[i]); mh_fputs(swork,fp);
1567: if (Sample && (M_n == 2) && (X0g == 0.3)) {
1568: sprintf(swork,"%%Iv[%d]-Iv2[%d]=%lg\n",i,i,Iv[i]-Iv2[i]); mh_fputs(swork,fp);
1569: }
1.3 takayama 1570: }
1571: sprintf(swork,"%%Ef=\n%lg\n",Ef); mh_fputs(swork,fp);
1572: sprintf(swork,"%%Hg=\n%lf\n",Hg); mh_fputs(swork,fp);
1573: sprintf(swork,"%%Dp=\n%d\n",Dp); mh_fputs(swork,fp);
1574: sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp);
1.13 ! takayama 1575: return(0);
1.1 takayama 1576: }
1577:
1.2 takayama 1578: static double gammam(double a,int n) {
1.1 takayama 1579: double v,v2;
1580: int i;
1581: v=mypower(sqrt(M_PI),(n*(n-1))/2);
1582: v2=0;
1583: for (i=1; i<=n; i++) {
1584: v2 += lgamma(a-((double)(i-1))/2.0); /* not for big n */
1585: }
1586: if (Debug) printf("gammam(%lg,%d)=%lg\n",a,n,v*exp(v2));
1587: return(v*exp(v2));
1588: }
1589:
1.2 takayama 1590: static double iv_factor(void) {
1.1 takayama 1591: double v1;
1592: double t;
1593: double b;
1594: double detSigma;
1595: double c;
1596: int i,n;
1597: n = (int) (*Ng);
1598: v1= mypower(sqrt(X0g),n*M_n);
1599: t = 0.0;
1600: for (i=0; i<M_n; i++) t += -X0g*Beta[i];
1601: v1 = v1*exp(t);
1602:
1603: b = 1; for (i=0; i<M_n; i++) b *= Beta[i];
1604: detSigma = 1.0/(b*mypower(2.0,M_n));
1605:
1606: c = gammam(((double)(M_n+1))/2.0, M_n)/
1.12 takayama 1607: ( mypower(sqrt(2), M_n*n)*mypower(sqrt(detSigma),n)*gammam(((double)(M_n+n+1))/2.0,M_n) );
1.1 takayama 1608: return( c*v1);
1609: }
1610:
1611:
1612:
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