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