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