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