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