Annotation of OpenXM/src/asir-contrib/testing/noro/wishartd.rr, Revision 1.2
1.2 ! noro 1: /* $OpenXM$ */
1.1 noro 2: /* A package for computing PDEs satisfied by a matrix 1F1 on diagonal regions */
3: if (version()<20151202) {error("The version of Risa/Asir must be after 20151202 to run this package.");} else{}
4: module n_wishartd$
5: localf compf1$
6: localf my_comp_f$
7: localf compc1, compc, compd, compt, addpf, subpf, addc, ftorat, ctorat, mulf$
8: localf lcmf, mulf, divf, mulc, mulcpf, diffc, muldpf, normalizef1, normalizec1, normalizepf$
9: localf wsetup, mtot, wishartpf, degf, removef, subst_f, lopitalpf, simpop, simppf$
10: localf eliminatetop, diagpf, diag0pf, adjop1, adjop, reducepf, reduceotherspf$
11: localf print_f, printc, printt, printpf, ftop, ctord, comppfrd, multpf, spolypf, nfpf$
12: localf nfpf0$
13: localf substc$
14: localf dpf,dpf2,dpf3,abs$
15: localf pftord$
16: localf lopital_others_pf$
17: localf pftop$
18: localf dnc,pftozpf,addzpf,zpftopf$
19: localf mul_lopitalpf$
20: localf indep_eqs$
21: localf mulpf$
22: localf pftoeuler$
23: localf gammam,prc,prc2$
24: localf ps,psvalue,act_op,decomp_op,create_homo,bpsvalue$
25: localf pfaffian,gen_basis$
26: localf evalpf,eval_pfaffian,genrk45$
27: localf solve_leq,ptom$
28: localf msubst$
29: localf vton,encodef1,encodef,encodec1,encodec$
30: localf vton,ftotex,ctotex,pftotex,optotex,eqtotex$
31: localf genprc,prob_by_hgm,prob_by_ps$
32: localf partition, all_partition, coef_partition, count_dup,simp_power1$
33: localf lop1,lopn,loph,mul_lopitalpf2$
34: localf muldmpf$
35: localf diagcheck,eliminatecheck$
36: localf getchi$
37: localf message$
38: localf pfaffian2, eval_pfaffian2$
39:
40: /*
41: * pfpoly = [[C,<<...>>],...]
42: * C = [[N,L],...] L = [[F,M],...]
43: * lc(F) = 1
44: */
45:
46: static Print$
47: static PF_N,PF_V,PF_DV$
48: static Tlopital,Tred,Tlineq,Tactmul$
49: static Tp,Tps,Trk$
50:
51: /* XXX execute ord([y1,y2,...,dy1,dy2,...]) */
52:
53: /* F1>F2 <=> var(F1)>var(F2) || var(F1)=var(F2) && rest(F1)>rest(F2) */
54: /* F1^M1 > F2^M2 <=> F1>F2 || F1=F2 && M1>M2 */
55:
56: def abs(A) { return A>0 ? A : -A; }
57:
58: def compf1(E1,E2)
59: {
60: F1 = E1[0]; F2 = E2[0];
61: if ( F1>F2 ) return 1;
62: if ( F1<F2 ) return -1;
63: M1 = E1[1]; M2 = E2[1];
64: if ( M1 > M2 ) return 1;
65: if ( M1 < M2 ) return -1;
66: return 0;
67: }
68:
69: def compc1(ND1,ND2)
70: {
71: if ( R = comp_f(ND1[1],ND2[1]) ) return R;
72: N1 = ND1[0]; N2 = ND2[0];
73: if ( N1 > N2 ) return 1;
74: if ( N1 < N2 ) return -1;
75: return 0;
76: }
77:
78: /* compare ND lists */
79: def compc(L1,L2)
80: {
81: for ( ; L1 != [] && L2 != []; L1 = cdr(L1), L2 = cdr(L2) ) {
82: E1 = car(L1); E2 = car(L2);
83: if ( R = compc1(E1,E2) ) return R;
84: }
85: if ( L1 != [] ) return 1;
86: if ( L2 != [] ) return -1;
87: return 0;
88: }
89:
90: def compd(D1,D2)
91: {
92: if ( D1 > D2 ) return 1;
93: if ( D1 < D2 ) return -1;
94: return 0;
95: }
96:
97: /* compare terms */
98: def compt(T1,T2)
99: {
100: if ( R = compd(T1[1],T2[1]) ) return R;
101: if ( R = compc(T1[0],T2[0]) ) return R;
102: return 0;
103: }
104:
105: def addpf(F1,F2)
106: {
107: R = [];
108: while ( F1 != [] && F2 != [] ) {
109: T1 = car(F1); T2 = car(F2);
110: C = compd(T1[1],T2[1]);
111: if ( C > 0 ) {
112: R = cons(T1,R); F1 = cdr(F1);
113: } else if ( C < 0 ) {
114: R = cons(T2,R); F2 = cdr(F2);
115: } else {
116: S = addc(T1[0],T2[0]);
117: if ( S != [] )
118: R = cons([S,T1[1]],R);
119: F1 = cdr(F1); F2 = cdr(F2);
120: }
121: }
122: R = reverse(R);
123: if ( F1 != [] ) R = append(R,F1);
124: else if ( F2 != [] ) R = append(R,F2);
125: return R;
126: }
127:
128: def subpf(F1,F2)
129: {
130: T = mulcpf([[-1,[]]],F2);
131: T = addpf(F1,T);
132: return T;
133: }
134:
135: def addc(F1,F2)
136: {
137: R = [];
138: while ( F1 != [] && F2 != [] ) {
139: T1 = car(F1); T2 = car(F2);
140: C = comp_f(T1[1],T2[1]);
141: if ( !T1[0] || !T2[0] )
142: error("afo");
143: if ( C > 0 ) {
144: R = cons(T1,R); F1 = cdr(F1);
145: } else if ( C < 0 ) {
146: R = cons(T2,R); F2 = cdr(F2);
147: } else {
148: S = T1[0]+T2[0];
149: if ( S )
150: R = cons([S,T1[1]],R);
151: F1 = cdr(F1); F2 = cdr(F2);
152: }
153: }
154: R = reverse(R);
155: if ( F1 != [] ) R = append(R,F1);
156: else if ( F2 != [] ) R = append(R,F2);
157: return R;
158: }
159:
160: def ftorat(F)
161: {
162: R = 1;
163: for ( ; F != []; F = cdr(F) ) {
164: F0 = car(F);
165: R *= F0[0]^F0[1];
166: }
167: return R;
168: }
169:
170: def ctorat(C)
171: {
172: R = 0;
173: for ( ; C != []; C = cdr(C) ) {
174: C0 = car(C);
175: R += C0[0]/ftorat(C0[1]);
176: R = red(R);
177: }
178: return R;
179: }
180:
181: def mulf(L1,L2)
182: {
183: R = [];
184: while ( L1 != [] && L2 != [] ) {
185: A1 = car(L1); A2 = car(L2);
186: if ( A1[0] > A2[0] ) {
187: R = cons(A1,R); L1 = cdr(L1);
188: } else if ( A1[0] < A2[0] ) {
189: R = cons(A2,R); L2 = cdr(L2);
190: } else {
191: R = cons([A1[0],A1[1]+A2[1]],R);
192: L1 = cdr(L1); L2 = cdr(L2);
193: }
194: }
195: R = reverse(R);
196: if ( L2 == [] ) return append(R,L1);
197: else if ( L1 == [] ) return append(R,L2);
198: else return R;
199: }
200:
201: def lcmf(L1,L2)
202: {
203: R = [];
204: while ( L1 != [] && L2 != [] ) {
205: A1 = car(L1); A2 = car(L2);
206: if ( A1[0] > A2[0] ) {
207: R = cons(A1,R); L1 = cdr(L1);
208: } else if ( A1[0] < A2[0] ) {
209: R = cons(A2,R); L2 = cdr(L2);
210: } else {
211: M = A1[1]>A2[1]?A1[1]:A2[1];
212: R = cons([A1[0],M],R);
213: L1 = cdr(L1); L2 = cdr(L2);
214: }
215: }
216: R = reverse(R);
217: if ( L2 == [] ) return append(R,L1);
218: else if ( L1 == [] ) return append(R,L2);
219: else return R;
220: }
221:
222: def mulf(L1,L2)
223: {
224: R = [];
225: while ( L1 != [] && L2 != [] ) {
226: A1 = car(L1); A2 = car(L2);
227: if ( A1[0] > A2[0] ) {
228: R = cons(A1,R); L1 = cdr(L1);
229: } else if ( A1[0] < A2[0] ) {
230: R = cons(A2,R); L2 = cdr(L2);
231: } else {
232: R = cons([A1[0],A1[1]+A2[1]],R);
233: L1 = cdr(L1); L2 = cdr(L2);
234: }
235: }
236: R = reverse(R);
237: if ( L2 == [] ) return append(R,L1);
238: else if ( L1 == [] ) return append(R,L2);
239: else return R;
240: }
241:
242: def divf(L1,L2)
243: {
244: R = [];
245: while ( L1 != [] && L2 != [] ) {
246: A1 = car(L1); A2 = car(L2);
247: if ( A1[0] > A2[0] ) {
248: R = cons(A1,R); L1 = cdr(L1);
249: } else if ( A1[0] < A2[0] ) {
250: error("divf : cannot happen");
251: } else {
252: M = A1[1]-A2[1];
253: if ( M < 0 )
254: error("divf : cannot happen");
255: if ( M > 0 )
256: R = cons([A1[0],M],R);
257: L1 = cdr(L1); L2 = cdr(L2);
258: }
259: }
260: R = reverse(R);
261: if ( L2 == [] ) return append(R,L1);
262: else if ( L1 == [] ) return append(R,L2);
263: else return R;
264: }
265:
266: def mulc(C1,C2)
267: {
268: R = [];
269: C1 = reverse(C1);
270: for ( ; C1 != []; C1 = cdr(C1) ) {
271: S = [];
272: A1 = car(C1);
273: for ( T = C2 ; T != []; T = cdr(T) ) {
274: A2 = car(T);
275: S = cons([A1[0]*A2[0],mulf(A1[1],A2[1])],S);
276: }
277: S = reverse(S);
278: R = addc(R,S);
279: }
280: return R;
281: }
282:
283: def vton(V)
284: {
285: for ( I = 1; I <= PF_N; I++ )
286: if ( V == PF_V[I] ) return I;
287: error("vton : no such variable");
288: }
289:
290: def encodef1(F)
291: {
292: A = F[0];
293: V1 = var(A);
294: N1 = vton(V1);
295: R = A-V1;
296: if ( !R )
297: return [N1,F[1]];
298: else
299: return [N1*65536+vton(var(R)),F[1]];
300: }
301:
302: def encodef(F)
303: {
304: return map(encodef1,F);
305: }
306:
307: def encodec1(C)
308: {
309: return cons(C[0],encodef(C[1]));
310: }
311:
312: def encodec(C)
313: {
314: R = map(encodec1,C);
315: }
316:
317: def mulcpf(C,F)
318: {
319: R = [];
320: for ( ; F != []; F = cdr(F) ) {
321: F0 = car(F);
322: C1 = mulc(C,F0[0]);
323: R = cons([C1,F0[1]],R);
324: }
325: return reverse(R);
326: }
327:
328: def mulpf(F1,F2)
329: {
330: R = [];
331: N = length(PF_V);
332: for ( T = reverse(F1); T != []; T = cdr(T) ) {
333: T0 = car(T); C = T0[0]; Op = T0[1];
334: E = dp_etov(Op);
335: S = F2;
336: for ( I = 0; I < N; I++ )
337: for ( J = 0; J < E[I]; J++ ) S = muldpf(PF_V[I],S);
338: S = mulcpf(C,S);
339: R = addpf(R,S);
340: }
341: return S;
342: }
343:
344: def diffc(X,C)
345: {
346: R = [];
347: for ( T = C; T != []; T = cdr(T) ) {
348: T0 = car(T);
349: N = T0[0];
350: S = [];
351: for ( L = T0[1]; L != []; S = cons(car(L),S), L = cdr(L) ) {
352: L0 = car(L);
353: F = L0[0]; M = L0[1];
354: DF = diff(F,X);
355: if ( DF ) {
356: V = cons([F,M+1],cdr(L));
357: for ( U = S; U != []; U = cdr(U) ) V = cons(car(U),V);
358: R = addc([[-M*N*DF,V]],R);
359: }
360: }
361: }
362: return R;
363: }
364:
365: def muldpf(X,F)
366: {
367: R = [];
368: DX = dp_ptod(strtov("d"+rtostr(X)),PF_DV);
369: for ( T = F; T != []; T = cdr(T) ) {
370: T0 = car(T);
371: C = diffc(X,T0[0]);
372: if ( C != [] )
373: R = addpf(R,[[T0[0],T0[1]*DX],[C,T0[1]]]);
374: else
375: R = addpf(R,[[T0[0],T0[1]*DX]]);
376: }
377: return R;
378: }
379:
380: /* d^m*c = sum_{i=0}^m m!/i!/(m-i)!*c^(i)*d^(m-i) */
381: def muldmpf(X,M,F)
382: {
383: R = [];
384: DX = strtov("d"+rtostr(X));
385: FM = fac(M);
386: for ( T = F; T != []; T = cdr(T) ) {
387: T0 = car(T);
388: C = T0[0]; Op = T0[1];
389: U = [];
390: for ( I = 0; I <= M; I++ ) {
391: if ( C == [] ) break;
392: C0 = [[FM/fac(I)/fac(M-I),[]]];
393: U = cons([mulc(C0,C),dp_ptod(DX^(M-I),PF_DV)*Op],U);
394: C = diffc(X,C);
395: }
396: U = reverse(U);
397: R = addpf(U,R);
398: }
399: return R;
400: }
401:
402: def normalizef1(F)
403: {
404: if ( coef(F,1,var(F)) < 0 ) F = -F;
405: return F;
406: }
407:
408: def normalizec1(C)
409: {
410: N = C[0]; L = C[1];
411: S = 1;
412: R = [];
413: for ( ; L != []; L = cdr(L) ) {
414: L0 = car(L);
415: F = L0[0]; M = L0[1];
416: if ( coef(F,1,var(F)) < 0 ) {
417: F = -F;
418: if ( M%2 ) S = -S;
419: }
420: R = cons([F,M],R);
421: }
422: return [S*N,reverse(qsort(R,n_wishartd.compf1))];
423: }
424:
425: def normalizepf(F)
426: {
427: R = [];
428: for ( ; F != []; F = cdr(F) ) {
429: F0 = car(F);
430: C = map(normalizec1,F[0]);
431: R = cons([C,F[1]],R);
432: }
433: return reverse(R);
434: }
435:
436: def substc(C,Y2,Y1)
437: {
438: R = [];
439: for ( T = C; T != []; T = cdr(T) ) {
440: T0 = car(T); N = T0[0]; D = T0[1];
441: Z = subst_f(D,Y2,Y1);
442: R = addc([[Z[0]*N,Z[1]]],R);
443: }
444: return R;
445: }
446:
447: /* dy0 : dummy variable for dy */
448:
449: def wsetup(N)
450: {
451: V = [];
452: DV = [];
453: for ( I = N; I>= 0; I-- ) {
454: YI = strtov("y"+rtostr(I));
455: DYI = strtov("dy"+rtostr(I));
456: V = cons(YI,V);
457: DV = cons(DYI,DV);
458: }
459: PF_N = N;
460: PF_V = V;
461: PF_DV = DV;
462: ord(append(V,DV));
463: }
464:
465: def mtot(M,Dn)
466: {
467: D = dp_ptod(M,PF_DV);
468: F = fctr(Dn); C = car(F)[0];
469: Dn = reverse(qsort(cdr(F),n_wishartd.compf1));
470: return [[[dp_hc(D)/C,Dn]],dp_ht(D)];
471: }
472:
473: def wishartpf(N,I)
474: {
475: YI = PF_V[I]; DYI = PF_DV[I];
476: R = [mtot(DYI^2,1)];
477: R = addpf([mtot(c*DYI,YI)],R);
478: R = addpf([mtot(-DYI,1)],R);
479: for ( J = 1; J <= N; J++ ) {
480: if ( J == I ) continue;
481: YJ = PF_V[J]; DYJ = PF_DV[J];
482: R = addpf([mtot(1/2*DYI,YI-YJ)],R);
483: R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
484: R = addpf([mtot(-1/2*DYI,YI)],R);
485: R = addpf([mtot(1/2*DYJ,YI)],R);
486: }
487: R = addpf([mtot(-a,YI)],R);
488: return R;
489: }
490:
491: def degf(F,D)
492: {
493: for ( ; F != []; F = cdr(F) ) {
494: F0 = car(F);
495: if ( F0[0] == D ) return F0[1];
496: }
497: return 0;
498: }
499:
500: def removef(F,D)
501: {
502: R = [];
503: for ( ; F != []; F = cdr(F) ) {
504: F0 = car(F);
505: if ( F0[0] != D ) R = cons(F0,R);
506: }
507: return reverse(R);
508: }
509:
510: def subst_f(F,Y2,Y1)
511: {
512: R = [];
513: Sgn = 1;
514: for ( ; F != []; F = cdr(F) ) {
515: F0 = car(F);
516: T = subst(F0[0],Y2,Y1);
517: if ( coef(T,1,var(T)) < 0 ) {
518: T = -T;
519: if ( F0[1]%2 ) Sgn = -Sgn;
520: }
521: R = cons([T,F0[1]],R);
522: }
523: if ( R == [] ) return [Sgn,R];
524: R = qsort(R,n_wishartd.compf1);
525: S = [car(R)]; R = cdr(R);
526: while ( R != [] ) {
527: R0 = car(R);
528: S0 = car(S);
529: if ( R0[0] == S0[0] )
530: S = cons([S0[0],S0[1]+R0[1]],cdr(S));
531: else
532: S = cons(R0,S);
533: R = cdr(R);
534: }
535: return [Sgn,S];
536: }
537:
538: /* Y2 -> Y1 */
539: def lopitalpf(P,Y1,Y2)
540: {
541: // if ( !member(Y2,vars(P)) ) return P;
542: D = normalizef1(Y2-Y1);
543: if ( D == Y2-Y1 ) Sgn = 1;
544: else Sgn = -1;
545: DY2 = strtov("d"+rtostr(Y2));
546: R = [];
547: for ( T = reverse(P); T != []; T = cdr(T) ) {
548: T0 = car(T);
549: C = T0[0]; Op = T0[1];
550: for ( U = reverse(C); U != []; U = cdr(U) ) {
551: U0 = car(U);
552: Nm = U0[0]; Dn = U0[1];
553: K = degf(Dn,D);
554: if ( !K ) {
555: Z = subst_f(Dn,Y2,Y1);
556: Sgn1 = Z[0]; Dn1 = Z[1];
557: R = addpf([[[[Sgn1*Nm,Dn1]],Op]],R);
558: } else {
559: Dn1 = removef(Dn,D);
560: C1 = [[1,Dn1]];
561: Z = [];
562: for ( J = 0; J <= K; J++ ) {
563: B = fac(K)/fac(J)/fac(K-J);
564: C1s = substc(C1,Y2,Y1);
565: if ( C1s != [] ) {
566: W = [[C1s,dp_ptod(DY2^(K-J),PF_DV)*Op]];
567: W = mulcpf([[B,[]]],W);
568: Z = addpf(W,Z);
569: }
570: C1 = diffc(Y2,C1);
571: if ( C1 == [] ) break;
572: }
573: Z = mulcpf([[Sgn^K*Nm/fac(K),[]]],Z);
574: R = addpf(Z,R);
575: }
576: }
577: }
578: return R;
579: }
580:
581: def lopital_others_pf(P,L) {
582: Q = P;
583: for ( T = L; T != []; T = cdr(T) ) {
584: T0 = car(T);
585: I0 = T0[0]; I1 = T0[1]; I = I1-I0+1;
586: for ( M = I0; M <= I1; M++ ) {
587: Q = lopitalpf(Q,PF_V[I0],PF_V[M]);
588: }
589: Q = simppf(Q,I0,I);
590: Q = adjop(Q,I0,I);
591: }
592: return Q;
593: }
594:
595: def simpop(Op,I0,NV)
596: {
597: E = dp_etov(Op);
598: R = [];
599: for ( J = 0, I = I0; J < NV; I++, J++ ) R = cons(E[I],R);
600: R = reverse(qsort(R));
601: for ( J = 0, I = I0; J < NV; I++, J++ ) E[I] = R[J];
602: return dp_vtoe(E);
603: }
604:
605: def simppf(P,I0,NV)
606: {
607: R = [];
608: for ( T = P; T != []; T = cdr(T) ) {
609: T0 = car(T);
610: R = addpf([[T0[0],simpop(T0[1],I0,NV)]],R);
611: }
612: return R;
613: }
614:
615: /* simplify (dy1+...+dyn)^k */
616:
617: def simp_power1(K,I0,NV)
618: {
619: P = all_partition(K,NV);
620: M = map(coef_partition,P,K,NV,I0);
621: for ( R = 0, T = M; T != []; T = cdr(T) )
622: R += car(T);
623: return R;
624: }
625:
626: def indep_eqs(Eq)
627: {
628: M = length(Eq);
629: D = 0;
630: for ( I = 0; I < M; I++ )
631: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
632: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
633: Terms = vector(N);
634: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
635: A = matrix(M,N);
636: for ( I = 0; I < M; I++ ) {
637: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
638: T = car(H)[1];
639: for ( K = 0; K < N; K++ )
640: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
641: }
642: }
643: for ( I = 0; I < M; I++ ) {
644: Dn = 1;
645: for ( J = 0; J < N; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
646: for ( J = 0; J < N; J++ ) A[I][J] *= Dn;
647: }
648: R = indep_rows_mod(A,lprime(0));
649: if ( length(R) == N ) {
650: L = [];
651: for ( I = N-1; I >= 0; I-- )
652: L = cons(Eq[R[I]],L);
653: return L;
654: } else
655: return 0;
656: }
657:
658: def eliminatetop(Eq)
659: {
660: Eq = indep_eqs(Eq);
661: if ( !Eq ) return 0;
662:
663: M = length(Eq);
664: R = vector(M);
665: D = 0;
666: for ( I = 0; I < M; I++ )
667: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
668: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
669: if ( N != M )
670: return 0;
671: N2 = 2*N;
672: Terms = vector(N);
673: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
674: A = matrix(N,N2);
675: for ( I = 0; I < N; I++ ) {
676: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
677: T = car(H)[1];
678: for ( K = 0; K < N; K++ )
679: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
680: }
681: A[I][N+I] = 1;
682: }
683: for ( I = 0; I < N; I++ ) {
684: Dn = 1;
685: for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
686: for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
687: }
688: Ret = generic_gauss_elim(A);
689: Num = Ret[0]; Den = Ret[1]; Ind0 = Ret[2]; Ind1 = Ret[3];
690: if ( length(Ind0) != N || Ind0[N-1] != N-1 )
691: return 0;
692: R = [];
693: if ( Print ) print(["N=",N]);
694: for ( I = 0; I < N; I++ ) {
695: if ( Print ) print(["I=",I],2);
696: T = [];
697: for ( L = 0; L < N; L++ )
698: if ( Num[I][L] ) T = addpf(T,mulcpf([[Num[I][L]/Den,[]]],Eq[L][1]));
699: R = cons([Terms[I],T],R);
700: }
701: if ( Print ) print("");
702: R = qsort(R,n_wishartd.compt);
703: return reverse(R);
704: }
705:
706: def eliminatecheck(Eq,Top)
707: {
708: Eq = indep_eqs(Eq);
709: if ( !Eq ) return 0;
710:
711: M = length(Eq);
712: R = vector(M);
713: D = 0;
714: for ( I = 0; I < M; I++ )
715: for ( H = Eq[I][0]; H != []; H = cdr(H) ) D += car(H)[1];
716: for ( N = 0, T = D; T; T = dp_rest(T), N++ );
717: if ( N != M )
718: return 0;
719: N2 = 2*N;
720: Terms = vector(N);
721: for ( I = 0, T = D; T; T = dp_rest(T), I++ ) Terms[I] = dp_ht(T);
722: A = matrix(N,N2);
723: for ( I = 0; I < N; I++ ) {
724: for ( H = Eq[I][0]; H != []; H = cdr(H) ) {
725: T = car(H)[1];
726: for ( K = 0; K < N; K++ )
727: if ( Terms[K] == T ) A[I][K] = ctorat(car(H)[0]);
728: }
729: A[I][N+I] = 1;
730: }
731: for ( I = 0; I < N; I++ ) {
732: Dn = 1;
733: for ( J = 0; J < N2; J++ ) if ( A[I][J] ) Dn = ilcm(dn(A[I][J]),Dn);
734: for ( J = 0; J < N2; J++ ) A[I][J] *= Dn;
735: }
736: if ( Top ) {
737: for ( I = 0; I < N; I++ ) {
738: for ( T = J = 0; J < N; J++ ) if ( J != I ) T += abs(A[I][J]);
739: if ( abs(A[I][I]) < T )
740: if ( Print ) print(["ng",I]);
741: }
742: } else {
743: for ( I = 1; I < N; I++ ) {
744: for ( T = J = 0; J < N-2; J++ ) if ( J != I-1 ) T += abs(A[I][J]);
745: if ( abs(A[I][I-1]) < T )
746: if ( Print ) print(["ng",I]);
747: }
748: }
749:
750: Ret = generic_gauss_elim_mod(A,lprime(0));
751: Num = Ret[0]; Ind0 = Ret[1]; Ind1 = Ret[2];
752: if ( length(Ind0) != N || Ind0[N-1] != N-1 )
753: return 0;
754: else
755: return 1;
756: }
757:
758: def mul_lopitalpf(Z,N,K,I0,I1,E)
759: {
760: I = I1-I0+1;
761: for ( J = I0; J <= N; J++ )
762: for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
763: for ( J = I0+1; J <= I1; J++ ) {
764: Z = lopitalpf(Z,PF_V[I0],PF_V[J]);
765: }
766: S = simppf(Z,I0,I);
767: H = []; L = [];
768: for ( ; S != []; S = cdr(S) ) {
769: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
770: else L = cons(car(S),L);
771: }
772: return [reverse(H),reverse(L)];
773: }
774:
775: /* Blocks = [B1,B2,...] */
776: /* Bi=[s,e] ys=...=ye f */
777:
778: def diag0pf(N,Blocks)
779: {
780: Tlopital = Tred = Tlineq = 0;
781: wsetup(N);
782: P = vector(N+1);
783: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
784: Simp = [];
785: Done = [];
786: Len = length(Blocks);
787: for ( A = 0; A < Len; A++ ) {
788: B = Blocks[A];
789: Blocks0 = setminus(Blocks,[B]);
790: I0 = B[0];
791: I1 = B[1];
792: I = I1-I0+1;
793: for ( K = I0; K <= I1; K++ )
794: P[K] = lopital_others_pf(P[K],Blocks0);
795: Rall = [];
796: for ( K = I+1; K >= 2; K-- ) {
797: if ( Print ) print(["K=",K],2);
798: DK = simp_power1(K,I0,I);
799: Eq = [];
800: TlopitalK = 0;
801: for ( T = DK; T; T = dp_rest(T) ) {
802: E = dp_etov(T);
803: for ( U = I0; U <= I1; U++ )
804: if ( E[U] >= 2 )
805: break;
806: if ( U > I1 ) continue;
807: E[U] -= 2;
808: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
809: Eq = cons(Ret,Eq);
810: }
811: Eq = reverse(Eq);
812:
813: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
814: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
815:
816: DY = mtot(-dy0^K,1);
817: if ( K == I+1 ) {
818: EqTop = addpf([DY],Eq0);
819: } else {
820: H = []; L = [];
821: for ( S = Eq0 ; S != []; S = cdr(S) ) {
822: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
823: else L = cons(car(S),L);
824: }
825: L = reverse(L); H = reverse(H);
826: Eq = cons([H,addpf([DY],L)],Eq);
827: }
828: T0 = time();
829: R = eliminatetop(Eq);
830: if ( R )
831: Rall = cons(R,Rall);
832: else
833: return [];
834: T1 = time(); Tlineq += T1[0]-T0[0];
835: if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
836: }
837: Rall = reverse(Rall);
838:
839: /* EqTop : dyi0 -> dy -> dyi0 */
840: for ( T = Rall; T != []; T = cdr(T) ) {
841: if ( Print ) print(["red",length(T)+1],2);
842: T0 = time();
843: EqTop = reducepf(EqTop,car(T));
844: T1 = time(); Tred += T1[0]-T0[0];
845: if ( Print ) print([T1[0]-T0[0],"sec"]);
846: }
847: EqTop = adjop(EqTop,I0,I);
848: Done = cons(EqTop,Done);
849:
850: }
851: if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
852: Done = ltov(Done);
853: Len = length(Done);
854: for ( I = 0; I < Len; I++ ) {
855: for ( G = [], J = 0; J < Len; J++ )
856: if ( J != I ) G = cons(Done[J],G);
857: Done[I] = nfpf(Done[I],G);
858: }
859: return vtol(Done);
860: }
861:
862: def diagpf(N,Blocks)
863: "usage : n_wishartd.diagpf(M,[[S1,E1],[S2,E2],...]),\n where M is the number of variables and [Si,Ei] is a diagonal block and S(i+1)=Ei + 1. For example, n_wishartd.diagpf(10,[[1,9],[10,10]) returns a system of PDEs satisfied by 1F1(y1,...,y1,y10)."
864: {
865: Tlopital = Tred = Tlineq = 0;
866: wsetup(N);
867: P = vector(N+1);
868: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
869: Simp = [];
870: Done = [];
871: Len = length(Blocks);
872: for ( A = 0; A < Len; A++ ) {
873: B = Blocks[A];
874: Blocks0 = setminus(Blocks,[B]);
875: I0 = B[0];
876: I1 = B[1];
877: I = I1-I0+1;
878: for ( K = I0; K <= I1; K++ )
879: P[K] = lopital_others_pf(P[K],Blocks0);
880: Rall = [];
881: for ( K = I+1; K >= 2; K-- ) {
882: if ( Print ) print(["K=",K],2);
883: DK = simp_power1(K,I0,I);
884: Eq = [];
885: TlopitalK = 0;
886: for ( T = DK; T; T = dp_rest(T) ) {
887: E = dp_etov(T);
888: for ( U = I1; U >= I0; U-- )
889: if ( E[U] >= 2 )
890: break;
891: if ( U < I0 ) continue;
892: E[U] -= 2;
893: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E);
894: Eq = cons(Ret,Eq);
895: }
896: Eq = reverse(Eq);
897:
898: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
899: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
900:
901: DY = mtot(-dy0^K,1);
902: if ( K == I+1 ) {
903: EqTop = addpf([DY],Eq0);
904: } else {
905: H = []; L = [];
906: for ( S = Eq0 ; S != []; S = cdr(S) ) {
907: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
908: else L = cons(car(S),L);
909: }
910: L = reverse(L); H = reverse(H);
911: Eq = cons([H,addpf([DY],L)],Eq);
912: }
913: T0 = time();
914: R = eliminatetop(Eq);
915: if ( R )
916: Rall = cons(R,Rall);
917: else {
918: if ( Print ) print(["lineq retry.."],2);
919: Eq1 = [];
920: Terms = [];
921: for ( T = dp_rest(DK); T; T = dp_rest(T) )
922: Terms = cons(dp_ht(T),Terms);
923: while ( Terms != [] ) {
924: for ( Q = 0; Terms != [] && Q < 10; Q++, Terms = cdr(Terms) ) {
925: E = dp_etov(car(Terms));
926: if ( E[I0] >= 2 ) {
927: E[I0] -= 2;
928: Ret = mul_lopitalpf(P[I0],N,K,I0,I1,E);
929: Eq1 = cons(Ret,Eq1);
930: }
931: }
932: Eq = append(Eq,Eq1);
933: R = eliminatetop(Eq);
934: if ( R ) break;
935: }
936: if ( !R )
937: error("diagpf : failed to find a solvable linear system");
938: Rall = cons(R,Rall);
939: }
940: T1 = time(); Tlineq += T1[0]-T0[0];
941: if ( Print ) print(["lineq",T1[0]-T0[0],"sec"]);
942: }
943: Rall = reverse(Rall);
944:
945: /* EqTop : dyi0 -> dy -> dyi0 */
946: for ( T = Rall; T != []; T = cdr(T) ) {
947: if ( Print ) print(["red",length(T)+1],2);
948: T0 = time();
949: EqTop = reducepf(EqTop,car(T));
950: T1 = time(); Tred += T1[0]-T0[0];
951: if ( Print ) print([T1[0]-T0[0],"sec"]);
952: }
953: EqTop = adjop(EqTop,I0,I);
954: Done = cons(EqTop,Done);
955:
956: }
957: if ( Print ) print(["lopital",Tlopital,"lineq",Tlineq,"red",Tred]);
958: Done = ltov(Done);
959: Len = length(Done);
960: for ( I = 0; I < Len; I++ ) {
961: for ( G = [], J = 0; J < Len; J++ )
962: if ( J != I ) G = cons(Done[J],G);
963: Done[I] = nfpf(Done[I],G);
964: }
965: return vtol(Done);
966: }
967:
968: def diagcheck(N)
969: {
970: Tmuld = Tlopital = 0;
971: MHist = [];
972: wsetup(N);
973: P = vector(N+1);
974: for ( J = 1; J <= N; J++ ) P[J] = wishartpf(N,J);
975: Simp = [];
976: Done = [];
977: B = [1,N];
978: I0 = B[0];
979: I1 = B[1];
980: I = I1-I0+1;
981: for ( K = 2; K <= I+1; K++ ) {
982: if ( Print ) print(["K=",K],2);
983: DK = simp_power1(K,I0,I);
984: Eq = [];
985: TlopitalK = 0;
986: for ( T = DK; T; T = dp_rest(T) ) {
987: E = dp_etov(T);
988: for ( U = I0; U <= I1; U++ )
989: if ( E[U] >= 2 )
990: break;
991: if ( U > I1 ) continue;
992: E[U] -= 2;
993: Ret = mul_lopitalpf2(P,U,N,K,I0,I1,E|high=1);
994: Eq = cons(Ret,Eq);
995: }
996: Eq = reverse(Eq);
997:
998: for ( Eq0 = [], T = DK; T; T = dp_rest(T) )
999: Eq0 = addpf([[[[dp_hc(T),[]]],dp_ht(T)]],Eq0);
1000:
1001: DY = mtot(-dy0^K,1);
1002: if ( K == I+1 ) {
1003: EqTop = addpf([DY],Eq0);
1004: } else {
1005: H = []; L = [];
1006: for ( S = Eq0 ; S != []; S = cdr(S) ) {
1007: if ( dp_td(car(S)[1]) == K ) H = cons(car(S),H);
1008: else L = cons(car(S),L);
1009: }
1010: L = reverse(L); H = reverse(H);
1011: Eq = cons([H,addpf([DY],L)],Eq);
1012: }
1013: R = eliminatecheck(Eq,K==I+1);
1014: if ( !R )
1015: return 0;
1016: }
1017: if ( Print ) print(["muld",Tmuld,"lopital",Tlopital]);
1018: return 1;
1019: }
1020:
1021: /* dyi -> dyi/K, dy->dyi */
1022: def adjop1(T,I,K)
1023: {
1024: C = T[0];
1025: Op = dp_etov(T[1]);
1026: if ( Op[I] ) C = mulc([[1/K,[]]],C);
1027: Op[I] += Op[0];
1028: Op[0] = 0;
1029: return [C,dp_vtoe(Op)];
1030: }
1031:
1032: def adjop(P,I,K)
1033: {
1034: P = map(adjop1,P,I,K);
1035: P = qsort(P,n_wishartd.compt);
1036: return reverse(P);
1037: }
1038:
1039: def dnc(C)
1040: {
1041: D = 1;
1042: for ( T = C; T != []; T = cdr(T) ) {
1043: T0 = car(T); N = T0[0];
1044: M = sdiv(ptozp(N),N);
1045: D = ilcm(D,M);
1046: }
1047: return D;
1048: }
1049:
1050: def pftozpf(F)
1051: {
1052: D = 1;
1053: for ( T = F; T != []; T = cdr(T) ) {
1054: T0 = car(T);
1055: D = ilcm(D,dnc(T0[0]));
1056: }
1057: return [D,mulcpf([[D,[]]],F)];
1058: }
1059:
1060: def zpftopf(F)
1061: {
1062: return mulcpf([[1/F[0],[]]],F[1]);
1063: }
1064:
1065: def addzpf(F1,F2)
1066: {
1067: D1 = F1[0]; D2 = F2[0];
1068: L = ilcm(D1,D2); C1 = L/D1; C2 = L/D2;
1069: N = addpf(mulcpf([[L/D1,[]]],F1[1]),mulcpf([[L/D2,[]]],F2[1]));
1070: return [L,N];
1071: }
1072:
1073: /* R : reducers */
1074: def reducepf(Eq,R)
1075: {
1076: Ret = pftozpf(Eq); Eq = Ret[1]; Dn = Ret[0];
1077: S = [1,[]];
1078: for ( T = R, U = []; T != []; T = cdr(T) )
1079: U = cons([car(T)[0],pftozpf(car(T)[1])],U);
1080: R = reverse(U);
1081: for ( T = reverse(Eq); T != []; T = cdr(T) ) {
1082: T0 = car(T); C = T0[0]; Op = T0[1];
1083: for ( U = (R); U != []; U = cdr(U) ) {
1084: U0 = car(U);
1085: if ( dp_redble(Op,U0[0]) ) {
1086: E = dp_etov(dp_subd(Op,U0[0]));
1087: Red = U0[1];
1088: N = length(E);
1089: for ( J = 1; J < N; J++ )
1090: for ( L = 0; L < E[J]; L++ ) Red = [Red[0],muldpf(PF_V[J],Red[1])];
1091: Red = [Red[0],mulcpf(C,Red[1])];
1092: Red = [Red[0],mulcpf([[-1,[]]],Red[1])];
1093: S = addzpf(S,Red);
1094: break;
1095: }
1096: }
1097: if ( U == [] ) S = addzpf([1,[T0]],S);
1098: }
1099: S = [S[0]*Dn,S[1]];
1100: return zpftopf(S);
1101: }
1102:
1103: def reduceotherspf(Eq,P,I1,N)
1104: {
1105: S = [];
1106: Z = Eq;
1107: while ( Z != [] ) {
1108: T0 = car(Z); C = T0[0]; Op = T0[1];
1109: for ( I = I1; I <= N; I++ ) {
1110: U0 = P[I];
1111: M = car(U0)[0][0][0];
1112: H = car(U0)[1];
1113: if ( dp_redble(Op,H) ) {
1114: E = dp_etov(dp_subd(Op,H));
1115: Red = U0;
1116: Len = length(E);
1117: for ( J = 0; J < Len; J++ )
1118: for ( L = 0; L < E[J]; L++ ) Red = muldpf(PF_V[J],Red);
1119: C1 = mulc([[1/M,[]]],C);
1120: Red = mulcpf(C1,Red);
1121: Z = subpf(Z,Red);
1122: break;
1123: }
1124: }
1125: if ( I > N ) {
1126: S = cons(T0,S);
1127: Z = cdr(Z);
1128: }
1129: }
1130: return reverse(S);
1131: }
1132:
1133: def print_f(F)
1134: {
1135: for ( ; F != []; F = cdr(F) ) {
1136: F0 = car(F);
1137: print("(",0); print(F0[0],0); print(")",0);
1138: if ( F0[1] > 1 ) {
1139: print("^",0); print(F0[1],0);
1140: }
1141: if ( cdr(F) != [] ) print("*",0);
1142: }
1143: }
1144:
1145: def printc(C)
1146: {
1147: print("(",0);
1148: for ( ; C != []; C = cdr(C) ) {
1149: C0 = car(C);
1150: print("(",0); print(C0[0],0); print(")",0);
1151: if ( C0[1] != [] ) {
1152: print("/(",0);
1153: print_f(C0[1]);
1154: print(")",0);
1155: }
1156: if ( cdr(C) != [] ) print("+",0);
1157: }
1158: print(")",0);
1159: }
1160:
1161: def printt(T)
1162: {
1163: printc(T[0]); print("*",0); print(dp_dtop(T[1],PF_DV),0);
1164: }
1165:
1166: def printpf(F)
1167: {
1168: for ( ; F != []; F = cdr(F) ) {
1169: printt(car(F));
1170: if ( cdr(F) != [] ) print("+",0);
1171: }
1172: print("");
1173: }
1174:
1175: def ftop(F)
1176: {
1177: P = 1;
1178: for ( ; F != []; F = cdr(F) ) {
1179: F0 = car(F);
1180: P *= F0[0]^F0[1];
1181: }
1182: return P;
1183: }
1184:
1185: def pftord(F)
1186: {
1187: R = [];
1188: for ( T = F; T != []; T = cdr(T) ) {
1189: T0 = car(T); C = T0[0]; Op = T0[1];
1190: R = cons([ctord(C),Op],R);
1191: }
1192: return reverse(R);
1193: }
1194:
1195: def pftop(F)
1196: {
1197: R = pftord(F);
1198: Op = F[0][1];
1199: N = length(dp_etov(Op));
1200: DY = [];
1201: for ( I = N; I >= 0; I-- ) DY = cons(strtov("dy"+rtostr(I)),DY);
1202: Lcm = [];
1203: for ( T = R; T != []; T = cdr(T) )
1204: Lcm = lcmf(Lcm,T[0][0][1]);
1205: S = 0;
1206: for ( T = R; T != []; T = cdr(T) ) {
1207: N = T[0][0][0]*ftop(divf(Lcm,T[0][0][1]));
1208: Op = dp_dtop(T[0][1],DY);
1209: S += N*Op;
1210: }
1211: return S;
1212: }
1213:
1214: def ctord(C)
1215: {
1216: N = 0; D = [];
1217: for ( T = reverse(C); T != []; T = cdr(T) ) {
1218: T0 = car(T); N0 = T0[0]; D0 = T0[1];
1219: L = lcmf(D,D0);
1220: /* N/D+N0/D0 = (N*(L/D)+N0*(L/D0))/L */
1221: N = N*ftop(divf(L,D))+N0*ftop(divf(L,D0));
1222: D = L;
1223: }
1224: R = [];
1225: for ( T = D; T != []; T = cdr(T) ) {
1226: T0 = car(T); F = T0[0]; M = T0[1];
1227: while ( M > 0 && (Q = tdiv(N,F)) ) {
1228: N = Q;
1229: M--;
1230: }
1231: if ( M ) R = cons([F,M],R);
1232: }
1233: D = reverse(R);
1234: return [N,D];
1235: }
1236:
1237:
1238: def comppfrd(P,R)
1239: {
1240: P = qsort(P,n_wishartd.compt); P=reverse(P);
1241: R = qsort(R,n_wishartd.compt); R=reverse(R);
1242: if ( length(P) != length(R) ) return 0;
1243: for ( ; P != []; P = cdr(P), R = cdr(R) ) {
1244: P0 = car(P); R0 = car(R);
1245: C0 = ctord(P0[0]); C1 = R0[0];
1246: D0 = ftop(C0[1]); D1 = ftop(C1[1]);
1247: if ( (D0 == D1) && C0[0] == -C1[0] ) continue;
1248: if ( (D0 == -D1) && C0[0] == C1[0] ) continue;
1249: return 0;
1250: }
1251: return 1;
1252: }
1253:
1254: def multpf(T,F)
1255: {
1256: E = dp_etov(T[1]);
1257: N = length(E);
1258: Z = F;
1259: for ( J = 1; J < N; J++ )
1260: for ( L = 0; L < E[J]; L++ ) Z = muldpf(PF_V[J],Z);
1261: Z = mulcpf(T[0],Z);
1262: return Z;
1263: }
1264:
1265: def spolypf(F,G)
1266: {
1267: F0 = car(F); G0 = car(G);
1268: DF = F0[1]; DG = G0[1];
1269: L = dp_lcm(DF,DG);
1270: F1 = multpf([F0[0],dp_subd(L,DF)],F);
1271: G1 = multpf([G0[0],dp_subd(L,DG)],G);
1272: S = subpf(F1,G1);
1273: return S;
1274: }
1275:
1276: def nfpf(F,G)
1277: {
1278: Z = F;
1279: R = [];
1280: while ( Z != [] ) {
1281: Z0 = car(Z); C = Z0[0]; D = Z0[1];
1282: for ( T = G; T != []; T = cdr(T) ) {
1283: T0 = car(T);
1284: if ( dp_redble(D,T0[0][1]) ) break;
1285: }
1286: if ( T != [] ) {
1287: CG = ctorat(T0[0][0]);
1288: C = mulc(C,[[-1/CG,[]]]);
1289: S = multpf([C,dp_subd(D,T0[0][1])],T0);
1290: Z = addpf(Z,S);
1291: } else {
1292: R = cons(Z0,R);
1293: Z = cdr(Z);
1294: }
1295: }
1296: S = [];
1297: for ( T = R; T != []; T = cdr(T) ) {
1298: U = ctord(T[0][0]);
1299: if ( U[0] )
1300: S = cons(T[0],S);
1301: }
1302: return S;
1303: }
1304:
1305: def nfpf0(F,G)
1306: {
1307: Z = F;
1308: R = [];
1309: while ( Z != [] ) {
1310: Z0 = car(Z); C = Z0[0]; D = Z0[1];
1311: for ( T = G; T != []; T = cdr(T) ) {
1312: T0 = car(T);
1313: if ( dp_redble(D,T0[0][1]) ) break;
1314: }
1315: if ( T != [] ) {
1316: CG = ctorat(T0[0][0]);
1317: C = mulc(C,[[-1/CG,[]]]);
1318: S = multpf([C,dp_subd(D,T0[0][1])],T0);
1319: Z = addpf(Z,S);
1320: } else {
1321: R = cons(Z0,R);
1322: Z = cdr(Z);
1323: }
1324: }
1325: S = [];
1326: for ( T = R; T != []; T = cdr(T) ) {
1327: U = ctord(T[0][0]);
1328: if ( U[0] )
1329: S = cons(T[0],S);
1330: }
1331: return S;
1332: }
1333:
1334: def pftoeuler(F)
1335: {
1336: VDV = [y1,dy1];
1337: P = pftop(F);
1338: Z = dp_ptod(P,VDV);
1339: E = dp_etov(dp_ht(Z));
1340: S = E[1]-E[0];
1341: if ( S < 0 )
1342: error("pftoeuler : invalid input");
1343: P *= y1^S;
1344: Z = dp_ptod(P,VDV);
1345: E = dp_etov(dp_ht(Z));
1346: D = E[0];
1347: C = vector(D+1);
1348: C[0] = 1;
1349: for ( I = 1; I <= D; I++ )
1350: C[I] = C[I-1]*(t-I+1);
1351: R = 0;
1352: for ( T = Z; T; T = dp_rest(T) ) {
1353: E = dp_etov(dp_ht(T));
1354: S = E[0]-E[1];
1355: if ( S < 0 )
1356: error("pftoeuler : invalid input");
1357: R += dp_hc(T)*y1^S*C[E[1]];
1358: }
1359: return R;
1360: }
1361:
1362: def gammam(M,A)
1363: {
1364: R = 1;
1365: for ( I = 1; I <= M; I++ )
1366: R *= mpfr_gamma(A-(I-1)/2);
1367: R *= eval(@pi^(1/4*M*(M-1)));
1368: return R;
1369: }
1370:
1371: def genprc(M,N,PL,SL,X)
1372: {
1373: A = (M+1)/2; C = (N+M+1)/2;
1374: DetS = 1;
1375: TrIS = 0;
1376: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
1377: DetS *= car(S)^car(T);
1378: TrIS += car(T)/car(S);
1379: }
1380: C = 1/eval(DetS^(1/2*N)*2^(1/2*N*M));
1381: return C;
1382: }
1383:
1384: /*
1385: * PL=[P1,P2,...]: sizes of blocks
1386: * SL=[S1,S2,...]: [S1xP1,S2xP2,..]
1387: */
1388:
1389: def prob_by_hgm(M,N,PL,SL,X)
1390: "usage : n_wishartd.prob_by_hgm(M,N,[P1,P2,...],[S1,S2,...]|options) where M is the number of variables, N is the degrees of freedom, Pi is the size of i-th block and Si is the value of i-th (repeated) eigenvalue of Sigma. options : rk5=1 => use 5-th order Runge-Kutta method (default : rk4) double=1 => use machine double float (default : MPFR) step=k => set the number of steps=k (default : 10^4) init=x => set the maximal coordinate of the initial point=x (default : 1). For example, n_wishartd.prob_by_hgm(3,10,[2,1],[1/10,1],10|rk5=1,double=1) computes Pr[l1<10|diag(1/10,1/10,1)] by RK5 and machine double float."
1391: {
1392: A = (M+1)/2; C = (N+M+1)/2;
1393:
1394: B = []; V = []; Beta = []; SBeta = 0;
1395: Prev = 1;
1396: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
1397: V = cons(strtov("y"+rtostr(Prev)),V);
1398: B = cons([Prev,Prev+car(T)-1],B);
1399: Prev += car(T);
1400: Beta = cons(1/(2*car(S)),Beta);
1401: SBeta += car(T)/(2*car(S));
1402: }
1403: if ( Prev != M+1 )
1404: error("invalid block specification");
1405: B = reverse(B); V = reverse(V);
1406: Beta = reverse(Beta);
1407:
1408: T0 = time();
1409: if ( type(Z=getopt(eq)) == -1 )
1410: Z = diagpf(M,B);
1411: T1 = time(); Tdiag = (T1[0]-T0[0])+(T1[1]-T0[1]);
1412:
1413: if ( type(Step=getopt(step)) == -1 )
1414: Step = 10000;
1415:
1416: if ( type(Double=getopt(double)) == -1 )
1417: Double = 0;
1418:
1419: Z = subst(Z,a,A,c,C);
1420: Init=getopt(init);
1421: Rk5 = getopt(rk5);
1422: if ( type(Rk5) == -1 ) Rk5 = 0;
1423: F = genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,Rk5,Step|double=Double,init=Init);
1424: if ( Print ) print(["diag",Tdiag,"pfaffian",Tp,"ps",Tps,"rk",Trk]);
1425: return F;
1426: }
1427:
1428: def prob_by_ps(M,N,PL,SL,X)
1429: {
1430: A = (M+1)/2; C = (N+M+1)/2;
1431:
1432: B = []; V = []; Beta = []; SBeta = 0;
1433: Prev = 1;
1434: for ( T = PL, S = SL; T != []; T = cdr(T), S = cdr(S) ) {
1435: V = cons(strtov("y"+rtostr(Prev)),V);
1436: B = cons([Prev,Prev+car(T)-1],B);
1437: Prev += car(T);
1438: Beta = cons(1/(2*car(S)),Beta);
1439: SBeta += car(T)/(2*car(S));
1440: }
1441: if ( Prev != M+1 )
1442: error("invalid block specification");
1443: B = reverse(B); V = reverse(V);
1444: Beta = reverse(Beta);
1445:
1446: if ( type(Z=getopt(eq)) == -1 )
1447: Z = diagpf(M,B);
1448:
1449: Z = subst(Z,a,A,c,C);
1450: GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2);
1451: C = GM*genprc(M,N,PL,SL,X);
1452:
1453: Beta = ltov(Beta)*eval(exp(0));
1454: X *= eval(exp(0));
1455: L = psvalue(Z,V,Beta*X); PS = L[0];
1456: MN2 = M*N/2;
1457: ExpF = eval(X^(MN2)*exp(-SBeta*X));
1458: return C*L[1]*ExpF;
1459: }
1460:
1461: def ps(Z,V,TD)
1462: {
1463: Tact = Tgr = Tactmul = 0;
1464: G = map(ptozp,map(pftop,Z));
1465: M = length(V);
1466: N = length(G);
1467: for ( I = 0, DV = []; I < M; I++ )
1468: DV = cons(strtov("d"+rtostr(V[I])),DV);
1469: DV = reverse(DV);
1470: VDV = append(V,DV);
1471: DG = map(decomp_op,G,V,DV);
1472: F = [1];
1473: R = 1;
1474: Den = 1;
1475: for ( I = 1 ; I <= TD; I++ ) {
1476: L = create_homo(V,I);
1477: FI = L[0]; CV = L[1];
1478: Eq = [];
1479: for ( K = 0; K < N; K++ ) {
1480: P = DG[K]; Len = length(P);
1481: /* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
1482: T0=time();
1483: Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
1484: for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
1485: Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
1486: T1=time(); Tact += T1[0]-T0[0];
1487: if ( Lhs ) {
1488: for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
1489: Eq = cons(dp_hc(T),Eq);
1490: }
1491: }
1492: }
1493: T0=time();
1494: #if 0
1495: Sol = solve_leq(Eq,CV);
1496: #else
1497: Sol = nd_f4(Eq,CV,0,0);
1498: #endif
1499: L = p_true_nf(FI,Sol,CV,0);
1500: Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
1501: FI = L[0]*(Den1/L[1]);
1502: T1=time(); Tgr += T1[0]-T0[0];
1503:
1504: MJ = Den1/Den; Den = Den1;
1505: for ( S = [], T = F; T != []; T = cdr(T) )
1506: S = cons(MJ*car(T),S);
1507: F = cons(FI,reverse(S));
1508: R = R*MJ+car(F);
1509: }
1510: return R/Den;
1511: }
1512:
1513: def msubst(F,V,X)
1514: {
1515: M = length(V);
1516: for ( J = 0, F0 = F*eval(exp(0)); J < M; J++ )
1517: F0 = subst(F0,V[J],X[J]);
1518: return F0;
1519: }
1520:
1521: def psvalue(Z,V,X)
1522: {
1523: Tact = Tgr = Tactmul = 0;
1524: G = map(ptozp,map(pftop,Z));
1525: M = length(V);
1526: N = length(G);
1527: for ( I = 0, DV = []; I < M; I++ )
1528: DV = cons(strtov("d"+rtostr(V[I])),DV);
1529: DV = reverse(DV);
1530: VDV = append(V,DV);
1531: DG = map(decomp_op,G,V,DV);
1532: Prev = getopt(prev);
1533: if ( type(Prev) == -1 ) {
1534: F = [1];
1535: R = 1;
1536: Val = eval(exp(0));
1537: Den = 1;
1538: I = 1;
1539: } else {
1540: /* Prev = [R/Den,Val1,F,Den] */
1541: Den = Prev[3];
1542: F = Prev[2];
1543: R = Prev[0]*Den;
1544: I = length(F);
1545: Val = msubst(Prev[0],V,X);
1546: ValI = msubst(F[0],V,X)/Den;
1547: if ( Val-ValI == Val )
1548: return [Prev[0],Val,F,Den];
1549: }
1550: for ( ; ; I++ ) {
1551: L = create_homo(V,I);
1552: FI = L[0]; CV = L[1];
1553: Eq = [];
1554: for ( K = 0; K < N; K++ ) {
1555: P = DG[K]; Len = length(P);
1556: /* P(0)F(I)*Den + (P(1)F(I-1)+...+P(I)F(0)) = 0 */
1557: T0=time();
1558: Lhs = dp_dtop(dp_weyl_act(P[0],dp_ptod(FI,V)),V)*Den;
1559: for ( T = F, J = 1; T != [] && J < Len; T = cdr(T), J++ )
1560: Lhs += dp_dtop(dp_weyl_act(P[J],dp_ptod(car(T),V)),V);
1561: T1=time(); Tact += T1[0]-T0[0];
1562: if ( Lhs ) {
1563: for ( T = dp_ptod(Lhs,V); T; T = dp_rest(T) ) {
1564: Eq = cons(dp_hc(T),Eq);
1565: }
1566: }
1567: }
1568: T0=time();
1569: Sol = solve_leq(Eq,CV);
1570: L = p_true_nf(FI,Sol,CV,0);
1571: Den1 = ilcm(Den,L[1]); MI = Den1/L[1];
1572: FI = L[0]*(Den1/L[1]);
1573: T1=time(); Tgr += T1[0]-T0[0];
1574:
1575: MJ = Den1/Den; Den = Den1;
1576: for ( S = [], T = F; T != []; T = cdr(T) )
1577: S = cons(MJ*car(T),S);
1578: F = cons(FI,reverse(S));
1579: R = R*MJ+car(F);
1580: F0 = msubst(FI,V,X)/Den;
1581: Val1 = Val+F0;
1582: if ( Val1 == Val ) {
1583: if ( Print ) print(["I=",I,"act",Tact,"actmul",Tactmul,"gr",Tgr]);
1584: delete_uc();
1585: return [R/Den,Val1,F,Den];
1586: } else {
1587: if ( Print ) print([F0],2);
1588: Val = Val1;
1589: }
1590: }
1591: }
1592:
1593: /* P -> [P0,P1,...] Pi = y^a dy^b, |a|-|b|=i */
1594:
1595: def decomp_op(P,V,DV)
1596: {
1597: M = length(V);
1598: VDV = append(V,DV);
1599: D = dp_ptod(P,VDV);
1600: for ( I = 0, T = D; T; T = dp_rest(T), I++ );
1601: for ( T = D, NotSet = 1; T; T = dp_rest(T) ) {
1602: E = dp_etov(T);
1603: for ( K = 0, J = 0; J < M; J++ )
1604: K += E[J]-E[M+J];
1605: if ( NotSet ) {
1606: Max = Min = K; NotSet = 0;
1607: } else if ( K > Max ) Max = K;
1608: else if ( K < Min ) Min = K;
1609: }
1610: W = vector(Max-Min+1);
1611: for ( T = D; T; T = dp_rest(T) ) {
1612: E = dp_etov(T);
1613: for ( K = 0, J = 0; J < M; J++ )
1614: K += E[J]-E[M+J];
1615: W[K-Min] += dp_hm(T);
1616: }
1617: W = map(dp_ptod,map(dp_dtop,W,VDV),DV);
1618: return W;
1619: }
1620:
1621: def create_homo(V,TD)
1622: {
1623: M = length(V);
1624: for ( S = 0, I = 0; I < M; I++ ) S += V[I];
1625: Tmp = 0;
1626: Nc = 0;
1627: for ( RI = 0, D = dp_ptod(S^TD,V); D; D = dp_rest(D), Nc++ ) {
1628: E = dp_etov(D);
1629: T = uc();
1630: U = dp_dtop(dp_ht(D),V);
1631: RI += T*U;
1632: Tmp += T;
1633: }
1634: CV = vector(Nc);
1635: for ( I = 0; I < Nc; I++ ) {
1636: CV[I] = var(Tmp); Tmp -= CV[I];
1637: }
1638: return [RI,vtol(CV)];
1639: }
1640:
1641: def act_op(P,V,F)
1642: {
1643: N = length(V);
1644: for ( R = 0; P; P = dp_rest(P) ) {
1645: C = dp_hc(P); E = dp_etov(P);
1646: for ( I = 0, S = F; I < N; I++ )
1647: for ( J = 0; J < E[I]; J++ ) S = diff(S,V[I]);
1648: T0 = time();
1649: R += C*S;
1650: T1 = time(); Tactmul += T1[0]-T0[0];
1651: }
1652: return R;
1653: }
1654:
1655: def gen_basis(D,DV)
1656: {
1657: N = length(D);
1658: B = [];
1659: E = vector(N);
1660: while ( 1 ) {
1661: B = cons(dp_dtop(dp_vtoe(E),DV),B);
1662: for ( I = N-1; I >= 0; I-- )
1663: if ( E[I] < D[I]-1 ) break;
1664: if ( I < 0 ) return reverse(B);
1665: E[I]++;
1666: for ( J = I+1; J < N; J++ ) E[J] = 0;
1667: }
1668: }
1669:
1670: def pfaffian(Z)
1671: {
1672: N = length(Z);
1673: D = vector(N);
1674: Mat = vector(N);
1675: V = vars(Z);
1676: DV = vector(N);
1677: for ( I = 0; I < N; I++ )
1678: DV[I] = strtov("d"+rtostr(V[I]));
1679: DV = vtol(DV);
1680: for ( I = 0; I < N; I++ ) {
1681: ZI = Z[I];
1682: HI = ZI[0][1];
1683: DI = dp_dtop(HI,PF_DV);
1684: for ( J = 0; J < N; J++ )
1685: if ( var(DI) == DV[J] )
1686: break;
1687: D[J] = deg(DI,DV[J]);
1688: }
1689: B = gen_basis(D,DV);
1690: Dim = length(B);
1691: Hist = [];
1692: for ( I = 0; I < N; I++ ) {
1693: if ( Print ) print(["I=",I]);
1694: A = matrix(Dim,Dim);
1695: for ( K = 0; K < Dim; K++ )
1696: for ( L = 0; L < Dim; L++ )
1697: A[K][L] = [];
1698: for ( J = 0; J < Dim; J++ ) {
1699: if ( Print ) print(["J=",J],2);
1700: Mon = DV[I]*B[J];
1701: for ( T = Hist; T != []; T = cdr(T) )
1702: if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
1703: if ( (T != []) ) {
1704: for ( L = 0; L < N; L++ )
1705: if ( Q == DV[L] ) break;
1706: Red1 = muldpf(V[L],car(T)[1]);
1707: Red = nfpf0(Red1,Z);
1708: } else
1709: Red = nfpf0([mtot(Mon,1)],Z);
1710: Hist = cons([Mon,Red],Hist);
1711: for ( T = Red; T != []; T = cdr(T) ) {
1712: T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
1713: for ( K = 0; K < Dim; K++ )
1714: if ( B[K] == Op ) break;
1715: if ( K == Dim )
1716: error("afo");
1717: else
1718: A[J][K] = C;
1719: }
1720: }
1721: if ( Print ) print("");
1722: A = map(ctord,A);
1723: Lcm = [];
1724: for ( K = 0; K < Dim; K++ )
1725: for ( L = 0; L < Dim; L++ )
1726: Lcm = lcmf(Lcm,A[K][L][1]);
1727: for ( K = 0; K < Dim; K++ )
1728: for ( L = 0; L < Dim; L++ ) {
1729: Num = A[K][L][0]; Den = A[K][L][1];
1730: A[K][L] = Num*ftorat(divf(Lcm,Den));
1731: }
1732: Lcm = ftorat(Lcm);
1733: if ( !Lcm ) Lcm = 1;
1734: Mat[I] = [A,Lcm];
1735: }
1736: return [Mat,B];
1737: }
1738:
1739: def pfaffian2(Z,V,Beta,XV)
1740: {
1741: N = length(Z);
1742: D = vector(N);
1743: Mat = vector(N);
1744: DV = vector(N);
1745: for ( I = 0; I < N; I++ )
1746: DV[I] = strtov("d"+rtostr(V[I]));
1747: DV = vtol(DV);
1748: for ( I = 0; I < N; I++ ) {
1749: ZI = Z[I];
1750: HI = ZI[0][1];
1751: DI = dp_dtop(HI,PF_DV);
1752: for ( J = 0; J < N; J++ )
1753: if ( var(DI) == DV[J] )
1754: break;
1755: D[J] = deg(DI,DV[J]);
1756: }
1757: B = gen_basis(D,DV);
1758: Dim = length(B);
1759: Hist = [];
1760: for ( I = 0; I < N; I++ ) {
1761: if ( Print ) print(["I=",I]);
1762: A = matrix(Dim,Dim);
1763: for ( K = 0; K < Dim; K++ )
1764: for ( L = 0; L < Dim; L++ )
1765: A[K][L] = [];
1766: for ( J = 0; J < Dim; J++ ) {
1767: if ( Print ) print(["J=",J],2);
1768: Mon = DV[I]*B[J];
1769: for ( T = Hist; T != []; T = cdr(T) )
1770: if ( (Q = tdiv(Mon,car(T)[0])) && Q==var(Q) ) break;
1771: if ( (T != []) ) {
1772: for ( L = 0; L < N; L++ )
1773: if ( Q == DV[L] ) break;
1774: Red1 = muldpf(V[L],car(T)[1]);
1775: Red = nfpf0(Red1,Z);
1776: } else
1777: Red = nfpf0([mtot(Mon,1)],Z);
1778: Hist = cons([Mon,Red],Hist);
1779: for ( T = Red; T != []; T = cdr(T) ) {
1780: T0 = car(T); C = T0[0]; Op = dp_dtop(T0[1],PF_DV);
1781: for ( K = 0; K < Dim; K++ )
1782: if ( B[K] == Op ) break;
1783: if ( K == Dim )
1784: error("afo");
1785: else
1786: A[J][K] = C;
1787: }
1788: }
1789: if ( Print ) print("");
1790: A = map(ctord,A);
1791: Lcm = [];
1792: for ( K = 0; K < Dim; K++ )
1793: for ( L = 0; L < Dim; L++ )
1794: Lcm = lcmf(Lcm,A[K][L][1]);
1795: for ( K = 0; K < Dim; K++ )
1796: for ( L = 0; L < Dim; L++ ) {
1797: Num = A[K][L][0]; Den = A[K][L][1];
1798: A[K][L] = Num*ftorat(divf(Lcm,Den));
1799: }
1800: Lcm = ftorat(Lcm);
1801: if ( !Lcm ) Lcm = 1;
1802: for ( K = 0; K < N; K++ ) {
1803: A = subst(A,V[K],Beta[K]*XV);
1804: Lcm = subst(Lcm,V[K],Beta[K]*XV);
1805: }
1806: A = map(red,A/Lcm);
1807: Lcm = 1;
1808: for ( K = 0; K < Dim; K++ )
1809: for ( L = 0; L < Dim; L++ )
1810: Lcm = lcm(dn(A[K][L]),Lcm);
1811: for ( K = 0; K < Dim; K++ )
1812: for ( L = 0; L < Dim; L++ )
1813: A[K][L] = nm(A[K][L])*sdiv(Lcm,dn(A[K][L]));
1814: Mat[I] = [A,Lcm];
1815: }
1816: Lcm = 1;
1817: for ( I = 0; I < N; I++ )
1818: Lcm = lcm(Mat[I][1],Lcm);
1819: R = matrix(Dim,Dim);
1820: for ( I = 0; I < N; I++ ) {
1821: A = Mat[I][0];
1822: for ( K = 0; K < Dim; K++ )
1823: for ( L = 0; L < Dim; L++ )
1824: R[K][L] += Beta[I]*nm(A[K][L])*sdiv(Lcm,Mat[I][1]);
1825: }
1826: Deg = 0;
1827: for ( K = 0; K < Dim; K++ )
1828: for ( L = 0; L < Dim; L++ ) {
1829: DegT = deg(R[K][L],XV);
1830: if ( DegT > Deg ) Deg = DegT;
1831: }
1832: RA = vector(Deg+1);
1833: for ( I = 0; I <= Deg; I++ )
1834: RA[I] = map(coef,R,I);
1835: return [[RA,Lcm],B];
1836: }
1837:
1838: def evalpf(F,V,X)
1839: {
1840: if ( !F ) return 0;
1841: F0 = F;
1842: N = length(V);
1843: for ( I = 0; I < N; I++ )
1844: F0 = subst(F0,V[I],X[I]);
1845: return F0;
1846: }
1847:
1848: def eval_pfaffian(Mat,Beta,SBeta,MN2,V,XI,F)
1849: {
1850: N = length(V);
1851: Dim = size(Mat[0][0])[0];
1852: R = vector(Dim);
1853: Mul = vector(N);
1854: X = XI*Beta;
1855: X = vtol(X);
1856: for ( K = 0; K < N; K++ )
1857: Mul[K] = Beta[K]/evalpf(Mat[K][1],V,X);
1858: for ( K = 0; K < N; K++ ) {
1859: MatK = Mat[K][0];
1860: for ( I = 0; I < Dim; I++ ) {
1861: for ( U = J = 0; J < Dim; J++ )
1862: U += substr2np(MatK[I][J],V,X)*F[J];
1863: R[I] += Mul[K]*U;
1864: }
1865: }
1866: for ( I = 0; I < Dim; I++ ) R[I] -= (SBeta-MN2/XI)*F[I];
1867: return R;
1868: }
1869:
1870: def eval_pfaffian2(Mat,Beta,SBeta,MN2,V,T,XI,F)
1871: {
1872: N = length(V);
1873: MA = Mat[0]; Den = subst(Mat[1],T,XI);
1874: Len = length(MA);
1875: Dim = size(MA[0])[0];
1876: R = vector(Dim);
1877: XII = 1;
1878: for ( I = 0; I < Len; I++ ) {
1879: R += MA[I]*(XII*F);
1880: XII *= XI;
1881: }
1882: R = R/Den - (SBeta-MN2/XI)*F;
1883: return R;
1884: }
1885:
1886:
1887: def genrk45(Z,M,N,V,PL,SL,Beta,SBeta,X,RK5,Step)
1888: {
1889: if ( type(Double=getopt(double)) == -1 ) {
1890: ctrl("bigfloat",1);
1891: Double = 0;
1892: } else
1893: ctrl("bigfloat",0);
1894: if ( type(Init=getopt(init)) == -1 )
1895: Init = 1;
1896: Len = length(V);
1897: DV = vector(Len);
1898: for ( I = 0; I < Len; I++ )
1899: DV[I] = strtov("d"+rtostr(V[I]));
1900: DV = vtol(DV);
1901: A = (1+M)/2; C = (1+M+N)/2;
1902: Z0 = subst(Z,a,A,c,C);
1903: T0 = time();
1904: Q = pfaffian2(Z0,V,Beta,t); Mat = Q[0]; Base = Q[1];
1905: MatLen = length(Q[0][0]);
1906: for ( I = 0; I < MatLen; I++ )
1907: Mat[0][I] *= eval(exp(0));
1908: T1 = time(); Tp = (T1[0]-T0[0])+(T1[1]-T0[1]);
1909: T0 = time();
1910: Beta = ltov(Beta)*eval(exp(0));
1911: X *= eval(exp(0));
1912: X1 = Beta*X;
1913: X0 = 0;
1914: for ( I = 0; I < Len; I++ )
1915: if ( Beta[I] > X0 ) X0 = Beta[I];
1916: X0 /= Init;
1917: /* Max Beta[I] is set to Init */
1918: Beta0 = Beta/X0;
1919: X0 = 1/X0;
1920: /* now Beta0 = X0*Beta */
1921: #if 0
1922: Prec = setbprec();
1923: setbprec(2*Prec);
1924: #endif
1925: L = psvalue(Z0,V,Beta0); PS = L[0];
1926: T1 = time(); Tps = (T1[0]-T0[0])+(T1[1]-T0[1]);
1927: T0 = time();
1928: #if 0
1929: setbprec(Prec);
1930: #endif
1931: Dim = length(Base);
1932: MN2 = M*N/2;
1933: ExpF = eval(X0^(MN2)*exp(-SBeta*X0));
1934: F = vector(Dim);
1935: for ( I = 0; I < Dim; I++ ) {
1936: DPS = act_op(dp_ptod(Base[I],DV),V,PS);
1937: for ( J = 0; J < Len; J++ )
1938: DPS = subst(DPS,V[J],Beta0[J]);
1939: F[I] = DPS*ExpF;
1940: }
1941: F0 = F*1;
1942: while ( 1 ) {
1943: F = F0*1;
1944: H = eval((X-X0)/Step);
1945: if ( Double ) {
1946: ctrl("bigfloat",0);
1947: X0 = todouble(X0);
1948: H = todouble(H);
1949: Beta = map(todouble,Beta);
1950: SBeta = todouble(SBeta);
1951: F = map(todouble,F);
1952: }
1953: R = [];
1954: GM = gammam(M,(M+1)/2)/gammam(M,(N+M+1)/2);
1955: O = eval(exp(0));
1956: if ( RK5 ) {
1957: A2 = 1/5*O; B21 = 1/5*O;
1958: A3 = 3/10*O;B31 = 3/40*O; B32 = 9/40*O;
1959: A4 = 3/5*O; B41 = 3/10*O; B42 = -9/10*O; B43 = 6/5*O;
1960: A5 = 1*O; B51 = -11/54*O; B52 = 5/2*O; B53 = -70/27*O; B54 = 35/27*O;
1961: A6 = 7/8*O; B61 = 1631/55296*O; B62 = 175/512*O; B63 = 575/13824*O; B64 = 44275/110592*O; B65 = 253/4096*O;
1962: C1 = 37/378*O; C2 = 0*O; C3 = 250/621*O; C4 = 125/594*O; C5 = 0*O; C6 = 512/1771*O;
1963: D1 = 2825/27648*O; D2 = 0*O; D3 = 18575/48384*O; D4 = 13525/55296*O; D5 = 277/14336*O; D6 = 1/4*O;
1964: for ( I = 0; I < Step; I++ ) {
1965: if ( Print ) {
1966: if ( !(I%1000) ) print([I],2);
1967: if ( !(I%10000) ) print("");
1968: }
1969: XI = X0+I*H;
1970: K1 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI,F);
1971: K2 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A2*H,F+B21*K1);
1972: K3 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A3*H,F+B31*K1+B32*K2);
1973: K4 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A4*H,F+B41*K1+B42*K2+B43*K3);
1974: K5 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A5*H,F+B51*K1+B52*K2+B53*K3+B54*K4);
1975: K6 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+A6*H,F+B61*K1+B62*K2+B63*K3+B64*K4+B65*K5);
1976: F += C1*K1+C2*K2+C3*K3+C4*K4+C5*K5+C6*K6;
1977: T = GM*F[0]*genprc(M,N,PL,SL,XI+H);
1978: if ( T < 0 || T > 1 ) break;
1979: R = cons([XI+H,T],R);
1980: }
1981: } else {
1982: for ( I = 0; I < Step; I++ ) {
1983: if ( Print ) {
1984: if ( !(I%1000) ) print([I],2);
1985: if ( !(I%10000) ) print("");
1986: }
1987: XI = X0+I*H;
1988: K1 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI,F);
1989: K2 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+H/2,F+1/2*K1);
1990: K3 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+H/2,F+1/2*K2);
1991: K4 = H*eval_pfaffian2(Mat,Beta,SBeta,MN2,V,t,XI+H,F+K3);
1992: F += 1/6*K1+1/3*K2+1/3*K3+1/6*K4;
1993: T = GM*F[0]*genprc(M,N,PL,SL,XI+H);
1994: if ( T < 0 || T > 1 ) break;
1995: R = cons([XI+H,T],R);
1996: }
1997: }
1998: if ( Print ) print("");
1999: if ( I == Step ) {
2000: T1 = time(); Trk = (T1[0]-T0[0])+(T1[1]-T0[1]);
2001: return reverse(R);
2002: } else {
2003: Step *= 2;
2004: if ( Print ) print("Step=",0); if ( Print ) print(Step);
2005: }
2006: }
2007: }
2008:
2009: def solve_leq(L,V)
2010: {
2011: N = length(V);
2012: B = [];
2013: for ( I = 0; I < N; I++, L = cdr(L) ) B = cons(car(L),B);
2014: while ( 1 ) {
2015: Sol = nd_f4(B,V,0,0);
2016: if ( length(Sol) != N ) {
2017: B = cons(car(L),Sol); L = cdr(L);
2018: } else {
2019: return Sol;
2020: }
2021: }
2022: }
2023:
2024:
2025: def ptom(P,V)
2026: {
2027: M = length(P); N = length(V);
2028: A = matrix(M,N+1);
2029: for ( I = 0; I < M; I++ ) {
2030: F = ptozp(P[I]);
2031: VT = V;
2032: J = 0;
2033: while ( F )
2034: if ( type(F) <= 1 ) {
2035: A[I][N] = F; break;
2036: } else {
2037: for ( FV = var(F); FV != car(VT); VT = cdr(VT), J++ );
2038: A[I][J] = coef(F,1);
2039: F -= A[I][J]*FV;
2040: }
2041: }
2042: return A;
2043: }
2044:
2045: def vton(V1,V)
2046: {
2047: N = length(V);
2048: for ( I = 0; I < N; I++ )
2049: if ( V1 == V[I] ) return I;
2050: error("vton : no such variable");
2051: }
2052:
2053: def ftotex(F,V)
2054: {
2055: if ( F == [] ) return "1";
2056: R = "";
2057: for ( T = F; T != []; T = cdr(T) ) {
2058: Fac = car(T)[0]; M = car(T)[1];
2059: V1 = var(Fac); NV1 = vton(V1,V);
2060: if ( Fac == V1 ) SFac = "y_{"+rtostr(NV1)+"}";
2061: else {
2062: V2 = var(Fac-V1);
2063: NV2 = vton(V2,V);
2064: SFac = "(y_{"+rtostr(NV1)+"}-y_{"+rtostr(NV2)+"})";
2065: }
2066: if ( M == 1 ) R += SFac;
2067: else R += SFac+"^{"+rtostr(M)+"}";
2068: }
2069: return R;
2070: }
2071:
2072: def ctotex(C,V)
2073: {
2074: R = "";
2075: for ( T = C; T != []; T = cdr(T) ) {
2076: if ( T != C ) R += "+";
2077: N = quotetotex(objtoquote(car(T)[0]));
2078: if ( car(T)[1] == [] ) {
2079: R += "("+N+")";
2080: } else {
2081: D = ftotex(car(T)[1],V);
2082: R += "\\frac{"+N+"}{"+D+"}";
2083: }
2084: }
2085: return R;
2086: }
2087:
2088: def optotex(Op,DV)
2089: {
2090: E = dp_etov(Op);
2091: N = length(E);
2092: R = "";
2093: for ( I = 0; I < N; I++ )
2094: if ( E[I] )
2095: if ( E[I] == 1 )
2096: R += "\\partial_{"+rtostr(I)+"}";
2097: else
2098: R += "\\partial_{"+rtostr(I)+"}^{"+rtostr(E[I])+"}";
2099: return R;
2100: }
2101:
2102: def pftotex(P,V,DV)
2103: {
2104: R = "";
2105: for ( T = P; T != []; T = cdr(T) ) {
2106: if ( T != P ) R += "+";
2107: C = ctotex(car(T)[0],V); Op = optotex(car(T)[1],DV);
2108: R += "("+C+")"+Op+"\\\\";
2109: }
2110: return R;
2111: }
2112:
2113: def eqtotex(Z)
2114: {
2115: shell("rm -f __tmp__.tex");
2116: output("__tmp__.tex");
2117: N = length(Z);
2118: print("\\documentclass[12pt]{jsarticle}");
2119: print("\\begin{document}");
2120: print("\\large\\noindent");
2121: for ( I = 0; I < N; I++ ) {
2122: print("$h_"+rtostr(I)+"=",0);
2123: T = mulcpf([[-1,[]]],Z[I]);
2124: print(pftotex(T,PF_V,PF_DV),0);
2125: print("$\\\\");
2126: }
2127: print("\\end{document}");
2128: output();
2129: shell("latex __tmp__");
2130: shell("xdvi __tmp__");
2131: }
2132:
2133: /*
2134: * for a partition L=[N1,...,Nl]
2135: * compute the coefficient of x1^N1*...*xl^Nl in phi((x1+...+xM)^N)
2136: * where phi(x(s(1))^n1*...*x(s(k))^nk)=x1^n1*...*xk^nk (n1>=...>=nk)
2137: * if count_dup(L)=[I1,...,Is], then the coefficient is
2138: * C=N!/(N1!*...*Nl!)*M!/((M-l)!*I1!*...*Is!)
2139: * return C*<<0,...0,N1,...,Nl,0,...,0>>
2140: position 0 Base
2141: */
2142:
2143: def coef_partition(L,N,M,Base)
2144: {
2145: R = fac(N);
2146: for ( T = L; T != []; T = cdr(T) )
2147: R = sdiv(R,fac(car(T)));
2148: S = count_dup(L);
2149: R *= fac(M)/fac(M-length(L));
2150: for ( T = S; T != []; T = cdr(T) )
2151: R = sdiv(R,fac(car(T)));
2152: E = vector(length(PF_DV));
2153: for ( I = Base, T = L; T != []; T = cdr(T), I++ )
2154: E[I] = car(T);
2155: return R*dp_vtoe(E);
2156: }
2157:
2158: /*
2159: * for a partition L=[N1,...,Nk], return [I1,...,Is]
2160: * s.t. N1=...=N(I1)>N(I1+1)=...=N(I1+I2)>...
2161: */
2162:
2163: def count_dup(L)
2164: {
2165: D = [];
2166: T = L;
2167: while ( T != [] ) {
2168: T0 = car(T); T = cdr(T);
2169: for ( I = 1; T != [] && car(T) == T0; T = cdr(T) ) I++;
2170: D = cons(I,D);
2171: }
2172: return D;
2173: }
2174:
2175: /* generate (N_1,...,N_L) s.t. N_1+...+N_L=N, N_1>=...>=N_L>=1, L<= K */
2176:
2177: def all_partition(N,K)
2178: {
2179: R = [];
2180: for ( I = 1; I <= K; I++ )
2181: R = append(partition(N,I),R);
2182: return map(reverse,R);
2183: }
2184:
2185: /* generate (N_1,...,N_K) s.t. N_1+...+N_K=N, 1<=N_1<=...<=N_K */
2186:
2187: def partition(N,K)
2188: {
2189: if ( N < K ) return [];
2190: else if ( N == K ) {
2191: S = [];
2192: for ( I = 0; I < K; I++ )
2193: S = cons(1,S);
2194: return [S];
2195: } else if ( K == 0 )
2196: return [];
2197: else {
2198: R = [];
2199: L = partition(N-K,K);
2200: for ( T = L; T != []; T = cdr(T) ) {
2201: S = [];
2202: for ( U = car(T); U != []; U = cdr(U) )
2203: S = cons(car(U)+1,S);
2204: R = cons(reverse(S),R);
2205: }
2206: L = partition(N-1,K-1);
2207: for ( T = L; T != []; T = cdr(T) )
2208: R = cons(cons(1,car(T)),R);
2209: return R;
2210: }
2211: }
2212:
2213: def mul_lopitalpf2(P,I,N,K,I0,I1,E)
2214: {
2215: YI = PF_V[I]; DYI = PF_DV[I];
2216: R = [mtot(DYI^2,1)];
2217: for ( J = I0; J <= I1; J++ ) {
2218: if ( J == I ) continue;
2219: YJ = PF_V[J]; DYJ = PF_DV[J];
2220: R = addpf([mtot(1/2*DYI,YI-YJ)],R);
2221: R = addpf([mtot(-1/2*DYJ,YI-YJ)],R);
2222: }
2223: S = subpf(P[I],R);
2224: R = loph(E,I,I0,I1);
2225: if ( type(getopt(high)) == -1 )
2226: S = mul_lopitalpf(S,N,K,I0,I1,E)[1];
2227: else
2228: S = 0;
2229: return [R,S];
2230: }
2231:
2232: /* the highest part of lim_{Y(I+1),...,Y(I+N-1)->YI} E(PI) */
2233: /* where E1 = E+(0,...,2,...,0) */
2234:
2235: def loph(E,I,I0,I1)
2236: {
2237: E1 = E*1;
2238: E1[I] += 2;
2239: R = [[[[1,[]]],dp_vtoe(E1)]];
2240: S = lopn(E,I,I0,I1);
2241: S = mulcpf([[1/2,[]]],S);
2242: R = addpf(R,S);
2243: R = simppf(R,I0,I1-I0+1);
2244: return R;
2245: }
2246:
2247: /* lim_{Y(I+1),...,Y(I+N-1)->YI} dy^E(1/(YI-Y(I+1))+...+1/(YI-Y(I+N-1))) */
2248:
2249: def lopn(E,I,I0,I1)
2250: {
2251: R = [];
2252: for ( K = I0; K <= I1; K++ ) {
2253: if ( K != I ) {
2254: L = lop1(E,I,K);
2255: R = addpf(R,L);
2256: }
2257: }
2258: return R;
2259: }
2260:
2261: /* lim_{YJ->YI} dy^E(1/(YI-YJ) */
2262: def lop1(E,I,J)
2263: {
2264: YI = PF_V[I]; YJ = PF_V[J];
2265: DYI = PF_DV[I]; DYJ = PF_DV[J];
2266: R = [];
2267: R = addpf([mtot(DYI,YI-YJ)],R);
2268: R = addpf([mtot(-DYJ,YI-YJ)],R);
2269: N = length(PF_V);
2270: BI = E[I]; BJ = E[J];
2271: for ( K = 1, D = 1; K < N; K++ )
2272: if ( K != I && K != J )
2273: D *= PF_DV[K]^E[K];
2274: S = [mtot(-1/(BJ+1)*DYI^(BI+1)*DYJ^(BJ+1)*D,1)];
2275: for ( K = 0; K <= BI; K++ )
2276: if ( -BI+BJ+2*K+2 )
2277: S = addpf([mtot(fac(BI)*fac(BJ)*(-BI+BJ+2*K+2)/fac(BI-K)/fac(BJ+K+2)*DYI^(BI-K)*DYJ^(BJ+K+2)*D,1)],S);
2278: return S;
2279: }
2280:
2281: def getchi(N)
2282: {
2283: CHI=[
2284: [5 ,11.07049769,1.145476226],
2285: [10 ,18.30703805,3.940299136],
2286: [20 ,31.41043284,10.85081139],
2287: [50 ,67.50480655,34.76425168],
2288: [100,124.3421134,77.92946517],
2289: [500,553.1268089,449.1467564],
2290: [1000, 1074.679449,927.594363]];
2291: for ( T = CHI; T != []; T = cdr(T) )
2292: if ( car(T)[0] == N ) return car(T);
2293: return [];
2294: }
2295:
2296: def message(S) {
2297: dp_gr_print(S);
2298: Print = S;
2299: }
2300: endmodule$
2301: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to wishartd.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / testing / noro